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

Comparing branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp (file contents):
Revision 1725 by tim, Wed Nov 10 22:01:06 2004 UTC vs.
Revision 1895 by tim, Thu Dec 16 19:20:49 2004 UTC

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

Diff Legend

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