ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents):
Revision 1726 by tim, Wed Nov 10 22:50:03 2004 UTC vs.
Revision 1886 by tim, Wed Dec 15 16:13:35 2004 UTC

# Line 31 | Line 31
31   */
32  
33   #include <algorithm>
34 <
34 > #include <set>
35   #include "brains/SimInfo.hpp"
36 + #include "primitives/Molecule.hpp"
37 + #include "UseTheForce/doForces_interface.h"
38 + #include "UseTheForce/notifyCutoffs_interface.h"
39   #include "utils/MemoryUtils.hpp"
40 + #include "utils/simError.h"
41  
42 + #ifdef IS_MPI
43 + #include "UseTheForce/mpiComponentPlan.h"
44 + #include "UseTheForce/DarkSide/simParallel_interface.h"
45 + #endif
46 +
47   namespace oopse {
48  
49 < SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
50 <        nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){
49 > SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
50 >                                ForceField* ff, Globals* simParams) :
51 >                                forceField_(ff), simParams_(simParams),
52 >                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
53 >                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
54 >                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
55 >                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
56 >                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
57 >                                sman_(NULL), fortranInitialized_(false) {
58 >
59 >            
60 >    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
61 >    MoleculeStamp* molStamp;
62 >    int nMolWithSameStamp;
63 >    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
64 >    int nGroups = 0;          //total cutoff groups defined in meta-data file
65 >    CutoffGroupStamp* cgStamp;    
66 >    RigidBodyStamp* rbStamp;
67 >    int nRigidAtoms = 0;
68 >    
69 >    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
70 >        molStamp = i->first;
71 >        nMolWithSameStamp = i->second;
72 >        
73 >        addMoleculeStamp(molStamp, nMolWithSameStamp);
74 >
75 >        //calculate atoms in molecules
76 >        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
77 >
78 >
79 >        //calculate atoms in cutoff groups
80 >        int nAtomsInGroups = 0;
81 >        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
82 >        
83 >        for (int j=0; j < nCutoffGroupsInStamp; j++) {
84 >            cgStamp = molStamp->getCutoffGroup(j);
85 >            nAtomsInGroups += cgStamp->getNMembers();
86 >        }
87 >
88 >        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
89 >        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
90 >
91 >        //calculate atoms in rigid bodies
92 >        int nAtomsInRigidBodies = 0;
93 >        int nRigidBodiesInStamp = molStamp->getNCutoffGroups();
94 >        
95 >        for (int j=0; j < nRigidBodiesInStamp; j++) {
96 >            rbStamp = molStamp->getRigidBody(j);
97 >            nAtomsInRigidBodies += rbStamp->getNMembers();
98 >        }
99 >
100 >        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
101 >        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
102 >        
103 >    }
104 >
105 >    //every free atom (atom does not belong to cutoff groups) is a cutoff group
106 >    //therefore the total number of cutoff groups in the system is equal to
107 >    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
108 >    //file plus the number of cutoff groups defined in meta-data file
109 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
110 >
111 >    //every free atom (atom does not belong to rigid bodies) is an integrable object
112 >    //therefore the total number of  integrable objects in the system is equal to
113 >    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
114 >    //file plus the number of  rigid bodies defined in meta-data file
115 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
116  
117 +    nGlobalMols_ = molStampIds_.size();
118 +
119 + #ifdef IS_MPI    
120 +    molToProcMap_.resize(nGlobalMols_);
121 + #endif
122 +    
123   }
124  
125   SimInfo::~SimInfo() {
126      //MemoryUtils::deleteVectorOfPointer(molecules_);
127 +
128 +    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
129 +    
130      delete sman_;
131 +    delete simParams_;
132 +    delete forceField_;
133  
134   }
135  
136  
137   bool SimInfo::addMolecule(Molecule* mol) {
138      MoleculeIterator i;
54 –    i = std::find(molecules_.begin(), molecules_.end(), mol);
55 –    if (i != molecules_.end() ) {
139  
140 <        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
140 >    i = molecules_.find(mol->getGlobalIndex());
141 >    if (i == molecules_.end() ) {
142 >
143 >        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
144          
145          nAtoms_ += mol->getNAtoms();
146          nBonds_ += mol->getNBonds();
# Line 65 | Line 151 | bool SimInfo::addMolecule(Molecule* mol) {
151          nCutoffGroups_ += mol->getNCutoffGroups();
152          nConstraints_ += mol->getNConstraints();
153  
154 +        addExcludePairs(mol);
155 +        
156          return true;
157      } else {
158          return false;
# Line 73 | Line 161 | bool SimInfo::removeMolecule(Molecule* mol) {
161  
162   bool SimInfo::removeMolecule(Molecule* mol) {
163      MoleculeIterator i;
164 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
164 >    i = molecules_.find(mol->getGlobalIndex());
165  
166      if (i != molecules_.end() ) {
167  
168 +        assert(mol == i->second);
169 +        
170          nAtoms_ -= mol->getNAtoms();
171          nBonds_ -= mol->getNBonds();
172          nBends_ -= mol->getNBends();
# Line 86 | Line 176 | bool SimInfo::removeMolecule(Molecule* mol) {
176          nCutoffGroups_ -= mol->getNCutoffGroups();
177          nConstraints_ -= mol->getNConstraints();
178  
179 +        removeExcludePairs(mol);
180          molecules_.erase(mol->getGlobalIndex());
181  
182          delete mol;
# Line 101 | Line 192 | Molecule* SimInfo::beginMolecule(MoleculeIterator& i)
192          
193   Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
194      i = molecules_.begin();
195 <    return i == molecules_.end() ? NULL : *i;
195 >    return i == molecules_.end() ? NULL : i->second;
196   }    
197  
198   Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
199      ++i;
200 <    return i == molecules_.end() ? NULL : *i;    
200 >    return i == molecules_.end() ? NULL : i->second;    
201   }
202  
203  
# Line 145 | Line 236 | void SimInfo::calcNdf() {
236      ndf_ = ndf_local;
237   #endif
238  
239 <    // nZconstraints is global, as are the 3 COM translations for the
239 >    // nZconstraints_ is global, as are the 3 COM translations for the
240      // entire system:
241 <    ndf_ = ndf_ - 3 - nZconstraints;
241 >    ndf_ = ndf_ - 3 - nZconstraint_;
242  
243   }
244  
# Line 198 | Line 289 | void SimInfo::calcNdfTrans() {
289      ndfTrans_ = ndfTrans_local;
290   #endif
291  
292 <    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
292 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
293  
294   }
295  
# Line 230 | Line 321 | void SimInfo::addExcludePairs(Molecule* mol) {
321          exclude_.addPair(b, c);        
322      }
323  
324 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
324 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
325          a = torsion->getAtomA()->getGlobalIndex();
326          b = torsion->getAtomB()->getGlobalIndex();        
327          c = torsion->getAtomC()->getGlobalIndex();        
# Line 275 | Line 366 | void SimInfo::removeExcludePairs(Molecule* mol) {
366          exclude_.removePair(b, c);        
367      }
368  
369 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
369 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
370          a = torsion->getAtomA()->getGlobalIndex();
371          b = torsion->getAtomB()->getGlobalIndex();        
372          c = torsion->getAtomC()->getGlobalIndex();        
# Line 299 | Line 390 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
390      curStampId = molStampIds_.size();
391  
392      moleculeStamps_.push_back(molStamp);
393 <    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
393 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
394   }
395  
396 < std::ostream& operator <<(ostream& o, SimInfo& info) {
396 > void SimInfo::update() {
397 >
398 >    setupSimType();
399 >
400 > #ifdef IS_MPI
401 >    setupFortranParallel();
402 > #endif
403 >
404 >    setupFortranSim();
405 >
406 >    //setup fortran force field
407 >    /** @deprecate */    
408 >    int isError = 0;
409 >    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
410 >    if(isError){
411 >        sprintf( painCave.errMsg,
412 >         "ForceField error: There was an error initializing the forceField in fortran.\n" );
413 >        painCave.isFatal = 1;
414 >        simError();
415 >    }
416 >  
417 >    
418 >    setupCutoff();
419 >
420 >    calcNdf();
421 >    calcNdfRaw();
422 >    calcNdfTrans();
423 >
424 >    fortranInitialized_ = true;
425 > }
426 >
427 > std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
428 >    SimInfo::MoleculeIterator mi;
429 >    Molecule* mol;
430 >    Molecule::AtomIterator ai;
431 >    Atom* atom;
432 >    std::set<AtomType*> atomTypes;
433 >
434 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
435 >
436 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
437 >            atomTypes.insert(atom->getAtomType());
438 >        }
439 >        
440 >    }
441 >
442 >    return atomTypes;        
443 > }
444 >
445 > void SimInfo::setupSimType() {
446 >    std::set<AtomType*>::iterator i;
447 >    std::set<AtomType*> atomTypes;
448 >    atomTypes = getUniqueAtomTypes();
449 >    
450 >    int useLennardJones = 0;
451 >    int useElectrostatic = 0;
452 >    int useEAM = 0;
453 >    int useCharge = 0;
454 >    int useDirectional = 0;
455 >    int useDipole = 0;
456 >    int useGayBerne = 0;
457 >    int useSticky = 0;
458 >    int useShape = 0;
459 >    int useFLARB = 0; //it is not in AtomType yet
460 >    int useDirectionalAtom = 0;    
461 >    int useElectrostatics = 0;
462 >    //usePBC and useRF are from simParams
463 >    int usePBC = simParams_->getPBC();
464 >    int useRF = simParams_->getUseRF();
465 >
466 >    //loop over all of the atom types
467 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
468 >        useLennardJones |= (*i)->isLennardJones();
469 >        useElectrostatic |= (*i)->isElectrostatic();
470 >        useEAM |= (*i)->isEAM();
471 >        useCharge |= (*i)->isCharge();
472 >        useDirectional |= (*i)->isDirectional();
473 >        useDipole |= (*i)->isDipole();
474 >        useGayBerne |= (*i)->isGayBerne();
475 >        useSticky |= (*i)->isSticky();
476 >        useShape |= (*i)->isShape();
477 >    }
478 >
479 >    if (useSticky || useDipole || useGayBerne || useShape) {
480 >        useDirectionalAtom = 1;
481 >    }
482 >
483 >    if (useCharge || useDipole) {
484 >        useElectrostatics = 1;
485 >    }
486 >
487 > #ifdef IS_MPI    
488 >    int temp;
489 >
490 >    temp = usePBC;
491 >    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
492 >
493 >    temp = useDirectionalAtom;
494 >    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
495 >
496 >    temp = useLennardJones;
497 >    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
498 >
499 >    temp = useElectrostatics;
500 >    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
501 >
502 >    temp = useCharge;
503 >    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
504 >
505 >    temp = useDipole;
506 >    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
507 >
508 >    temp = useSticky;
509 >    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
510 >
511 >    temp = useGayBerne;
512 >    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
513 >
514 >    temp = useEAM;
515 >    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
516 >
517 >    temp = useShape;
518 >    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
519 >
520 >    temp = useFLARB;
521 >    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
522 >
523 >    temp = useRF;
524 >    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
525 >    
526 > #endif
527 >
528 >    fInfo_.SIM_uses_PBC = usePBC;    
529 >    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
530 >    fInfo_.SIM_uses_LennardJones = useLennardJones;
531 >    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
532 >    fInfo_.SIM_uses_Charges = useCharge;
533 >    fInfo_.SIM_uses_Dipoles = useDipole;
534 >    fInfo_.SIM_uses_Sticky = useSticky;
535 >    fInfo_.SIM_uses_GayBerne = useGayBerne;
536 >    fInfo_.SIM_uses_EAM = useEAM;
537 >    fInfo_.SIM_uses_Shapes = useShape;
538 >    fInfo_.SIM_uses_FLARB = useFLARB;
539 >    fInfo_.SIM_uses_RF = useRF;
540 >
541 >    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
542 >
543 >        if (simParams_->haveDielectric()) {
544 >            fInfo_.dielect = simParams_->getDielectric();
545 >        } else {
546 >            sprintf(painCave.errMsg,
547 >                    "SimSetup Error: No Dielectric constant was set.\n"
548 >                    "\tYou are trying to use Reaction Field without"
549 >                    "\tsetting a dielectric constant!\n");
550 >            painCave.isFatal = 1;
551 >            simError();
552 >        }
553 >        
554 >    } else {
555 >        fInfo_.dielect = 0.0;
556 >    }
557 >
558 > }
559 >
560 > void SimInfo::setupFortranSim() {
561 >    int isError;
562 >    int nExclude;
563 >    std::vector<int> fortranGlobalGroupMembership;
564 >    
565 >    nExclude = exclude_.getSize();
566 >    isError = 0;
567 >
568 >    //globalGroupMembership_ is filled by SimCreator    
569 >    for (int i = 0; i < nGlobalAtoms_; i++) {
570 >        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
571 >    }
572  
573 +    //calculate mass ratio of cutoff group
574 +    std::vector<double> mfact;
575 +    SimInfo::MoleculeIterator mi;
576 +    Molecule* mol;
577 +    Molecule::CutoffGroupIterator ci;
578 +    CutoffGroup* cg;
579 +    Molecule::AtomIterator ai;
580 +    Atom* atom;
581 +    double totalMass;
582 +
583 +    //to avoid memory reallocation, reserve enough space for mfact
584 +    mfact.reserve(getNCutoffGroups());
585 +    
586 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
587 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
588 +
589 +            totalMass = cg->getMass();
590 +            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
591 +                        mfact.push_back(atom->getMass()/totalMass);
592 +            }
593 +
594 +        }      
595 +    }
596 +
597 +    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
598 +    std::vector<int> identArray;
599 +
600 +    //to avoid memory reallocation, reserve enough space identArray
601 +    identArray.reserve(getNAtoms());
602 +    
603 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
604 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
605 +            identArray.push_back(atom->getIdent());
606 +        }
607 +    }    
608 +
609 +    //fill molMembershipArray
610 +    //molMembershipArray is filled by SimCreator    
611 +    std::vector<int> molMembershipArray(nGlobalAtoms_);
612 +    for (int i = 0; i < nGlobalAtoms_; i++) {
613 +        molMembershipArray[i] = globalMolMembership_[i] + 1;
614 +    }
615 +    
616 +    //setup fortran simulation
617 +    //gloalExcludes and molMembershipArray should go away (They are never used)
618 +    //why the hell fortran need to know molecule?
619 +    //OOPSE = Object-Obfuscated Parallel Simulation Engine
620 +    int nGlobalExcludes = 0;
621 +    int* globalExcludes = NULL;
622 +    int* excludeList = exclude_.getExcludeList();
623 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
624 +                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
625 +                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
626 +
627 +    if( isError ){
628 +
629 +        sprintf( painCave.errMsg,
630 +                 "There was an error setting the simulation information in fortran.\n" );
631 +        painCave.isFatal = 1;
632 +        painCave.severity = OOPSE_ERROR;
633 +        simError();
634 +    }
635 +
636 + #ifdef IS_MPI
637 +    sprintf( checkPointMsg,
638 +       "succesfully sent the simulation information to fortran.\n");
639 +    MPIcheckPoint();
640 + #endif // is_mpi
641 + }
642 +
643 +
644 + #ifdef IS_MPI
645 + void SimInfo::setupFortranParallel() {
646 +    
647 +    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
648 +    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
649 +    std::vector<int> localToGlobalCutoffGroupIndex;
650 +    SimInfo::MoleculeIterator mi;
651 +    Molecule::AtomIterator ai;
652 +    Molecule::CutoffGroupIterator ci;
653 +    Molecule* mol;
654 +    Atom* atom;
655 +    CutoffGroup* cg;
656 +    mpiSimData parallelData;
657 +    int isError;
658 +
659 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
660 +
661 +        //local index(index in DataStorge) of atom is important
662 +        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
663 +            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
664 +        }
665 +
666 +        //local index of cutoff group is trivial, it only depends on the order of travesing
667 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
668 +            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
669 +        }        
670 +        
671 +    }
672 +
673 +    //fill up mpiSimData struct
674 +    parallelData.nMolGlobal = getNGlobalMolecules();
675 +    parallelData.nMolLocal = getNMolecules();
676 +    parallelData.nAtomsGlobal = getNGlobalAtoms();
677 +    parallelData.nAtomsLocal = getNAtoms();
678 +    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
679 +    parallelData.nGroupsLocal = getNCutoffGroups();
680 +    parallelData.myNode = worldRank;
681 +    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
682 +
683 +    //pass mpiSimData struct and index arrays to fortran
684 +    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
685 +                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
686 +                    &localToGlobalCutoffGroupIndex[0], &isError);
687 +
688 +    if (isError) {
689 +        sprintf(painCave.errMsg,
690 +                "mpiRefresh errror: fortran didn't like something we gave it.\n");
691 +        painCave.isFatal = 1;
692 +        simError();
693 +    }
694 +
695 +    sprintf(checkPointMsg, " mpiRefresh successful.\n");
696 +    MPIcheckPoint();
697 +
698 +
699 + }
700 +
701 + #endif
702 +
703 + double SimInfo::calcMaxCutoffRadius() {
704 +
705 +
706 +    std::set<AtomType*> atomTypes;
707 +    std::set<AtomType*>::iterator i;
708 +    std::vector<double> cutoffRadius;
709 +
710 +    //get the unique atom types
711 +    atomTypes = getUniqueAtomTypes();
712 +
713 +    //query the max cutoff radius among these atom types
714 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
715 +        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
716 +    }
717 +
718 +    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
719 + #ifdef IS_MPI
720 +    //pick the max cutoff radius among the processors
721 + #endif
722 +
723 +    return maxCutoffRadius;
724 + }
725 +
726 + void SimInfo::setupCutoff() {
727 +    double rcut_;  //cutoff radius
728 +    double rsw_; //switching radius
729 +    
730 +    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
731 +        
732 +        if (!simParams_->haveRcut()){
733 +            sprintf(painCave.errMsg,
734 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
735 +                "\tOOPSE will use a default value of 15.0 angstroms"
736 +                "\tfor the cutoffRadius.\n");
737 +            painCave.isFatal = 0;
738 +            simError();
739 +            rcut_ = 15.0;
740 +        } else{
741 +            rcut_ = simParams_->getRcut();
742 +        }
743 +
744 +        if (!simParams_->haveRsw()){
745 +            sprintf(painCave.errMsg,
746 +                "SimCreator Warning: No value was set for switchingRadius.\n"
747 +                "\tOOPSE will use a default value of\n"
748 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
749 +            painCave.isFatal = 0;
750 +            simError();
751 +            rsw_ = 0.95 * rcut_;
752 +        } else{
753 +            rsw_ = simParams_->getRsw();
754 +        }
755 +
756 +    } else {
757 +        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
758 +        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
759 +        
760 +        if (simParams_->haveRcut()) {
761 +            rcut_ = simParams_->getRcut();
762 +        } else {
763 +            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
764 +            rcut_ = calcMaxCutoffRadius();
765 +        }
766 +
767 +        if (simParams_->haveRsw()) {
768 +            rsw_  = simParams_->getRsw();
769 +        } else {
770 +            rsw_ = rcut_;
771 +        }
772 +    
773 +    }
774 +        
775 +    double rnblist = rcut_ + 1; // skin of neighbor list
776 +
777 +    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
778 +    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
779 + }
780 +
781 + void SimInfo::addProperty(GenericData* genData) {
782 +    properties_.addProperty(genData);  
783 + }
784 +
785 + void SimInfo::removeProperty(const std::string& propName) {
786 +    properties_.removeProperty(propName);  
787 + }
788 +
789 + void SimInfo::clearProperties() {
790 +    properties_.clearProperties();
791 + }
792 +
793 + std::vector<std::string> SimInfo::getPropertyNames() {
794 +    return properties_.getPropertyNames();  
795 + }
796 +      
797 + std::vector<GenericData*> SimInfo::getProperties() {
798 +    return properties_.getProperties();
799 + }
800 +
801 + GenericData* SimInfo::getPropertyByName(const std::string& propName) {
802 +    return properties_.getPropertyByName(propName);
803 + }
804 +
805 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
806 +    sman_ = sman;
807 +
808 +    Molecule* mol;
809 +    RigidBody* rb;
810 +    Atom* atom;
811 +    SimInfo::MoleculeIterator mi;
812 +    Molecule::RigidBodyIterator rbIter;
813 +    Molecule::AtomIterator atomIter;;
814 +
815 +    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
816 +        
817 +        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
818 +            atom->setSnapshotManager(sman_);
819 +        }
820 +        
821 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
822 +            rb->setSnapshotManager(sman_);
823 +        }
824 +    }    
825 +    
826 + }
827 +
828 + Vector3d SimInfo::getComVel(){
829 +    SimInfo::MoleculeIterator i;
830 +    Molecule* mol;
831 +
832 +    Vector3d comVel(0.0);
833 +    double totalMass = 0.0;
834 +    
835 +
836 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
837 +        double mass = mol->getMass();
838 +        totalMass += mass;
839 +        comVel += mass * mol->getComVel();
840 +    }  
841 +
842 +    comVel /= totalMass;
843 +
844 +    return comVel;
845 + }
846 +
847 + Vector3d SimInfo::getCom(){
848 +    SimInfo::MoleculeIterator i;
849 +    Molecule* mol;
850 +
851 +    Vector3d com(0.0);
852 +    double totalMass = 0.0;
853 +    
854 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
855 +        double mass = mol->getMass();
856 +        totalMass += mass;
857 +        com += mass * mol->getCom();
858 +    }  
859 +
860 +    com /= totalMass;
861 +
862 +    return com;
863 +
864 + }        
865 +
866 + std::ostream& operator <<(std::ostream& o, SimInfo& info) {
867 +
868      return o;
869   }
870  
871   }//end namespace oopse
872 +

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines