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

Diff Legend

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