ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ForceFields.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/ForceFields.cpp (file contents):
Revision 1158 by tim, Tue May 11 21:14:26 2004 UTC vs.
Revision 1187 by chrisfen, Sat May 22 18:16:18 2004 UTC

# Line 81 | Line 81 | void ForceFields::doForces( int calcPot, int calcStres
81    double* u_l;
82    double* rc;
83    double* massRatio;
84 +  double factor;
85    SimState* config;
86  
87    Molecule* myMols;
# Line 92 | Line 93 | void ForceFields::doForces( int calcPot, int calcStres
93    int numCutoffGroups;
94    CutoffGroup* myCutoffGroup;
95    vector<CutoffGroup*>::iterator iterCutoff;
95 –  Atom* cutoffAtom;
96 –  vector<Atom*>::iterator iterAtom;  
96    double com[3];
97 <  double tempPos[3];
99 <  int atomIndex;
97 >  vector<double> rcGroup;
98    
99    short int passedCalcPot = (short int)calcPot;
100    short int passedCalcStress = (short int)calcStress;
# Line 131 | Line 129 | void ForceFields::doForces( int calcPot, int calcStres
129    u_l = config->getUlArray();
130  
131    if(entry_plug->haveCutoffGroups){
134 –    //if
132      myMols = entry_plug->molecules;
133      numMol = entry_plug->n_mol;
134      for(int i  = 0; i < numMol; i++){
138 –      numAtom = myMols[i].getNAtoms();
139 –      myAtoms = myMols[i].getMyAtoms();
140 –
141 –      
142 –      for(int j = 0; j < numAtom; j++){
143 – #ifdef IS_MPI
144 –        atomIndex = myAtoms[j]->getGlobalIndex();
145 – #else
146 –        atomIndex = myAtoms[j]->getIndex();
147 – #endif
148 –
149 –        if(myMols[i].belongToCutoffGroup(atomIndex))
150 –          continue;
151 –        else{
152 –          myAtoms[j]->getPos(tempPos);
153 –          myAtoms[j]->setRc(tempPos);
154 –        }
155 –          
156 –      }
135          
136        numCutoffGroups = myMols[i].getNCutoffGroups();
137        for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL;
138                                                      myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
139          //get center of mass of the cutoff group
140 <        myCutoffGroup->getCOM(com);
140 >        myCutoffGroup->getCOM(com);
141  
142 <        for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL;
143 <                                             cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
144 <          cutoffAtom->setRc(com);
145 <        }  
168 <                                
142 >        rcGroup.push_back(com[0]);
143 >        rcGroup.push_back(com[1]);
144 >        rcGroup.push_back(com[2]);
145 >        
146        }// end for(myCutoffGroup)
147        
148      }//end for(int i = 0)
149  
150 <    rc = config->getRcArray();
150 >    rc = &rcGroup[0];
151    }
152    else{
153      // center of mass of the group is the same as position of the atom  if cutoff group does not exist
# Line 221 | Line 198 | void ForceFields::doForces( int calcPot, int calcStres
198    }
199  
200  
201 +  if (entry_plug->useThermInt) {
202 +
203 +    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
204 +    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
205 +      for (j=0; j< 3; j++)
206 +        frc[3*i + j] *= factor;
207 +      if (entry_plug->integrableObjects[i]->isDirectional()) {
208 +        for (j=0; j< 3; j++)
209 +          trq[3*i + j] *= factor;
210 +      }
211 +    }
212 +    entry_plug->vRaw = entry_plug->lrPot;
213 +    entry_plug->lrPot *= factor;
214 +    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
215 +    entry_plug->vHarm = entry_plug->restraint->getVharm();
216 +  }
217 +
218   #ifdef IS_MPI
219    sprintf( checkPointMsg,
220             "returned from the force calculation.\n" );
# Line 252 | Line 246 | void ForceFields::initFortran(int ljMixPolicy, int use
246   #endif // is_mpi
247    
248   }
249 +
250 +
251 + void ForceFields::initRestraints(){
252 +  int i;
253 +  // store the initial info.
254 +  // set the omega values to zero
255 +  for (i=0; i<entry_plug->integrableObjects.size(); i++)
256 +    entry_plug->integrableObjects[i]->setZangle( 0.0 );
257 +
258 +  entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects);
259 +
260 + }
261 +
262 + void ForceFields::dumpzAngle(){
263 +
264 +  // store the initial info.
265 +  entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects);
266 +
267 + }

Diff Legend

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