ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/brains/SimInfo.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/brains/SimInfo.cpp (file contents):
Revision 1735 by tim, Fri Nov 12 17:40:03 2004 UTC vs.
Revision 1804 by tim, Tue Nov 30 19:58:25 2004 UTC

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

Diff Legend

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