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 1807 by tim, Tue Nov 30 22:43:51 2004 UTC

# Line 33 | Line 33
33   #include <algorithm>
34  
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;
48 –    int nCutoffAtoms; // number of atoms belong to cutoff groups
49 –    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;
57 –    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 206 | 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 257 | 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 289 | 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 334 | 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 358 | 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() {
# Line 371 | Line 406 | void SimInfo::update() {
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  
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 –    
423      calcNdf();
424      calcNdfRaw();
425      calcNdfTrans();
# Line 392 | Line 428 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
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 432 | 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 506 | 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 528 | 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 573 | 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 601 | 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 657 | 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 669 | 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
# Line 719 | Line 768 | void SimInfo::setupCutoff() {
768          }
769  
770          if (globals_->haveRsw()) {
771 <            rsw_  = globals_->getRsw()
771 >            rsw_  = globals_->getRsw();
772          } else {
773              rsw_ = rcut_;
774          }
# Line 756 | 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