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 1157 by tim, Tue May 11 20:33:41 2004 UTC vs.
Revision 1212 by chrisfen, Tue Jun 1 17:15:43 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->beginAtom(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 218 | Line 195 | void ForceFields::doForces( int calcPot, int calcStres
195  
196    for(i=0; i<entry_plug->n_mol; i++ ){
197      entry_plug->molecules[i].atoms2rigidBodies();
198 +  }
199 +
200 +
201 +  if (entry_plug->useSolidThermInt && !entry_plug->useLiquidThermInt) {
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 +  if (entry_plug->useLiquidThermInt) {
219  
220 +    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
221 +    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
222 +      for (j=0; j< 3; j++)
223 +        frc[3*i + j] *= factor;
224 +      if (entry_plug->integrableObjects[i]->isDirectional()) {
225 +        for (j=0; j< 3; j++)
226 +          trq[3*i + j] *= factor;
227 +      }
228 +    }
229 +    entry_plug->vRaw = entry_plug->lrPot;
230 +    entry_plug->lrPot *= factor;
231 +    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
232 +    entry_plug->vHarm = entry_plug->restraint->getVharm();
233 +  }
234 +
235   #ifdef IS_MPI
236    sprintf( checkPointMsg,
237             "returned from the force calculation.\n" );
# Line 252 | Line 263 | void ForceFields::initFortran(int ljMixPolicy, int use
263   #endif // is_mpi
264    
265   }
266 +
267 +
268 + void ForceFields::initRestraints(){
269 +  int i;
270 +  // store the initial info.
271 +  // set the omega values to zero
272 +  for (i=0; i<entry_plug->integrableObjects.size(); i++)
273 +    entry_plug->integrableObjects[i]->setZangle( 0.0 );
274 +
275 +  entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects);
276 +
277 + }
278 +
279 + void ForceFields::dumpzAngle(){
280 +
281 +  // store the initial info.
282 +  entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects);
283 +
284 + }

Diff Legend

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