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 1735 by tim, Fri Nov 12 17:40:03 2004 UTC vs.
Revision 1739 by tim, Mon Nov 15 18:02:15 2004 UTC

# Line 41 | Line 41 | SimInfo::SimInfo(const std::vector<std::pair<MoleculeS
41                                  ForceField* ff, Globals* globals) :
42                                  forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
43                                  nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
44 <                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) {
44 >                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL),
45 >                                fortranInitialized_(false) {
46  
47      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
48      MoleculeStamp* molStamp;
49      int nMolWithSameStamp;
50 +    int nCutoffAtoms; // number of atoms belong to cutoff groups
51 +    int nGroups;          //total cutoff groups defined in meta-data file
52      CutoffGroupStamp* cgStamp;
53 <    int nAtomsIngroups;
54 <    int nCutoffGroupsInStamp;    
55 <
53 >    int nAtomsInGroups;
54 >    int nCutoffGroupsInStamp;
55 >    
56 >    RigidBodyStamp* rbStamp;
57 >    int nAtomsInRigidBodies;
58 >    int nRigidBodiesInStamp;
59 >    int nRigidAtoms;
60 >    int nRigidBodies;
61 >        
62      nGlobalAtoms_ =  0;
56 –    ngroups = 0;
63      
64 +    nGroups = 0;
65 +    nCutoffAtoms = 0;
66 +
67 +    nRigidBodies
68 +    nRigidBodies = 0;
69 +    
70      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
71          molStamp = i->first;
72          nMolWithSameStamp = i->second;
73          
74          addMoleculeStamp(molStamp, nMolWithSameStamp);
75 <        
75 >
76 >        //calculate atoms in molecules
77          nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
78 <        
79 <        nAtomsIngroups = 0;
78 >
79 >
80 >        //calculate atoms in cutoff groups
81 >        nAtomsInGroups = 0;
82          nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
83          
84          for (int j=0; j < nCutoffGroupsInStamp; j++) {
85              cgStamp = molStamp->getCutoffGroup(j);
86 <            nAtomsIngroups += cgStamp->getNMembers();
86 >            nAtomsInGroups += cgStamp->getNMembers();
87          }
88  
89 <        ngroups += *nMolWithSameStamp;
90 <        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
89 >        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
90 >        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
91 >
92 >        //calculate atoms in rigid bodies
93 >        nAtomsInRigidBodies = 0;
94 >        nRigidBodiesInStamp = molStamp->getNCutoffGroups();
95 >        
96 >        for (int j=0; j < nRigidBodiesInStamp; j++) {
97 >            rbStamp = molStamp->getRigidBody(j);
98 >            nRigidBodiesInStamp += rbStamp->getNMembers();
99 >        }
100 >
101 >        nRigidBodies += nRigidBodiesInStamp * nMolWithSameStamp;
102 >        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
103 >        
104      }
105  
106      //every free atom (atom does not belong to cutoff groups) is a cutoff group
107      //therefore the total number of cutoff groups in the system is equal to
108      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
109      //file plus the number of cutoff groups defined in meta-data file
110 <    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups;
110 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
111 >
112 >    //every free atom (atom does not belong to rigid bodies) is a rigid body
113 >    //therefore the total number of cutoff groups in the system is equal to
114 >    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
115 >    //file plus the number of  rigid bodies defined in meta-data file
116 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nRigidBodies;
117  
118      //initialize globalGroupMembership_, every element of this array will be 0
119      globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
# Line 159 | Line 193 | Molecule* SimInfo::beginMolecule(MoleculeIterator& i)
193          
194   Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
195      i = molecules_.begin();
196 <    return i == molecules_.end() ? NULL : *i;
196 >    return i == molecules_.end() ? NULL : i->second;
197   }    
198  
199   Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
200      ++i;
201 <    return i == molecules_.end() ? NULL : *i;    
201 >    return i == molecules_.end() ? NULL : i->second;    
202   }
203  
204  
# Line 362 | Line 396 | void SimInfo::update() {
396  
397   void SimInfo::update() {
398  
399 +    setupSimType();
400 +
401   #ifdef IS_MPI
402      setupFortranParallel();
403   #endif
404  
405      setupFortranSim();
406 +
407 +    setupCutoff();
408  
409 +    //notify fortran whether reaction field is used or not. It is deprecated now
410 +    //int isError = 0;
411 +    //initFortranFF( &useReactionField, &isError );
412 +
413 +    //if(isError){
414 +    //    sprintf( painCave.errMsg,
415 +    //    "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" );
416 +    //    painCave.isFatal = 1;
417 +    //    simError();
418 +    //}
419 +    
420      calcNdf();
421      calcNdfRaw();
422      calcNdfTrans();
423 +
424 +    fortranInitialized_ = true;
425   }
426  
427   std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
# Line 657 | Line 708 | double SimInfo::calcMaxCutoffRadius() {
708   #endif
709  
710      return maxCutoffRadius;
711 + }
712 +
713 + void SimInfo::setupCutoff() {
714 +    double rcut_;  //cutoff radius
715 +    double rsw_; //switching radius
716 +    
717 +    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
718 +        
719 +        if (!globals_->haveRcut()){
720 +            sprintf(painCave.errMsg,
721 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
722 +                "\tOOPSE will use a default value of 15.0 angstroms"
723 +                "\tfor the cutoffRadius.\n");
724 +            painCave.isFatal = 0;
725 +            simError();
726 +            rcut_ = 15.0;
727 +        } else{
728 +            rcut_ = globals_->getRcut();
729 +        }
730 +
731 +        if (!globals_->haveRsw()){
732 +            sprintf(painCave.errMsg,
733 +                "SimCreator Warning: No value was set for switchingRadius.\n"
734 +                "\tOOPSE will use a default value of\n"
735 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
736 +            painCave.isFatal = 0;
737 +            simError();
738 +            rsw_ = 0.95 * rcut_;
739 +        } else{
740 +            rsw_ = globals_->getRsw();
741 +        }
742 +
743 +    } else {
744 +        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
745 +        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
746 +        
747 +        if (globals_->haveRcut()) {
748 +            rcut_ = globals_->getRcut();
749 +        } else {
750 +            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
751 +            rcut_ = calcMaxCutoffRadius();
752 +        }
753 +
754 +        if (globals_->haveRsw()) {
755 +            rsw_  = globals_->getRsw()
756 +        } else {
757 +            rsw_ = rcut_;
758 +        }
759 +    
760 +    }
761 +        
762 +    double rnblist = rcut_ + 1; // skin of neighbor list
763 +
764 +    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
765 +    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
766   }
767  
768   void SimInfo::addProperty(GenericData* genData) {

Diff Legend

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