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 1738 by tim, Sat Nov 13 05:08:12 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),
49 <                                fortranInitialized_(false) {
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;
48 –    int nCutoffAtoms; // number of atoms belong to cutoff groups
49 –    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;
57 <    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          
67 –        nAtomsIngroups = 0;
68 –        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
69 –        
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 99 | Line 123 | SimInfo::~SimInfo() {
123      MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
124      
125      delete sman_;
126 <    delete globals_;
126 >    delete simParams_;
127      delete forceField_;
128  
129   }
# Line 109 | 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 122 | 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 145 | 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 206 | 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 257 | 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 289 | 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 334 | 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 358 | 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() {
# Line 371 | Line 398 | void SimInfo::update() {
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  
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 –    
415      calcNdf();
416      calcNdfRaw();
417      calcNdfTrans();
# Line 392 | Line 420 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
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 426 | 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 506 | 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 528 | 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 566 | Line 605 | void SimInfo::setupFortranSim() {
605      //molMembershipArray is filled by SimCreator    
606      std::vector<int> molMembershipArray(nGlobalAtoms_);
607      for (int i = 0; i < nGlobalAtoms_; i++) {
608 <        molMembershipArray.push_back(globalMolMembership_[i] + 1);
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 601 | 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 657 | Line 698 | double SimInfo::calcMaxCutoffRadius() {
698   double SimInfo::calcMaxCutoffRadius() {
699  
700  
701 <    std::vector<AtomType*> atomTypes;
702 <    std::vector<AtomType*>::iterator i;
701 >    std::set<AtomType*> atomTypes;
702 >    std::set<AtomType*>::iterator i;
703      std::vector<double> cutoffRadius;
704  
705      //get the unique atom types
# Line 669 | Line 710 | double SimInfo::calcMaxCutoffRadius() {
710          cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
711      }
712  
713 <    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
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
# Line 683 | Line 724 | void SimInfo::setupCutoff() {
724      
725      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
726          
727 <        if (!globals_->haveRcut()){
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"
# Line 692 | Line 733 | void SimInfo::setupCutoff() {
733              simError();
734              rcut_ = 15.0;
735          } else{
736 <            rcut_ = globals_->getRcut();
736 >            rcut_ = simParams_->getRcut();
737          }
738  
739 <        if (!globals_->haveRsw()){
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"
# Line 704 | Line 745 | void SimInfo::setupCutoff() {
745              simError();
746              rsw_ = 0.95 * rcut_;
747          } else{
748 <            rsw_ = globals_->getRsw();
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 (globals_->haveRcut()) {
756 <            rcut_ = globals_->getRcut();
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 (globals_->haveRsw()) {
763 <            rsw_  = globals_->getRsw()
762 >        if (simParams_->haveRsw()) {
763 >            rsw_  = simParams_->getRsw();
764          } else {
765              rsw_ = rcut_;
766          }
# Line 756 | 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