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 1719 by tim, Fri Nov 5 23:38:27 2004 UTC vs.
Revision 1883 by tim, Mon Dec 13 22:30:27 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_);
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 <    std::vector<Molecule*>::iterator i;
54 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
55 <    if (i != molecules_.end() ) {
56 <        molecules_.push_back(mol);
138 >    MoleculeIterator i;
139  
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();
147          nBends_ += mol->getNBends();
# Line 64 | 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 71 | Line 160 | bool SimInfo::removeMolecule(Molecule* mol) {
160   }
161  
162   bool SimInfo::removeMolecule(Molecule* mol) {
163 <    std::vector<Molecule*>::iterator i;
164 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
163 >    MoleculeIterator i;
164 >    i = molecules_.find(mol->getGlobalIndex());
165  
166      if (i != molecules_.end() ) {
167 <        molecules_.push_back(mol);
167 >
168 >        assert(mol == i->second);
169 >        
170          nAtoms_ -= mol->getNAtoms();
171          nBonds_ -= mol->getNBonds();
172          nBends_ -= mol->getNBends();
# Line 85 | 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;
183 +        
184          return true;
185      } else {
186          return false;
# Line 94 | Line 190 | Molecule* SimInfo::beginMolecule(std::vector<Molecule*
190   }    
191  
192          
193 < Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
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(std::vector<Molecule*>::iterator& i) {
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  
204 < void SimInfo::calcNDF() {
204 > void SimInfo::calcNdf() {
205      int ndf_local;
206 <    std::vector<Molecule*>::iterator i;
206 >    MoleculeIterator i;
207      std::vector<StuntDouble*>::iterator j;
208      Molecule* mol;
209      StuntDouble* integrableObject;
# Line 140 | 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  
245 < void SimInfo::calcNDFRaw() {
245 > void SimInfo::calcNdfRaw() {
246      int ndfRaw_local;
247  
248 <    std::vector<Molecule*>::iterator i;
248 >    MoleculeIterator i;
249      std::vector<StuntDouble*>::iterator j;
250      Molecule* mol;
251      StuntDouble* integrableObject;
# Line 181 | Line 277 | void SimInfo::calcNDFTrans() {
277   #endif
278   }
279  
280 < void SimInfo::calcNDFTrans() {
280 > void SimInfo::calcNdfTrans() {
281      int ndfTrans_local;
282  
283      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
# Line 193 | 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 225 | 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 270 | 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 282 | Line 378 | void SimInfo::removeExcludePairs(Molecule* mol) {
378          exclude_.removePair(b, c);
379          exclude_.removePair(b, d);
380          exclude_.removePair(c, d);        
381 +    }
382 +
383 + }
384 +
385 +
386 + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
387 +    int curStampId;
388 +
389 +    //index from 0
390 +    curStampId = molStampIds_.size();
391 +
392 +    moleculeStamps_.push_back(molStamp);
393 +    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
394 + }
395 +
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 +    bool usePBC = simParams_->getPBC();
464 +    bool 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