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 1842 by tim, Fri Dec 3 20:30:07 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 71 | Line 153 | bool SimInfo::removeMolecule(Molecule* mol) {
153   }
154  
155   bool SimInfo::removeMolecule(Molecule* mol) {
156 <    std::vector<Molecule*>::iterator i;
157 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
156 >    MoleculeIterator i;
157 >    i = molecules_.find(mol->getGlobalIndex());
158  
159      if (i != molecules_.end() ) {
160 <        molecules_.push_back(mol);
160 >
161 >        assert(mol == i->second);
162 >        
163          nAtoms_ -= mol->getNAtoms();
164          nBonds_ -= mol->getNBonds();
165          nBends_ -= mol->getNBends();
# Line 85 | Line 169 | bool SimInfo::removeMolecule(Molecule* mol) {
169          nCutoffGroups_ -= mol->getNCutoffGroups();
170          nConstraints_ -= mol->getNConstraints();
171  
172 +        molecules_.erase(mol->getGlobalIndex());
173 +
174 +        delete mol;
175 +        
176          return true;
177      } else {
178          return false;
# Line 94 | Line 182 | Molecule* SimInfo::beginMolecule(std::vector<Molecule*
182   }    
183  
184          
185 < Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
185 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
186      i = molecules_.begin();
187 <    return i == molecules_.end() ? NULL : *i;
187 >    return i == molecules_.end() ? NULL : i->second;
188   }    
189  
190 < Molecule* SimInfo::nextMolecule(std::vector<Molecule*>::iterator& i) {
190 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
191      ++i;
192 <    return i == molecules_.end() ? NULL : *i;    
192 >    return i == molecules_.end() ? NULL : i->second;    
193   }
194  
195  
196   void SimInfo::calcNdf() {
197      int ndf_local;
198 <    std::vector<Molecule*>::iterator i;
198 >    MoleculeIterator i;
199      std::vector<StuntDouble*>::iterator j;
200      Molecule* mol;
201      StuntDouble* integrableObject;
# Line 140 | Line 228 | void SimInfo::calcNdf() {
228      ndf_ = ndf_local;
229   #endif
230  
231 <    // nZconstraints is global, as are the 3 COM translations for the
231 >    // nZconstraints_ is global, as are the 3 COM translations for the
232      // entire system:
233 <    ndf_ = ndf_ - 3 - nZconstraints;
233 >    ndf_ = ndf_ - 3 - nZconstraint_;
234  
235   }
236  
237   void SimInfo::calcNdfRaw() {
238      int ndfRaw_local;
239  
240 <    std::vector<Molecule*>::iterator i;
240 >    MoleculeIterator i;
241      std::vector<StuntDouble*>::iterator j;
242      Molecule* mol;
243      StuntDouble* integrableObject;
# Line 193 | Line 281 | void SimInfo::calcNdfTrans() {
281      ndfTrans_ = ndfTrans_local;
282   #endif
283  
284 <    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
284 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
285  
286   }
287  
# Line 225 | Line 313 | void SimInfo::addExcludePairs(Molecule* mol) {
313          exclude_.addPair(b, c);        
314      }
315  
316 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
316 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
317          a = torsion->getAtomA()->getGlobalIndex();
318          b = torsion->getAtomB()->getGlobalIndex();        
319          c = torsion->getAtomC()->getGlobalIndex();        
# Line 270 | Line 358 | void SimInfo::removeExcludePairs(Molecule* mol) {
358          exclude_.removePair(b, c);        
359      }
360  
361 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
361 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
362          a = torsion->getAtomA()->getGlobalIndex();
363          b = torsion->getAtomB()->getGlobalIndex();        
364          c = torsion->getAtomC()->getGlobalIndex();        
# Line 282 | Line 370 | void SimInfo::removeExcludePairs(Molecule* mol) {
370          exclude_.removePair(b, c);
371          exclude_.removePair(b, d);
372          exclude_.removePair(c, d);        
373 +    }
374 +
375 + }
376 +
377 +
378 + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
379 +    int curStampId;
380 +
381 +    //index from 0
382 +    curStampId = molStampIds_.size();
383 +
384 +    moleculeStamps_.push_back(molStamp);
385 +    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
386 + }
387 +
388 + void SimInfo::update() {
389 +
390 +    setupSimType();
391 +
392 + #ifdef IS_MPI
393 +    setupFortranParallel();
394 + #endif
395 +
396 +    setupFortranSim();
397 +
398 +    //setup fortran force field
399 +    /** @deprecate */    
400 +    int isError = 0;
401 +    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
402 +    if(isError){
403 +        sprintf( painCave.errMsg,
404 +         "ForceField error: There was an error initializing the forceField in fortran.\n" );
405 +        painCave.isFatal = 1;
406 +        simError();
407 +    }
408 +  
409 +    
410 +    setupCutoff();
411 +
412 +    calcNdf();
413 +    calcNdfRaw();
414 +    calcNdfTrans();
415 +
416 +    fortranInitialized_ = true;
417 + }
418 +
419 + std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
420 +    SimInfo::MoleculeIterator mi;
421 +    Molecule* mol;
422 +    Molecule::AtomIterator ai;
423 +    Atom* atom;
424 +    std::set<AtomType*> atomTypes;
425 +
426 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
427 +
428 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
429 +            atomTypes.insert(atom->getAtomType());
430 +        }
431 +        
432 +    }
433 +
434 +    return atomTypes;        
435 + }
436 +
437 + void SimInfo::setupSimType() {
438 +    std::set<AtomType*>::iterator i;
439 +    std::set<AtomType*> atomTypes;
440 +    atomTypes = getUniqueAtomTypes();
441 +    
442 +    int useLennardJones = 0;
443 +    int useElectrostatic = 0;
444 +    int useEAM = 0;
445 +    int useCharge = 0;
446 +    int useDirectional = 0;
447 +    int useDipole = 0;
448 +    int useGayBerne = 0;
449 +    int useSticky = 0;
450 +    int useShape = 0;
451 +    int useFLARB = 0; //it is not in AtomType yet
452 +    int useDirectionalAtom = 0;    
453 +    int useElectrostatics = 0;
454 +    //usePBC and useRF are from simParams
455 +    bool usePBC = simParams_->getPBC();
456 +    bool useRF = simParams_->getUseRF();
457 +
458 +    //loop over all of the atom types
459 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
460 +        useLennardJones |= (*i)->isLennardJones();
461 +        useElectrostatic |= (*i)->isElectrostatic();
462 +        useEAM |= (*i)->isEAM();
463 +        useCharge |= (*i)->isCharge();
464 +        useDirectional |= (*i)->isDirectional();
465 +        useDipole |= (*i)->isDipole();
466 +        useGayBerne |= (*i)->isGayBerne();
467 +        useSticky |= (*i)->isSticky();
468 +        useShape |= (*i)->isShape();
469 +    }
470 +
471 +    if (useSticky || useDipole || useGayBerne || useShape) {
472 +        useDirectionalAtom = 1;
473 +    }
474 +
475 +    if (useCharge || useDipole) {
476 +        useElectrostatics = 1;
477 +    }
478 +
479 + #ifdef IS_MPI    
480 +    int temp;
481 +
482 +    temp = usePBC;
483 +    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
484 +
485 +    temp = useDirectionalAtom;
486 +    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
487 +
488 +    temp = useLennardJones;
489 +    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
490 +
491 +    temp = useElectrostatics;
492 +    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
493 +
494 +    temp = useCharge;
495 +    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
496 +
497 +    temp = useDipole;
498 +    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
499 +
500 +    temp = useSticky;
501 +    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
502 +
503 +    temp = useGayBerne;
504 +    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
505 +
506 +    temp = useEAM;
507 +    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
508 +
509 +    temp = useShape;
510 +    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
511 +
512 +    temp = useFLARB;
513 +    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
514 +
515 +    temp = useRF;
516 +    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
517 +    
518 + #endif
519 +
520 +    fInfo_.SIM_uses_PBC = usePBC;    
521 +    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
522 +    fInfo_.SIM_uses_LennardJones = useLennardJones;
523 +    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
524 +    fInfo_.SIM_uses_Charges = useCharge;
525 +    fInfo_.SIM_uses_Dipoles = useDipole;
526 +    fInfo_.SIM_uses_Sticky = useSticky;
527 +    fInfo_.SIM_uses_GayBerne = useGayBerne;
528 +    fInfo_.SIM_uses_EAM = useEAM;
529 +    fInfo_.SIM_uses_Shapes = useShape;
530 +    fInfo_.SIM_uses_FLARB = useFLARB;
531 +    fInfo_.SIM_uses_RF = useRF;
532 +
533 +    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
534 +
535 +        if (simParams_->haveDielectric()) {
536 +            fInfo_.dielect = simParams_->getDielectric();
537 +        } else {
538 +            sprintf(painCave.errMsg,
539 +                    "SimSetup Error: No Dielectric constant was set.\n"
540 +                    "\tYou are trying to use Reaction Field without"
541 +                    "\tsetting a dielectric constant!\n");
542 +            painCave.isFatal = 1;
543 +            simError();
544 +        }
545 +        
546 +    } else {
547 +        fInfo_.dielect = 0.0;
548      }
549  
550   }
551  
552 + void SimInfo::setupFortranSim() {
553 +    int isError;
554 +    int nExclude;
555 +    std::vector<int> fortranGlobalGroupMembership;
556 +    
557 +    nExclude = exclude_.getSize();
558 +    isError = 0;
559  
560 +    //globalGroupMembership_ is filled by SimCreator    
561 +    for (int i = 0; i < nGlobalAtoms_; i++) {
562 +        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
563 +    }
564 +
565 +    //calculate mass ratio of cutoff group
566 +    std::vector<double> mfact;
567 +    SimInfo::MoleculeIterator mi;
568 +    Molecule* mol;
569 +    Molecule::CutoffGroupIterator ci;
570 +    CutoffGroup* cg;
571 +    Molecule::AtomIterator ai;
572 +    Atom* atom;
573 +    double totalMass;
574 +
575 +    //to avoid memory reallocation, reserve enough space for mfact
576 +    mfact.reserve(getNCutoffGroups());
577 +    
578 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
579 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
580 +
581 +            totalMass = cg->getMass();
582 +            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
583 +                        mfact.push_back(atom->getMass()/totalMass);
584 +            }
585 +
586 +        }      
587 +    }
588 +
589 +    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
590 +    std::vector<int> identArray;
591 +
592 +    //to avoid memory reallocation, reserve enough space identArray
593 +    identArray.reserve(getNAtoms());
594 +    
595 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
596 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
597 +            identArray.push_back(atom->getIdent());
598 +        }
599 +    }    
600 +
601 +    //fill molMembershipArray
602 +    //molMembershipArray is filled by SimCreator    
603 +    std::vector<int> molMembershipArray(nGlobalAtoms_);
604 +    for (int i = 0; i < nGlobalAtoms_; i++) {
605 +        molMembershipArray.push_back(globalMolMembership_[i] + 1);
606 +    }
607 +    
608 +    //setup fortran simulation
609 +    //gloalExcludes and molMembershipArray should go away (They are never used)
610 +    //why the hell fortran need to know molecule?
611 +    //OOPSE = Object-Obfuscated Parallel Simulation Engine
612 +    int nGlobalExcludes = 0;
613 +    int* globalExcludes = NULL;
614 +    int* excludeList = exclude_.getExcludeList();
615 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
616 +                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
617 +                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
618 +
619 +    if( isError ){
620 +
621 +        sprintf( painCave.errMsg,
622 +                 "There was an error setting the simulation information in fortran.\n" );
623 +        painCave.isFatal = 1;
624 +        painCave.severity = OOPSE_ERROR;
625 +        simError();
626 +    }
627 +
628 + #ifdef IS_MPI
629 +    sprintf( checkPointMsg,
630 +       "succesfully sent the simulation information to fortran.\n");
631 +    MPIcheckPoint();
632 + #endif // is_mpi
633 + }
634 +
635 +
636 + #ifdef IS_MPI
637 + void SimInfo::setupFortranParallel() {
638 +    
639 +    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
640 +    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
641 +    std::vector<int> localToGlobalCutoffGroupIndex;
642 +    SimInfo::MoleculeIterator mi;
643 +    Molecule::AtomIterator ai;
644 +    Molecule::CutoffGroupIterator ci;
645 +    Molecule* mol;
646 +    Atom* atom;
647 +    CutoffGroup* cg;
648 +    mpiSimData parallelData;
649 +    int isError;
650 +
651 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
652 +
653 +        //local index(index in DataStorge) of atom is important
654 +        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
655 +            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
656 +        }
657 +
658 +        //local index of cutoff group is trivial, it only depends on the order of travesing
659 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
660 +            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
661 +        }        
662 +        
663 +    }
664 +
665 +    //fill up mpiSimData struct
666 +    parallelData.nMolGlobal = getNGlobalMolecules();
667 +    parallelData.nMolLocal = getNMolecules();
668 +    parallelData.nAtomsGlobal = getNGlobalAtoms();
669 +    parallelData.nAtomsLocal = getNAtoms();
670 +    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
671 +    parallelData.nGroupsLocal = getNCutoffGroups();
672 +    parallelData.myNode = worldRank;
673 +    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
674 +
675 +    //pass mpiSimData struct and index arrays to fortran
676 +    setFsimParallel(parallelData, &(parallelData->nAtomsLocal),
677 +                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
678 +                    &localToGlobalCutoffGroupIndex[0], &isError);
679 +
680 +    if (isError) {
681 +        sprintf(painCave.errMsg,
682 +                "mpiRefresh errror: fortran didn't like something we gave it.\n");
683 +        painCave.isFatal = 1;
684 +        simError();
685 +    }
686 +
687 +    sprintf(checkPointMsg, " mpiRefresh successful.\n");
688 +    MPIcheckPoint();
689 +
690 +
691 + }
692 +
693 + #endif
694 +
695 + double SimInfo::calcMaxCutoffRadius() {
696 +
697 +
698 +    std::set<AtomType*> atomTypes;
699 +    std::set<AtomType*>::iterator i;
700 +    std::vector<double> cutoffRadius;
701 +
702 +    //get the unique atom types
703 +    atomTypes = getUniqueAtomTypes();
704 +
705 +    //query the max cutoff radius among these atom types
706 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
707 +        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
708 +    }
709 +
710 +    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
711 + #ifdef IS_MPI
712 +    //pick the max cutoff radius among the processors
713 + #endif
714 +
715 +    return maxCutoffRadius;
716 + }
717 +
718 + void SimInfo::setupCutoff() {
719 +    double rcut_;  //cutoff radius
720 +    double rsw_; //switching radius
721 +    
722 +    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
723 +        
724 +        if (!simParams_->haveRcut()){
725 +            sprintf(painCave.errMsg,
726 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
727 +                "\tOOPSE will use a default value of 15.0 angstroms"
728 +                "\tfor the cutoffRadius.\n");
729 +            painCave.isFatal = 0;
730 +            simError();
731 +            rcut_ = 15.0;
732 +        } else{
733 +            rcut_ = simParams_->getRcut();
734 +        }
735 +
736 +        if (!simParams_->haveRsw()){
737 +            sprintf(painCave.errMsg,
738 +                "SimCreator Warning: No value was set for switchingRadius.\n"
739 +                "\tOOPSE will use a default value of\n"
740 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
741 +            painCave.isFatal = 0;
742 +            simError();
743 +            rsw_ = 0.95 * rcut_;
744 +        } else{
745 +            rsw_ = simParams_->getRsw();
746 +        }
747 +
748 +    } else {
749 +        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
750 +        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
751 +        
752 +        if (simParams_->haveRcut()) {
753 +            rcut_ = simParams_->getRcut();
754 +        } else {
755 +            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
756 +            rcut_ = calcMaxCutoffRadius();
757 +        }
758 +
759 +        if (simParams_->haveRsw()) {
760 +            rsw_  = simParams_->getRsw();
761 +        } else {
762 +            rsw_ = rcut_;
763 +        }
764 +    
765 +    }
766 +        
767 +    double rnblist = rcut_ + 1; // skin of neighbor list
768 +
769 +    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
770 +    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
771 + }
772 +
773 + void SimInfo::addProperty(GenericData* genData) {
774 +    properties_.addProperty(genData);  
775 + }
776 +
777 + void SimInfo::removeProperty(const std::string& propName) {
778 +    properties_.removeProperty(propName);  
779 + }
780 +
781 + void SimInfo::clearProperties() {
782 +    properties_.clearProperties();
783 + }
784 +
785 + std::vector<std::string> SimInfo::getPropertyNames() {
786 +    return properties_.getPropertyNames();  
787 + }
788 +      
789 + std::vector<GenericData*> SimInfo::getProperties() {
790 +    return properties_.getProperties();
791 + }
792 +
793 + GenericData* SimInfo::getPropertyByName(const std::string& propName) {
794 +    return properties_.getPropertyByName(propName);
795 + }
796 +
797 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
798 +    sman_ = sman;
799 +
800 +    Molecule* mol;
801 +    RigidBody* rb;
802 +    Atom* atom;
803 +    SimInfo::MoleculeIterator mi;
804 +    Molecule::RigidBodyIterator rbIter;
805 +    Molecule::AtomIterator atomIter;;
806 +
807 +    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
808 +        
809 +        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
810 +            atom->setSnapshotManager(sman_);
811 +        }
812 +        
813 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
814 +            rb->setSnapshotManager(sman_);
815 +        }
816 +    }    
817 +    
818 + }
819 +
820 + std::ostream& operator <<(std::ostream& o, SimInfo& info) {
821 +
822 +    return o;
823 + }
824 +
825   }//end namespace oopse
826 +

Diff Legend

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