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 1738 by tim, Sat Nov 13 05:08:12 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(const std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
50 <                                ForceField* ff, Globals* globals) :
51 <                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
52 <                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
53 <                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL),
54 <                                fortranInitialized_(false) {
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;
48 –    int nCutoffAtoms; // number of atoms belong to cutoff groups
49 –    int ngroups;          //total cutoff groups defined in meta-data file
61      MoleculeStamp* molStamp;
62      int nMolWithSameStamp;
63 <    CutoffGroupStamp* cgStamp;
64 <    int nAtomsIngroups;
65 <    int nCutoffGroupsInStamp;    
66 <
67 <    nGlobalAtoms_ =  0;
57 <    ngroups = 0;
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 <        
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          
67 –        nAtomsIngroups = 0;
68 –        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
69 –        
83          for (int j=0; j < nCutoffGroupsInStamp; j++) {
84              cgStamp = molStamp->getCutoffGroup(j);
85 <            nAtomsIngroups += cgStamp->getNMembers();
85 >            nAtomsInGroups += cgStamp->getNMembers();
86          }
87  
88 <        ngroups += *nMolWithSameStamp;
89 <        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
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;
109 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
110  
111 <    //initialize globalGroupMembership_, every element of this array will be 0
112 <    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
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  
# Line 99 | Line 128 | SimInfo::~SimInfo() {
128      MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
129      
130      delete sman_;
131 <    delete globals_;
131 >    delete simParams_;
132      delete forceField_;
133  
134   }
# Line 109 | Line 138 | bool SimInfo::addMolecule(Molecule* mol) {
138      MoleculeIterator i;
139  
140      i = molecules_.find(mol->getGlobalIndex());
141 <    if (i != molecules_.end() ) {
141 >    if (i == molecules_.end() ) {
142  
143 <        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
143 >        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
144          
145          nAtoms_ += mol->getNAtoms();
146          nBonds_ += mol->getNBonds();
# Line 122 | 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 145 | 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;
# Line 206 | Line 238 | void SimInfo::calcNdf() {
238  
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  
# Line 257 | 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 289 | 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 334 | 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 358 | Line 390 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
390      curStampId = molStampIds_.size();
391  
392      moleculeStamps_.push_back(molStamp);
393 <    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
393 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
394   }
395  
396   void SimInfo::update() {
# Line 371 | Line 403 | void SimInfo::update() {
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  
376 –    //notify fortran whether reaction field is used or not. It is deprecated now
377 –    //int isError = 0;
378 –    //initFortranFF( &useReactionField, &isError );
379 –
380 –    //if(isError){
381 –    //    sprintf( painCave.errMsg,
382 –    //    "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" );
383 –    //    painCave.isFatal = 1;
384 –    //    simError();
385 –    //}
386 –    
420      calcNdf();
421      calcNdfRaw();
422      calcNdfTrans();
# Line 392 | Line 425 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
425   }
426  
427   std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
428 <    typename SimInfo::MoleculeIterator mi;
428 >    SimInfo::MoleculeIterator mi;
429      Molecule* mol;
430 <    typename Molecule::AtomIterator ai;
430 >    Molecule::AtomIterator ai;
431      Atom* atom;
432      std::set<AtomType*> atomTypes;
433  
# Line 426 | Line 459 | void SimInfo::setupSimType() {
459      int useFLARB = 0; //it is not in AtomType yet
460      int useDirectionalAtom = 0;    
461      int useElectrostatics = 0;
462 <    //usePBC and useRF are from globals
463 <    bool usePBC = globals_->getPBC();
464 <    bool useRF = globals_->getUseRF();
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();
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) {
# Line 506 | Line 539 | void SimInfo::setupSimType() {
539      fInfo_.SIM_uses_RF = useRF;
540  
541      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
542 <        fInfo_.dielect = dielectric;
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      }
# Line 528 | Line 572 | void SimInfo::setupFortranSim() {
572  
573      //calculate mass ratio of cutoff group
574      std::vector<double> mfact;
575 <    typename SimInfo::MoleculeIterator mi;
575 >    SimInfo::MoleculeIterator mi;
576      Molecule* mol;
577 <    typename Molecule::CutoffGroupIterator ci;
577 >    Molecule::CutoffGroupIterator ci;
578      CutoffGroup* cg;
579 <    typename Molecule::AtomIterator ai;
579 >    Molecule::AtomIterator ai;
580      Atom* atom;
581      double totalMass;
582  
# Line 566 | Line 610 | void SimInfo::setupFortranSim() {
610      //molMembershipArray is filled by SimCreator    
611      std::vector<int> molMembershipArray(nGlobalAtoms_);
612      for (int i = 0; i < nGlobalAtoms_; i++) {
613 <        molMembershipArray.push_back(globalMolMembership_[i] + 1);
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 <    
621 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
622 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
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 ){
# Line 601 | Line 647 | void SimInfo::setupFortranParallel() {
647      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
648      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
649      std::vector<int> localToGlobalCutoffGroupIndex;
650 <    typename SimInfo::MoleculeIterator mi;
651 <    typename Molecule::AtomIterator ai;
652 <    typename Molecule::CutoffGroupIterator ci;
650 >    SimInfo::MoleculeIterator mi;
651 >    Molecule::AtomIterator ai;
652 >    Molecule::CutoffGroupIterator ci;
653      Molecule* mol;
654      Atom* atom;
655      CutoffGroup* cg;
# Line 632 | Line 678 | void SimInfo::setupFortranParallel() {
678      parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
679      parallelData.nGroupsLocal = getNCutoffGroups();
680      parallelData.myNode = worldRank;
681 <    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
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),
684 >    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
685 >                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
686                      &localToGlobalCutoffGroupIndex[0], &isError);
687  
688      if (isError) {
# Line 657 | Line 703 | double SimInfo::calcMaxCutoffRadius() {
703   double SimInfo::calcMaxCutoffRadius() {
704  
705  
706 <    std::vector<AtomType*> atomTypes;
707 <    std::vector<AtomType*>::iterator i;
706 >    std::set<AtomType*> atomTypes;
707 >    std::set<AtomType*>::iterator i;
708      std::vector<double> cutoffRadius;
709  
710      //get the unique atom types
# Line 669 | Line 715 | double SimInfo::calcMaxCutoffRadius() {
715          cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
716      }
717  
718 <    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
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
# Line 683 | Line 729 | void SimInfo::setupCutoff() {
729      
730      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
731          
732 <        if (!globals_->haveRcut()){
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"
# Line 692 | Line 738 | void SimInfo::setupCutoff() {
738              simError();
739              rcut_ = 15.0;
740          } else{
741 <            rcut_ = globals_->getRcut();
741 >            rcut_ = simParams_->getRcut();
742          }
743  
744 <        if (!globals_->haveRsw()){
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"
# Line 704 | Line 750 | void SimInfo::setupCutoff() {
750              simError();
751              rsw_ = 0.95 * rcut_;
752          } else{
753 <            rsw_ = globals_->getRsw();
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 (globals_->haveRcut()) {
761 <            rcut_ = globals_->getRcut();
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 (globals_->haveRsw()) {
768 <            rsw_  = globals_->getRsw()
767 >        if (simParams_->haveRsw()) {
768 >            rsw_  = simParams_->getRsw();
769          } else {
770              rsw_ = rcut_;
771          }
# Line 756 | Line 802 | GenericData* SimInfo::getPropertyByName(const std::str
802      return properties_.getPropertyByName(propName);
803   }
804  
805 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
806 +    sman_ = sman;
807  
808 < std::ostream& operator <<(ostream& o, SimInfo& info) {
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