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 1733 by tim, Fri Nov 12 06:19:04 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(const std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
45 <                                ForceField* ff, Globals* globals) :
46 <                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
47 <                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
48 <                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(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;
47 –    int nCutoffAtoms; // number of atoms belong to cutoff groups
48 –    int ngroups;          //total cutoff groups defined in meta-data file
56      MoleculeStamp* molStamp;
57      int nMolWithSameStamp;
58 <    CutoffGroupStamp* cgStamp;
59 <    int nAtomsIngroups;
60 <    int nCutoffGroupsInStamp;    
61 <
62 <    nGlobalAtoms_ =  0;
56 <    ngroups = 0;
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 <        
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          
66 –        nAtomsIngroups = 0;
67 –        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
68 –        
78          for (int j=0; j < nCutoffGroupsInStamp; j++) {
79              cgStamp = molStamp->getCutoffGroup(j);
80 <            nAtomsIngroups += cgStamp->getNMembers();
80 >            nAtomsInGroups += cgStamp->getNMembers();
81          }
82  
83 <        ngroups += *nMolWithSameStamp;
84 <        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
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;
104 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
105  
106 <    //initialize globalGroupMembership_, every element of this array will be 0
107 <    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
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  
# Line 98 | Line 123 | SimInfo::~SimInfo() {
123      MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
124      
125      delete sman_;
126 <    delete globals_;
126 >    delete simParams_;
127      delete forceField_;
128  
129   }
# Line 108 | Line 133 | bool SimInfo::addMolecule(Molecule* mol) {
133      MoleculeIterator i;
134  
135      i = molecules_.find(mol->getGlobalIndex());
136 <    if (i != molecules_.end() ) {
136 >    if (i == molecules_.end() ) {
137  
138 <        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
138 >        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
139          
140          nAtoms_ += mol->getNAtoms();
141          nBonds_ += mol->getNBonds();
# Line 121 | 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 144 | 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;
# Line 159 | Line 187 | Molecule* SimInfo::beginMolecule(MoleculeIterator& i)
187          
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(MoleculeIterator& i) {
194      ++i;
195 <    return i == molecules_.end() ? NULL : *i;    
195 >    return i == molecules_.end() ? NULL : i->second;    
196   }
197  
198  
# Line 205 | Line 233 | void SimInfo::calcNdf() {
233  
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  
# Line 256 | 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 288 | 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 333 | 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 357 | Line 385 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
385      curStampId = molStampIds_.size();
386  
387      moleculeStamps_.push_back(molStamp);
388 <    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
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 <    typename SimInfo::MoleculeIterator mi;
423 >    SimInfo::MoleculeIterator mi;
424      Molecule* mol;
425 <    typename Molecule::AtomIterator ai;
425 >    Molecule::AtomIterator ai;
426      Atom* atom;
427      std::set<AtomType*> atomTypes;
428  
# Line 408 | Line 454 | void SimInfo::setupSimType() {
454      int useFLARB = 0; //it is not in AtomType yet
455      int useDirectionalAtom = 0;    
456      int useElectrostatics = 0;
457 <    //usePBC and useRF are from globals
458 <    bool usePBC = globals_->getPBC();
459 <    bool useRF = globals_->getUseRF();
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();
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) {
# Line 488 | Line 534 | void SimInfo::setupSimType() {
534      fInfo_.SIM_uses_RF = useRF;
535  
536      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
537 <        fInfo_.dielect = dielectric;
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      }
# Line 510 | Line 567 | void SimInfo::setupFortranSim() {
567  
568      //calculate mass ratio of cutoff group
569      std::vector<double> mfact;
570 <    typename SimInfo::MoleculeIterator mi;
570 >    SimInfo::MoleculeIterator mi;
571      Molecule* mol;
572 <    typename Molecule::CutoffGroupIterator ci;
572 >    Molecule::CutoffGroupIterator ci;
573      CutoffGroup* cg;
574 <    typename Molecule::AtomIterator ai;
574 >    Molecule::AtomIterator ai;
575      Atom* atom;
576      double totalMass;
577  
# Line 544 | Line 601 | void SimInfo::setupFortranSim() {
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 <    
616 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
617 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
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 ){
# Line 576 | Line 642 | void SimInfo::setupFortranParallel() {
642      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
643      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
644      std::vector<int> localToGlobalCutoffGroupIndex;
645 <    typename SimInfo::MoleculeIterator mi;
646 <    typename Molecule::AtomIterator ai;
647 <    typename Molecule::CutoffGroupIterator ci;
645 >    SimInfo::MoleculeIterator mi;
646 >    Molecule::AtomIterator ai;
647 >    Molecule::CutoffGroupIterator ci;
648      Molecule* mol;
649      Atom* atom;
650      CutoffGroup* cg;
# Line 629 | Line 695 | void SimInfo::addProperty(GenericData* genData) {
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   }
# Line 653 | Line 797 | GenericData* SimInfo::getPropertyByName(const std::str
797      return properties_.getPropertyByName(propName);
798   }
799  
800 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
801 +    sman_ = sman;
802  
803 < std::ostream& operator <<(ostream& o, SimInfo& info) {
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