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 1854 by tim, Sun Dec 5 22:02:42 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 206 | Line 230 | void SimInfo::calcNdf() {
230  
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  
# Line 257 | 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 289 | 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 334 | 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 358 | Line 382 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
382      curStampId = molStampIds_.size();
383  
384      moleculeStamps_.push_back(molStamp);
385 <    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
385 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
386   }
387  
388   void SimInfo::update() {
# Line 371 | Line 395 | void SimInfo::update() {
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  
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 –    
412      calcNdf();
413      calcNdfRaw();
414      calcNdfTrans();
# Line 392 | Line 417 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
417   }
418  
419   std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
420 <    typename SimInfo::MoleculeIterator mi;
420 >    SimInfo::MoleculeIterator mi;
421      Molecule* mol;
422 <    typename Molecule::AtomIterator ai;
422 >    Molecule::AtomIterator ai;
423      Atom* atom;
424      std::set<AtomType*> atomTypes;
425  
# Line 426 | Line 451 | void SimInfo::setupSimType() {
451      int useFLARB = 0; //it is not in AtomType yet
452      int useDirectionalAtom = 0;    
453      int useElectrostatics = 0;
454 <    //usePBC and useRF are from globals
455 <    bool usePBC = globals_->getPBC();
456 <    bool useRF = globals_->getUseRF();
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();
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) {
# Line 506 | Line 531 | void SimInfo::setupSimType() {
531      fInfo_.SIM_uses_RF = useRF;
532  
533      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
534 <        fInfo_.dielect = dielectric;
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      }
# Line 528 | Line 564 | void SimInfo::setupFortranSim() {
564  
565      //calculate mass ratio of cutoff group
566      std::vector<double> mfact;
567 <    typename SimInfo::MoleculeIterator mi;
567 >    SimInfo::MoleculeIterator mi;
568      Molecule* mol;
569 <    typename Molecule::CutoffGroupIterator ci;
569 >    Molecule::CutoffGroupIterator ci;
570      CutoffGroup* cg;
571 <    typename Molecule::AtomIterator ai;
571 >    Molecule::AtomIterator ai;
572      Atom* atom;
573      double totalMass;
574  
# Line 566 | Line 602 | void SimInfo::setupFortranSim() {
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);
605 >        molMembershipArray[i] = 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 <    
613 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
614 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
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 ){
# Line 601 | Line 639 | void SimInfo::setupFortranParallel() {
639      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
640      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
641      std::vector<int> localToGlobalCutoffGroupIndex;
642 <    typename SimInfo::MoleculeIterator mi;
643 <    typename Molecule::AtomIterator ai;
644 <    typename Molecule::CutoffGroupIterator ci;
642 >    SimInfo::MoleculeIterator mi;
643 >    Molecule::AtomIterator ai;
644 >    Molecule::CutoffGroupIterator ci;
645      Molecule* mol;
646      Atom* atom;
647      CutoffGroup* cg;
# Line 657 | Line 695 | double SimInfo::calcMaxCutoffRadius() {
695   double SimInfo::calcMaxCutoffRadius() {
696  
697  
698 <    std::vector<AtomType*> atomTypes;
699 <    std::vector<AtomType*>::iterator i;
698 >    std::set<AtomType*> atomTypes;
699 >    std::set<AtomType*>::iterator i;
700      std::vector<double> cutoffRadius;
701  
702      //get the unique atom types
# Line 669 | Line 707 | double SimInfo::calcMaxCutoffRadius() {
707          cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
708      }
709  
710 <    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
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
# Line 683 | Line 721 | void SimInfo::setupCutoff() {
721      
722      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
723          
724 <        if (!globals_->haveRcut()){
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"
# Line 692 | Line 730 | void SimInfo::setupCutoff() {
730              simError();
731              rcut_ = 15.0;
732          } else{
733 <            rcut_ = globals_->getRcut();
733 >            rcut_ = simParams_->getRcut();
734          }
735  
736 <        if (!globals_->haveRsw()){
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"
# Line 704 | Line 742 | void SimInfo::setupCutoff() {
742              simError();
743              rsw_ = 0.95 * rcut_;
744          } else{
745 <            rsw_ = globals_->getRsw();
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 (globals_->haveRcut()) {
753 <            rcut_ = globals_->getRcut();
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 (globals_->haveRsw()) {
760 <            rsw_  = globals_->getRsw()
759 >        if (simParams_->haveRsw()) {
760 >            rsw_  = simParams_->getRsw();
761          } else {
762              rsw_ = rcut_;
763          }
# Line 756 | Line 794 | GenericData* SimInfo::getPropertyByName(const std::str
794      return properties_.getPropertyByName(propName);
795   }
796  
797 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
798 +    sman_ = sman;
799  
800 < std::ostream& operator <<(ostream& o, SimInfo& info) {
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 + Vector3d SimInfo::getComVel(){
821 +    SimInfo::MoleculeIterator i;
822 +    Molecule* mol;
823 +
824 +    Vector3d comVel(0.0);
825 +    double totalMass = 0.0;
826 +    
827 +
828 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
829 +        double mass = mol->getMass();
830 +        totalMass += mass;
831 +        comVel += mass * mol->getComVel();
832 +    }  
833 +
834 +    comVel /= totalMass;
835 +
836 +    return comVel;
837 + }
838 +
839 + Vector3d SimInfo::getCom(){
840 +    SimInfo::MoleculeIterator i;
841 +    Molecule* mol;
842 +
843 +    Vector3d com(0.0);
844 +    double totalMass = 0.0;
845 +    
846 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
847 +        double mass = mol->getMass();
848 +        totalMass += mass;
849 +        com += mass * mol->getCom();
850 +    }  
851 +
852 +    com /= totalMass;
853 +
854 +    return com;
855 +
856 + }        
857 +
858 + std::ostream& operator <<(std::ostream& o, SimInfo& info) {
859 +
860      return o;
861   }
862  
863   }//end namespace oopse
864 +

Diff Legend

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