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 1258 by chrisfen, Thu Jun 10 17:09:16 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 214 | Line 191 | void ForceFields::doForces( int calcPot, int calcStres
191               "Error returned from the fortran force calculation.\n" );
192      painCave.isFatal = 1;
193      simError();
194 +  }
195 +
196 +  if (entry_plug->useSolidThermInt && !entry_plug->useLiquidThermInt) {
197 +    
198 +    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
199 +    for (i=0; i < entry_plug->n_atoms; i++) {
200 +      for (j=0; j< 3; j++)
201 +        frc[3*i + j] *= factor;
202 +      if (entry_plug->atoms[i]->isDirectional()) {
203 +        for (j=0; j< 3; j++)
204 +          trq[3*i + j] *= factor;
205 +      }
206 +    }
207 +    entry_plug->vRaw = entry_plug->lrPot;
208 +    entry_plug->lrPot *= factor;
209 +    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
210 +    entry_plug->vHarm = entry_plug->restraint->getVharm();
211    }
212 +  
213 +  if (entry_plug->useLiquidThermInt) {
214 +    
215 +    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
216 +    for (i=0; i < entry_plug->n_atoms; i++) {
217 +      for (j=0; j< 3; j++)
218 +        frc[3*i + j] *= factor;
219 +      if (entry_plug->atoms[i]->isDirectional()) {
220 +        for (j=0; j< 3; j++)
221 +          trq[3*i + j] *= factor;
222 +      }
223 +    }
224 +    entry_plug->vRaw = entry_plug->lrPot;
225 +    entry_plug->lrPot *= factor;
226 +  }
227  
228    for(i=0; i<entry_plug->n_mol; i++ ){
229      entry_plug->molecules[i].atoms2rigidBodies();
# Line 252 | Line 261 | void ForceFields::initFortran(int ljMixPolicy, int use
261   #endif // is_mpi
262    
263   }
264 +
265 +
266 + void ForceFields::initRestraints(){
267 +  int i;
268 +  // store the initial info.
269 +  // set the omega values to zero
270 +  for (i=0; i<entry_plug->integrableObjects.size(); i++)
271 +    entry_plug->integrableObjects[i]->setZangle( 0.0 );
272 +
273 +  entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects);
274 +
275 + }
276 +
277 + void ForceFields::dumpzAngle(){
278 +
279 +  // store the initial info.
280 +  entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects);
281 +
282 + }

Diff Legend

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