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

Comparing trunk/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 772 by gezelter, Fri Sep 19 16:01:07 2003 UTC vs.
Revision 782 by mmeineke, Tue Sep 23 20:34:31 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];
69 <  double A[3][3], I[3][3];
70 <  double angle, mass;
71 <  double vel[3], pos[3], frc[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 rj[3];
85 <  double instaTemp, instaPress, instaVol;
86 <  double tt2, tb2;
87 <  double sc[3];
88 <  double eta2ij;
89 <  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
90 <  double bigScale, smallScale, offDiagMax;
91 <  double COM[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 <  tt2 = tauThermostat * tauThermostat;
97 <  tb2 = tauBarostat * tauBarostat;
96 > template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
97 >  int i,j;
98 >  double vScale[3][3];
99  
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 <      // use the angular velocities to propagate the rotation matrix a
136 <      // full time step
137 <
138 <      dAtom->getA(A);
139 <      dAtom->getI(I);
140 <    
141 <      // rotate about the x-axis      
142 <      angle = dt2 * ji[0] / I[0][0];
143 <      this->rotate( 1, 2, angle, ji, A );
144 <
145 <      // rotate about the y-axis
146 <      angle = dt2 * ji[1] / I[1][1];
147 <      this->rotate( 2, 0, angle, ji, A );
148 <      
149 <      // rotate about the z-axis
150 <      angle = dt * ji[2] / I[2][2];
151 <      this->rotate( 0, 1, angle, ji, A);
152 <      
153 <      // rotate about the y-axis
154 <      angle = dt2 * ji[1] / I[1][1];
155 <      this->rotate( 2, 0, angle, ji, A );
156 <      
157 <       // rotate about the x-axis
158 <      angle = dt2 * ji[0] / I[0][0];
159 <      this->rotate( 1, 2, angle, ji, A );
160 <      
161 <      dAtom->setJ( ji );
162 <      dAtom->setA( A  );    
163 <    }    
164 <  }
165 <
166 <  // advance chi half step
167 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
168 <
169 <  // calculate the integral of chidt
170 <  integralOfChidt += dt2*chi;
171 <
172 <  // advance eta half step
173 <
174 <  for(i = 0; i < 3; i ++)
175 <    for(j = 0; j < 3; j++){
176 <      if( i == j)
177 <        eta[i][j] += dt2 *  instaVol *
178 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
179 <      else
180 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
122 >      if (i == j) {
123 >        vScale[i][j] += chi;          
124 >      }              
125      }
182 –    
183 –  //save the old positions
184 –  for(i = 0; i < nAtoms; i++){
185 –    atoms[i]->getPos(pos);
186 –    for(j = 0; j < 3; j++)
187 –      oldPos[i*3 + j] = pos[j];
126    }
127    
128 <  //the first estimation of r(t+dt) is equal to  r(t)
129 <    
192 <  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];
201 <      
202 <      info->matVecMul3( eta, rj, sc );
203 <      
204 <      for(j = 0; j < 3; j++)
205 <        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 277 | Line 215 | template<typename T> void NPTf<T>::moveA() {
215      info->getBoxM(hm);
216      info->matMul3(hm, scaleMat, hmnew);
217      info->setBoxM(hmnew);
280 –  }
281 –  
282 – }
283 –
284 – template<typename T> void NPTf<T>::moveB( void ){
285 –
286 –  //new version of NPTf
287 –  int i, j, k;
288 –  DirectionalAtom* dAtom;
289 –  double Tb[3], ji[3];
290 –  double vel[3], myVel[3], frc[3];
291 –  double mass;
292 –
293 –  double instaTemp, instaPress, instaVol;
294 –  double tt2, tb2;
295 –  double sc[3];
296 –  double press[3][3], vScale[3][3];
297 –  double oldChi, prevChi;
298 –  double oldEta[3][3], prevEta[3][3], diffEta;
299 –  
300 –  tt2 = tauThermostat * tauThermostat;
301 –  tb2 = tauBarostat * tauBarostat;
302 –
303 –  // Set things up for the iteration:
304 –
305 –  oldChi = chi;
306 –  
307 –  for(i = 0; i < 3; i++)
308 –    for(j = 0; j < 3; j++)
309 –      oldEta[i][j] = eta[i][j];
310 –
311 –  for( i=0; i<nAtoms; i++ ){
312 –
313 –    atoms[i]->getVel( vel );
314 –
315 –    for (j=0; j < 3; j++)
316 –      oldVel[3*i + j]  = vel[j];
317 –
318 –    if( atoms[i]->isDirectional() ){
319 –
320 –      dAtom = (DirectionalAtom *)atoms[i];
321 –
322 –      dAtom->getJ( ji );
323 –
324 –      for (j=0; j < 3; j++)
325 –        oldJi[3*i + j] = ji[j];
326 –
327 –    }
218    }
329 –
330 –  // do the iteration:
331 –
332 –  instaVol = tStats->getVolume();
333 –  
334 –  for (k=0; k < 4; k++) {
335 –    
336 –    instaTemp = tStats->getTemperature();
337 –    tStats->getPressureTensor(press);
338 –
339 –    // evolve chi another half step using the temperature at t + dt/2
340 –
341 –    prevChi = chi;
342 –    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
343 –    
344 –    for(i = 0; i < 3; i++)
345 –      for(j = 0; j < 3; j++)
346 –        prevEta[i][j] = eta[i][j];
347 –
348 –    //advance eta half step and calculate scale factor for velocity
349 –
350 –    for(i = 0; i < 3; i ++)
351 –      for(j = 0; j < 3; j++){
352 –        if( i == j) {
353 –          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
354 –            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
355 –          vScale[i][j] = eta[i][j] + chi;
356 –        } else {
357 –          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
358 –          vScale[i][j] = eta[i][j];
359 –        }
360 –      }  
361 –    
362 –    for( i=0; i<nAtoms; i++ ){
363 –
364 –      atoms[i]->getFrc( frc );
365 –      atoms[i]->getVel(vel);
366 –      
367 –      mass = atoms[i]->getMass();
368 –    
369 –      for (j = 0; j < 3; j++)
370 –        myVel[j] = oldVel[3*i + j];
371 –      
372 –      info->matVecMul3( vScale, myVel, sc );
373 –      
374 –      // velocity half step
375 –      for (j=0; j < 3; j++) {
376 –        // velocity half step  (use chi from previous step here):
377 –        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
378 –      }
379 –      
380 –      atoms[i]->setVel( vel );
381 –      
382 –      if( atoms[i]->isDirectional() ){
383 –
384 –        dAtom = (DirectionalAtom *)atoms[i];
385 –  
386 –        // get and convert the torque to body frame      
387 –  
388 –        dAtom->getTrq( Tb );
389 –        dAtom->lab2Body( Tb );      
390 –            
391 –        for (j=0; j < 3; j++)
392 –          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
393 –      
394 –          dAtom->setJ( ji );
395 –      }
396 –    }
397 –
398 –    if (nConstrained) {
399 –      constrainB();
400 –    }
401 –    
402 –    diffEta = 0;
403 –    for(i = 0; i < 3; i++)
404 –      diffEta += pow(prevEta[i][i] - eta[i][i], 2);    
405 –    
406 –    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
407 –      break;
408 –  }
409 –
410 –  //calculate integral of chidt
411 –  integralOfChidt += dt2*chi;
412 –  
219   }
220  
221 < template<typename T> void NPTf<T>::resetIntegrator() {
222 <  int i,j;
223 <  
418 <  chi = 0.0;
221 > template<typename T> bool NPTf<T>::etaConverged() {
222 >  int i;
223 >  double diffEta, sumEta;
224  
225 <  for(i = 0; i < 3; i++)
226 <    for (j = 0; j < 3; j++)
227 <      eta[i][j] = 0.0;
423 <
424 < }
425 <
426 < template<typename T> int NPTf<T>::readyCheck() {
427 <
428 <  //check parent's readyCheck() first
429 <  if (T::readyCheck() == -1)
430 <    return -1;
431 <
432 <  // First check to see if we have a target temperature.
433 <  // Not having one is fatal.
225 >  sumEta = 0;
226 >  for(i = 0; i < 3; i++)
227 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
228    
229 <  if (!have_target_temp) {
436 <    sprintf( painCave.errMsg,
437 <             "NPTf error: You can't use the NPTf integrator\n"
438 <             "   without a targetTemp!\n"
439 <             );
440 <    painCave.isFatal = 1;
441 <    simError();
442 <    return -1;
443 <  }
444 <
445 <  if (!have_target_pressure) {
446 <    sprintf( painCave.errMsg,
447 <             "NPTf error: You can't use the NPTf integrator\n"
448 <             "   without a targetPressure!\n"
449 <             );
450 <    painCave.isFatal = 1;
451 <    simError();
452 <    return -1;
453 <  }
229 >  diffEta = sqrt( sumEta / 3.0 );
230    
231 <  // We must set tauThermostat.
456 <  
457 <  if (!have_tau_thermostat) {
458 <    sprintf( painCave.errMsg,
459 <             "NPTf error: If you use the NPTf\n"
460 <             "   integrator, you must set tauThermostat.\n");
461 <    painCave.isFatal = 1;
462 <    simError();
463 <    return -1;
464 <  }    
465 <
466 <  // We must set tauBarostat.
467 <  
468 <  if (!have_tau_barostat) {
469 <    sprintf( painCave.errMsg,
470 <             "NPTf error: If you use the NPTf\n"
471 <             "   integrator, you must set tauBarostat.\n");
472 <    painCave.isFatal = 1;
473 <    simError();
474 <    return -1;
475 <  }    
476 <
477 <  
478 <  // We need NkBT a lot, so just set it here: This is the RAW number
479 <  // of particles, so no subtraction or addition of constraints or
480 <  // orientational degrees of freedom:
481 <  
482 <  NkBT = (double)Nparticles * kB * targetTemp;
483 <  
484 <  // fkBT is used because the thermostat operates on more degrees of freedom
485 <  // than the barostat (when there are particles with orientational degrees
486 <  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
487 <  
488 <  fkBT = (double)info->ndf * kB * targetTemp;
489 <
490 <  return 1;
231 >  return ( diffEta <= etaTolerance );
232   }
233  
234   template<typename T> double NPTf<T>::getConservedQuantity(void){
235 <
235 >  
236    double conservedQuantity;
237 <  double Energy;
237 >  double totalEnergy;
238    double thermostat_kinetic;
239    double thermostat_potential;
240    double barostat_kinetic;
# Line 501 | Line 242 | template<typename T> double NPTf<T>::getConservedQuant
242    double trEta;
243    double a[3][3], b[3][3];
244  
245 <  Energy = tStats->getTotalE();
245 >  totalEnergy = 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 512 | 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) /
260      eConvert;
261  
262 <  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
262 >  conservedQuantity = totalEnergy + 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