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 1735 by tim, Fri Nov 12 17:40:03 2004 UTC vs.
Revision 1824 by tim, Thu Dec 2 03:12:25 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,
44 > SimInfo::SimInfo(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) {
48 >                                nCutoffGroups_(0), nConstraints_(0), nZconstraint_(0), sman_(NULL),
49 >                                fortranInitialized_(false) {
50  
51      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
52      MoleculeStamp* molStamp;
53      int nMolWithSameStamp;
54 +    int nCutoffAtoms; // number of atoms belong to cutoff groups
55 +    int nGroups;          //total cutoff groups defined in meta-data file
56      CutoffGroupStamp* cgStamp;
57 <    int nAtomsIngroups;
58 <    int nCutoffGroupsInStamp;    
59 <
57 >    int nAtomsInGroups;
58 >    int nCutoffGroupsInStamp;
59 >    
60 >    RigidBodyStamp* rbStamp;
61 >    int nAtomsInRigidBodies;
62 >    int nRigidBodiesInStamp;
63 >    int nRigidAtoms;
64 >    int nRigidBodies;
65 >        
66      nGlobalAtoms_ =  0;
56 –    ngroups = 0;
67      
68 +    nGroups = 0;
69 +    nCutoffAtoms = 0;
70 +    nRigidBodies = 0;
71 +    
72      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
73          molStamp = i->first;
74          nMolWithSameStamp = i->second;
75          
76          addMoleculeStamp(molStamp, nMolWithSameStamp);
77 <        
77 >
78 >        //calculate atoms in molecules
79          nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
80 <        
81 <        nAtomsIngroups = 0;
80 >
81 >
82 >        //calculate atoms in cutoff groups
83 >        nAtomsInGroups = 0;
84          nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
85          
86          for (int j=0; j < nCutoffGroupsInStamp; j++) {
87              cgStamp = molStamp->getCutoffGroup(j);
88 <            nAtomsIngroups += cgStamp->getNMembers();
88 >            nAtomsInGroups += cgStamp->getNMembers();
89          }
90  
91 <        ngroups += *nMolWithSameStamp;
92 <        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
91 >        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
92 >        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
93 >
94 >        //calculate atoms in rigid bodies
95 >        nAtomsInRigidBodies = 0;
96 >        nRigidBodiesInStamp = molStamp->getNCutoffGroups();
97 >        
98 >        for (int j=0; j < nRigidBodiesInStamp; j++) {
99 >            rbStamp = molStamp->getRigidBody(j);
100 >            nRigidBodiesInStamp += rbStamp->getNMembers();
101 >        }
102 >
103 >        nRigidBodies += nRigidBodiesInStamp * nMolWithSameStamp;
104 >        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
105 >        
106      }
107  
108      //every free atom (atom does not belong to cutoff groups) is a cutoff group
109      //therefore the total number of cutoff groups in the system is equal to
110      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
111      //file plus the number of cutoff groups defined in meta-data file
112 <    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups;
112 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
113  
114 +    //every free atom (atom does not belong to rigid bodies) is a rigid body
115 +    //therefore the total number of cutoff groups in the system is equal to
116 +    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
117 +    //file plus the number of  rigid bodies defined in meta-data file
118 +    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nRigidBodies;
119 +
120      //initialize globalGroupMembership_, every element of this array will be 0
121      globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
122  
# Line 159 | Line 195 | Molecule* SimInfo::beginMolecule(MoleculeIterator& i)
195          
196   Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
197      i = molecules_.begin();
198 <    return i == molecules_.end() ? NULL : *i;
198 >    return i == molecules_.end() ? NULL : i->second;
199   }    
200  
201   Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
202      ++i;
203 <    return i == molecules_.end() ? NULL : *i;    
203 >    return i == molecules_.end() ? NULL : i->second;    
204   }
205  
206  
# Line 205 | Line 241 | void SimInfo::calcNdf() {
241  
242      // nZconstraints_ is global, as are the 3 COM translations for the
243      // entire system:
244 <    ndf_ = ndf_ - 3 - nZconstraints_;
244 >    ndf_ = ndf_ - 3 - nZconstraint_;
245  
246   }
247  
# Line 256 | Line 292 | void SimInfo::calcNdfTrans() {
292      ndfTrans_ = ndfTrans_local;
293   #endif
294  
295 <    ndfTrans_ = ndfTrans_ - 3 - nZconstraints_;
295 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
296  
297   }
298  
# Line 288 | Line 324 | void SimInfo::addExcludePairs(Molecule* mol) {
324          exclude_.addPair(b, c);        
325      }
326  
327 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
327 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
328          a = torsion->getAtomA()->getGlobalIndex();
329          b = torsion->getAtomB()->getGlobalIndex();        
330          c = torsion->getAtomC()->getGlobalIndex();        
# Line 333 | Line 369 | void SimInfo::removeExcludePairs(Molecule* mol) {
369          exclude_.removePair(b, c);        
370      }
371  
372 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
372 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
373          a = torsion->getAtomA()->getGlobalIndex();
374          b = torsion->getAtomB()->getGlobalIndex();        
375          c = torsion->getAtomC()->getGlobalIndex();        
# Line 357 | Line 393 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
393      curStampId = molStampIds_.size();
394  
395      moleculeStamps_.push_back(molStamp);
396 <    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
396 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
397   }
398  
399   void SimInfo::update() {
400  
401 +    setupSimType();
402 +
403   #ifdef IS_MPI
404      setupFortranParallel();
405   #endif
406  
407      setupFortranSim();
408  
409 +    //setup fortran force field
410 +    /** @deprecate */    
411 +    int isError = 0;
412 +    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
413 +    if(isError){
414 +        sprintf( painCave.errMsg,
415 +         "ForceField error: There was an error initializing the forceField in fortran.\n" );
416 +        painCave.isFatal = 1;
417 +        simError();
418 +    }
419 +  
420 +    
421 +    setupCutoff();
422 +
423      calcNdf();
424      calcNdfRaw();
425      calcNdfTrans();
426 +
427 +    fortranInitialized_ = true;
428   }
429  
430   std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
431 <    typename SimInfo::MoleculeIterator mi;
431 >    SimInfo::MoleculeIterator mi;
432      Molecule* mol;
433 <    typename Molecule::AtomIterator ai;
433 >    Molecule::AtomIterator ai;
434      Atom* atom;
435      std::set<AtomType*> atomTypes;
436  
# Line 414 | Line 468 | void SimInfo::setupSimType() {
468  
469      //loop over all of the atom types
470      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
471 <        useLennardJones |= i->isLennardJones();
472 <        useElectrostatic |= i->isElectrostatic();
473 <        useEAM |= i->isEAM();
474 <        useCharge |= i->isCharge();
475 <        useDirectional |= i->isDirectional();
476 <        useDipole |= i->isDipole();
477 <        useGayBerne |= i->isGayBerne();
478 <        useSticky |= i->isSticky();
479 <        useShape |= i->isShape();
471 >        useLennardJones |= (*i)->isLennardJones();
472 >        useElectrostatic |= (*i)->isElectrostatic();
473 >        useEAM |= (*i)->isEAM();
474 >        useCharge |= (*i)->isCharge();
475 >        useDirectional |= (*i)->isDirectional();
476 >        useDipole |= (*i)->isDipole();
477 >        useGayBerne |= (*i)->isGayBerne();
478 >        useSticky |= (*i)->isSticky();
479 >        useShape |= (*i)->isShape();
480      }
481  
482      if (useSticky || useDipole || useGayBerne || useShape) {
# Line 488 | Line 542 | void SimInfo::setupSimType() {
542      fInfo_.SIM_uses_RF = useRF;
543  
544      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
545 <        fInfo_.dielect = dielectric;
545 >
546 >        if (globals_->haveDielectric()) {
547 >            fInfo_.dielect = globals_->getDielectric();
548 >        } else {
549 >            sprintf(painCave.errMsg,
550 >                    "SimSetup Error: No Dielectric constant was set.\n"
551 >                    "\tYou are trying to use Reaction Field without"
552 >                    "\tsetting a dielectric constant!\n");
553 >            painCave.isFatal = 1;
554 >            simError();
555 >        }
556 >        
557      } else {
558          fInfo_.dielect = 0.0;
559      }
# Line 510 | Line 575 | void SimInfo::setupFortranSim() {
575  
576      //calculate mass ratio of cutoff group
577      std::vector<double> mfact;
578 <    typename SimInfo::MoleculeIterator mi;
578 >    SimInfo::MoleculeIterator mi;
579      Molecule* mol;
580 <    typename Molecule::CutoffGroupIterator ci;
580 >    Molecule::CutoffGroupIterator ci;
581      CutoffGroup* cg;
582 <    typename Molecule::AtomIterator ai;
582 >    Molecule::AtomIterator ai;
583      Atom* atom;
584      double totalMass;
585  
# Line 555 | Line 620 | void SimInfo::setupFortranSim() {
620      //gloalExcludes and molMembershipArray should go away (They are never used)
621      //why the hell fortran need to know molecule?
622      //OOPSE = Object-Obfuscated Parallel Simulation Engine
623 <    
624 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
625 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
623 >    int nGlobalExcludes = 0;
624 >    int* globalExcludes = NULL;
625 >    int* excludeList = exclude_.getExcludeList();
626 >    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
627 >                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
628                    &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
629  
630      if( isError ){
# Line 583 | Line 650 | void SimInfo::setupFortranParallel() {
650      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
651      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
652      std::vector<int> localToGlobalCutoffGroupIndex;
653 <    typename SimInfo::MoleculeIterator mi;
654 <    typename Molecule::AtomIterator ai;
655 <    typename Molecule::CutoffGroupIterator ci;
653 >    SimInfo::MoleculeIterator mi;
654 >    Molecule::AtomIterator ai;
655 >    Molecule::CutoffGroupIterator ci;
656      Molecule* mol;
657      Atom* atom;
658      CutoffGroup* cg;
# Line 639 | Line 706 | double SimInfo::calcMaxCutoffRadius() {
706   double SimInfo::calcMaxCutoffRadius() {
707  
708  
709 <    std::vector<AtomType*> atomTypes;
710 <    std::vector<AtomType*>::iterator i;
709 >    std::set<AtomType*> atomTypes;
710 >    std::set<AtomType*>::iterator i;
711      std::vector<double> cutoffRadius;
712  
713      //get the unique atom types
# Line 651 | Line 718 | double SimInfo::calcMaxCutoffRadius() {
718          cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
719      }
720  
721 <    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
721 >    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
722   #ifdef IS_MPI
723      //pick the max cutoff radius among the processors
724   #endif
725  
726      return maxCutoffRadius;
727 + }
728 +
729 + void SimInfo::setupCutoff() {
730 +    double rcut_;  //cutoff radius
731 +    double rsw_; //switching radius
732 +    
733 +    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
734 +        
735 +        if (!globals_->haveRcut()){
736 +            sprintf(painCave.errMsg,
737 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
738 +                "\tOOPSE will use a default value of 15.0 angstroms"
739 +                "\tfor the cutoffRadius.\n");
740 +            painCave.isFatal = 0;
741 +            simError();
742 +            rcut_ = 15.0;
743 +        } else{
744 +            rcut_ = globals_->getRcut();
745 +        }
746 +
747 +        if (!globals_->haveRsw()){
748 +            sprintf(painCave.errMsg,
749 +                "SimCreator Warning: No value was set for switchingRadius.\n"
750 +                "\tOOPSE will use a default value of\n"
751 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
752 +            painCave.isFatal = 0;
753 +            simError();
754 +            rsw_ = 0.95 * rcut_;
755 +        } else{
756 +            rsw_ = globals_->getRsw();
757 +        }
758 +
759 +    } else {
760 +        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
761 +        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
762 +        
763 +        if (globals_->haveRcut()) {
764 +            rcut_ = globals_->getRcut();
765 +        } else {
766 +            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
767 +            rcut_ = calcMaxCutoffRadius();
768 +        }
769 +
770 +        if (globals_->haveRsw()) {
771 +            rsw_  = globals_->getRsw();
772 +        } else {
773 +            rsw_ = rcut_;
774 +        }
775 +    
776 +    }
777 +        
778 +    double rnblist = rcut_ + 1; // skin of neighbor list
779 +
780 +    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
781 +    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
782   }
783  
784   void SimInfo::addProperty(GenericData* genData) {
# Line 683 | Line 805 | GenericData* SimInfo::getPropertyByName(const std::str
805      return properties_.getPropertyByName(propName);
806   }
807  
808 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
809 +    sman_ = sman;
810  
811 +    Molecule* mol;
812 +    RigidBody* rb;
813 +    Atom* atom;
814 +    SimInfo::MoleculeIterator mi;
815 +    Molecule::RigidBodyIterator rbIter;
816 +    Molecule::AtomIterator atomIter;;
817 +
818 +    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
819 +        
820 +        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
821 +            atom->setSnapshotManager(sman_);
822 +        }
823 +        
824 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
825 +            rb->setSnapshotManager(sman_);
826 +        }
827 +    }    
828 +    
829 + }
830 +
831   std::ostream& operator <<(ostream& o, SimInfo& info) {
832  
833      return o;

Diff Legend

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