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

Comparing trunk/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 778 by mmeineke, Fri Sep 19 20:00:27 2003 UTC vs.
Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC

# Line 26 | Line 26 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
26   template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
27    T( theInfo, the_ff )
28   {
29 <  int i, j;
30 <  chi = 0.0;
31 <  integralOfChidt = 0.0;
32 <
33 <  for(i = 0; i < 3; i++)
34 <    for (j = 0; j < 3; j++)
29 >  
30 >  int i,j;
31 >  
32 >  for(i = 0; i < 3; i++){
33 >    for (j = 0; j < 3; j++){
34 >      
35        eta[i][j] = 0.0;
36 +      oldEta[i][j] = 0.0;
37 +    }
38 +  }
39 + }
40  
41 <  have_tau_thermostat = 0;
38 <  have_tau_barostat = 0;
39 <  have_target_temp = 0;
40 <  have_target_pressure = 0;
41 > template<typename T> NPTf<T>::~NPTf() {
42  
43 <  have_chi_tolerance = 0;
44 <  have_eta_tolerance = 0;
44 <  have_pos_iter_tolerance = 0;
43 >  // empty for now
44 > }
45  
46 <  oldPos = new double[3*nAtoms];
47 <  oldVel = new double[3*nAtoms];
48 <  oldJi = new double[3*nAtoms];
49 < #ifdef IS_MPI
50 <  Nparticles = mpiSim->getTotAtoms();
51 < #else
52 <  Nparticles = theInfo->n_atoms;
53 < #endif
54 <
46 > template<typename T> void NPTf<T>::resetIntegrator() {
47 >  
48 >  int i, j;
49 >  
50 >  for(i = 0; i < 3; i++)
51 >    for (j = 0; j < 3; j++)
52 >      eta[i][j] = 0.0;
53 >  
54 >  T::resetIntegrator();
55   }
56  
57 < template<typename T> NPTf<T>::~NPTf() {
58 <  delete[] oldPos;
59 <  delete[] oldVel;
60 <  delete[] oldJi;
57 > template<typename T> void NPTf<T>::evolveEtaA() {
58 >  
59 >  int i, j;
60 >  
61 >  for(i = 0; i < 3; i ++){
62 >    for(j = 0; j < 3; j++){
63 >      if( i == j)
64 >        eta[i][j] += dt2 *  instaVol *
65 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
66 >      else
67 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
68 >    }
69 >  }
70 >  
71 >  for(i = 0; i < 3; i++)
72 >    for (j = 0; j < 3; j++)
73 >      oldEta[i][j] = eta[i][j];
74   }
75  
76 < template<typename T> void NPTf<T>::moveA() {
76 > template<typename T> void NPTf<T>::evolveEtaB() {
77 >  
78 >  int i,j;
79  
80 <  // new version of NPTf
81 <  int i, j, k;
82 <  DirectionalAtom* dAtom;
68 <  double Tb[3], ji[3];
80 >  for(i = 0; i < 3; i++)
81 >    for (j = 0; j < 3; j++)
82 >      prevEta[i][j] = eta[i][j];
83  
84 <  double mass;
85 <  double vel[3], pos[3], frc[3];
84 >  for(i = 0; i < 3; i ++){
85 >    for(j = 0; j < 3; j++){
86 >      if( i == j) {
87 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
88 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
89 >      } else {
90 >        eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
91 >      }
92 >    }
93 >  }
94 > }
95  
96 <  double rj[3];
97 <  double instaTemp, instaPress, instaVol;
98 <  double tt2, tb2;
76 <  double sc[3];
77 <  double eta2ij;
78 <  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
79 <  double bigScale, smallScale, offDiagMax;
80 <  double COM[3];
96 > template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
97 >  int i,j;
98 >  double vScale[3][3];
99  
82 –  tt2 = tauThermostat * tauThermostat;
83 –  tb2 = tauBarostat * tauBarostat;
84 –
85 –  instaTemp = tStats->getTemperature();
86 –  tStats->getPressureTensor(press);
87 –  instaVol = tStats->getVolume();
88 –  
89 –  tStats->getCOM(COM);
90 –
91 –  //calculate scale factor of veloity
100    for (i = 0; i < 3; i++ ) {
101      for (j = 0; j < 3; j++ ) {
102        vScale[i][j] = eta[i][j];
# Line 99 | Line 107 | template<typename T> void NPTf<T>::moveA() {
107      }
108    }
109    
110 <  //evolve velocity half step
111 <  for( i=0; i<nAtoms; i++ ){
110 >  info->matVecMul3( vScale, vel, sc );
111 > }
112  
113 <    atoms[i]->getVel( vel );
114 <    atoms[i]->getFrc( frc );
113 > template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
114 >  int i,j;
115 >  double myVel[3];
116 >  double vScale[3][3];
117  
118 <    mass = atoms[i]->getMass();
119 <    
120 <    info->matVecMul3( vScale, vel, sc );
111 <
112 <    for (j=0; j < 3; j++) {
113 <      // velocity half step
114 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
115 <    }
116 <
117 <    atoms[i]->setVel( vel );
118 <  
119 <    if( atoms[i]->isDirectional() ){
120 <
121 <      dAtom = (DirectionalAtom *)atoms[i];
122 <
123 <      // get and convert the torque to body frame
118 >  for (i = 0; i < 3; i++ ) {
119 >    for (j = 0; j < 3; j++ ) {
120 >      vScale[i][j] = eta[i][j];
121        
122 <      dAtom->getTrq( Tb );
123 <      dAtom->lab2Body( Tb );
124 <      
128 <      // get the angular momentum, and propagate a half step
129 <
130 <      dAtom->getJ( ji );
131 <
132 <      for (j=0; j < 3; j++)
133 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
134 <      
135 <      this->rotationPropagation( dAtom, ji );
136 <  
137 <      dAtom->setJ( ji );
138 <    }    
139 <  }
140 <
141 <  // advance chi half step
142 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
143 <
144 <  // calculate the integral of chidt
145 <  integralOfChidt += dt2*chi;
146 <
147 <  // advance eta half step
148 <
149 <  for(i = 0; i < 3; i ++)
150 <    for(j = 0; j < 3; j++){
151 <      if( i == j)
152 <        eta[i][j] += dt2 *  instaVol *
153 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
154 <      else
155 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
122 >      if (i == j) {
123 >        vScale[i][j] += chi;          
124 >      }              
125      }
157 –    
158 –  //save the old positions
159 –  for(i = 0; i < nAtoms; i++){
160 –    atoms[i]->getPos(pos);
161 –    for(j = 0; j < 3; j++)
162 –      oldPos[i*3 + j] = pos[j];
126    }
127    
128 <  //the first estimation of r(t+dt) is equal to  r(t)
129 <    
167 <  for(k = 0; k < 4; k ++){
128 >  for (j = 0; j < 3; j++)
129 >    myVel[j] = oldVel[3*index + j];
130  
131 <    for(i =0 ; i < nAtoms; i++){
131 >  info->matVecMul3( vScale, myVel, sc );
132 > }
133  
134 <      atoms[i]->getVel(vel);
135 <      atoms[i]->getPos(pos);
134 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
135 >                                               int index, double sc[3]){
136 >  int j;
137 >  double rj[3];
138  
139 <      for(j = 0; j < 3; j++)
140 <        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
176 <      
177 <      info->matVecMul3( eta, rj, sc );
178 <      
179 <      for(j = 0; j < 3; j++)
180 <        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
139 >  for(j=0; j<3; j++)
140 >    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
141  
142 <      atoms[i]->setPos( pos );
142 >  info->matVecMul3( eta, rj, sc );
143 > }
144  
145 <    }
145 > template<typename T> void NPTf<T>::scaleSimBox( void ){
146  
147 <    if (nConstrained) {
148 <      constrainA();
149 <    }
150 <  }  
147 >  int i,j,k;
148 >  double scaleMat[3][3];
149 >  double eta2ij;
150 >  double bigScale, smallScale, offDiagMax;
151 >  double hm[3][3], hmnew[3][3];
152 >  
153  
154 <
154 >
155    // Scale the box after all the positions have been moved:
156    
157    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
# Line 253 | Line 216 | template<typename T> void NPTf<T>::moveA() {
216      info->matMul3(hm, scaleMat, hmnew);
217      info->setBoxM(hmnew);
218    }
256 –  
219   }
220  
221 < template<typename T> void NPTf<T>::moveB( void ){
222 <
223 <  //new version of NPTf
262 <  int i, j, k;
263 <  DirectionalAtom* dAtom;
264 <  double Tb[3], ji[3];
265 <  double vel[3], myVel[3], frc[3];
266 <  double mass;
221 > template<typename T> bool NPTf<T>::etaConverged() {
222 >  int i;
223 >  double diffEta, sumEta;
224  
225 <  double instaTemp, instaPress, instaVol;
269 <  double tt2, tb2;
270 <  double sc[3];
271 <  double press[3][3], vScale[3][3];
272 <  double oldChi, prevChi;
273 <  double oldEta[3][3], prevEta[3][3], diffEta;
274 <  
275 <  tt2 = tauThermostat * tauThermostat;
276 <  tb2 = tauBarostat * tauBarostat;
277 <
278 <  // Set things up for the iteration:
279 <
280 <  oldChi = chi;
281 <  
225 >  sumEta = 0;
226    for(i = 0; i < 3; i++)
227 <    for(j = 0; j < 3; j++)
284 <      oldEta[i][j] = eta[i][j];
285 <
286 <  for( i=0; i<nAtoms; i++ ){
287 <
288 <    atoms[i]->getVel( vel );
289 <
290 <    for (j=0; j < 3; j++)
291 <      oldVel[3*i + j]  = vel[j];
292 <
293 <    if( atoms[i]->isDirectional() ){
294 <
295 <      dAtom = (DirectionalAtom *)atoms[i];
296 <
297 <      dAtom->getJ( ji );
298 <
299 <      for (j=0; j < 3; j++)
300 <        oldJi[3*i + j] = ji[j];
301 <
302 <    }
303 <  }
304 <
305 <  // do the iteration:
306 <
307 <  instaVol = tStats->getVolume();
227 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
228    
229 <  for (k=0; k < 4; k++) {
310 <    
311 <    instaTemp = tStats->getTemperature();
312 <    tStats->getPressureTensor(press);
313 <
314 <    // evolve chi another half step using the temperature at t + dt/2
315 <
316 <    prevChi = chi;
317 <    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
318 <    
319 <    for(i = 0; i < 3; i++)
320 <      for(j = 0; j < 3; j++)
321 <        prevEta[i][j] = eta[i][j];
322 <
323 <    //advance eta half step and calculate scale factor for velocity
324 <
325 <    for(i = 0; i < 3; i ++)
326 <      for(j = 0; j < 3; j++){
327 <        if( i == j) {
328 <          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
329 <            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
330 <          vScale[i][j] = eta[i][j] + chi;
331 <        } else {
332 <          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
333 <          vScale[i][j] = eta[i][j];
334 <        }
335 <      }  
336 <    
337 <    for( i=0; i<nAtoms; i++ ){
338 <
339 <      atoms[i]->getFrc( frc );
340 <      atoms[i]->getVel(vel);
341 <      
342 <      mass = atoms[i]->getMass();
343 <    
344 <      for (j = 0; j < 3; j++)
345 <        myVel[j] = oldVel[3*i + j];
346 <      
347 <      info->matVecMul3( vScale, myVel, sc );
348 <      
349 <      // velocity half step
350 <      for (j=0; j < 3; j++) {
351 <        // velocity half step  (use chi from previous step here):
352 <        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
353 <      }
354 <      
355 <      atoms[i]->setVel( vel );
356 <      
357 <      if( atoms[i]->isDirectional() ){
358 <
359 <        dAtom = (DirectionalAtom *)atoms[i];
229 >  diffEta = sqrt( sumEta / 3.0 );
230    
231 <        // get and convert the torque to body frame      
362 <  
363 <        dAtom->getTrq( Tb );
364 <        dAtom->lab2Body( Tb );      
365 <            
366 <        for (j=0; j < 3; j++)
367 <          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
368 <      
369 <          dAtom->setJ( ji );
370 <      }
371 <    }
372 <
373 <    if (nConstrained) {
374 <      constrainB();
375 <    }
376 <    
377 <    diffEta = 0;
378 <    for(i = 0; i < 3; i++)
379 <      diffEta += pow(prevEta[i][i] - eta[i][i], 2);    
380 <    
381 <    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
382 <      break;
383 <  }
384 <
385 <  //calculate integral of chidt
386 <  integralOfChidt += dt2*chi;
387 <  
231 >  return ( diffEta <= etaTolerance );
232   }
233  
390 – template<typename T> void NPTf<T>::resetIntegrator() {
391 –  int i,j;
392 –  
393 –  chi = 0.0;
394 –
395 –  for(i = 0; i < 3; i++)
396 –    for (j = 0; j < 3; j++)
397 –      eta[i][j] = 0.0;
398 –
399 – }
400 –
401 – template<typename T> int NPTf<T>::readyCheck() {
402 –
403 –  //check parent's readyCheck() first
404 –  if (T::readyCheck() == -1)
405 –    return -1;
406 –
407 –  // First check to see if we have a target temperature.
408 –  // Not having one is fatal.
409 –  
410 –  if (!have_target_temp) {
411 –    sprintf( painCave.errMsg,
412 –             "NPTf error: You can't use the NPTf integrator\n"
413 –             "   without a targetTemp!\n"
414 –             );
415 –    painCave.isFatal = 1;
416 –    simError();
417 –    return -1;
418 –  }
419 –
420 –  if (!have_target_pressure) {
421 –    sprintf( painCave.errMsg,
422 –             "NPTf error: You can't use the NPTf integrator\n"
423 –             "   without a targetPressure!\n"
424 –             );
425 –    painCave.isFatal = 1;
426 –    simError();
427 –    return -1;
428 –  }
429 –  
430 –  // We must set tauThermostat.
431 –  
432 –  if (!have_tau_thermostat) {
433 –    sprintf( painCave.errMsg,
434 –             "NPTf error: If you use the NPTf\n"
435 –             "   integrator, you must set tauThermostat.\n");
436 –    painCave.isFatal = 1;
437 –    simError();
438 –    return -1;
439 –  }    
440 –
441 –  // We must set tauBarostat.
442 –  
443 –  if (!have_tau_barostat) {
444 –    sprintf( painCave.errMsg,
445 –             "NPTf error: If you use the NPTf\n"
446 –             "   integrator, you must set tauBarostat.\n");
447 –    painCave.isFatal = 1;
448 –    simError();
449 –    return -1;
450 –  }    
451 –
452 –  
453 –  // We need NkBT a lot, so just set it here: This is the RAW number
454 –  // of particles, so no subtraction or addition of constraints or
455 –  // orientational degrees of freedom:
456 –  
457 –  NkBT = (double)Nparticles * kB * targetTemp;
458 –  
459 –  // fkBT is used because the thermostat operates on more degrees of freedom
460 –  // than the barostat (when there are particles with orientational degrees
461 –  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
462 –  
463 –  fkBT = (double)info->ndf * kB * targetTemp;
464 –
465 –  return 1;
466 – }
467 –
234   template<typename T> double NPTf<T>::getConservedQuantity(void){
235 <
235 >  
236    double conservedQuantity;
237    double Energy;
238    double thermostat_kinetic;
# Line 478 | Line 244 | template<typename T> double NPTf<T>::getConservedQuant
244  
245    Energy = tStats->getTotalE();
246  
247 <  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
247 >  thermostat_kinetic = fkBT* tt2 * chi * chi /
248      (2.0 * eConvert);
249  
250    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 487 | Line 253 | template<typename T> double NPTf<T>::getConservedQuant
253    info->matMul3(a, eta, b);
254    trEta = info->matTrace3(b);
255  
256 <  barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta /
256 >  barostat_kinetic = NkBT * tb2 * trEta /
257      (2.0 * eConvert);
258    
259    barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
# Line 496 | Line 262 | template<typename T> double NPTf<T>::getConservedQuant
262    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
263      barostat_kinetic + barostat_potential;
264    
265 <  cout.width(8);
266 <  cout.precision(8);
265 > //   cout.width(8);
266 > //   cout.precision(8);
267  
268 <  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269 <      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270 <      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
268 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270 > //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
271  
272    return conservedQuantity;
273 +  
274   }

Diff Legend

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