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

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 586 by mmeineke, Wed Jul 9 22:14:06 2003 UTC vs.
Revision 769 by tim, Fri Sep 19 14:22:29 2003 UTC

# Line 9 | Line 9
9   #include "Integrator.hpp"
10   #include "simError.h"
11  
12 + #ifdef IS_MPI
13 + #include "mpiSimulation.hpp"
14 + #endif
15  
16 +
17   // Basic isotropic thermostating and barostating via the Melchionna
18   // modification of the Hoover algorithm:
19   //
# Line 20 | Line 24 | NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
24   //
25   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26  
27 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
28 <  Integrator( theInfo, the_ff )
27 > template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
28 >  T( theInfo, the_ff )
29   {
30    chi = 0.0;
31    eta = 0.0;
32 +  integralOfChidt = 0.0;
33    have_tau_thermostat = 0;
34    have_tau_barostat = 0;
35    have_target_temp = 0;
36    have_target_pressure = 0;
37 +  have_chi_tolerance = 0;
38 +  have_eta_tolerance = 0;
39 +  have_pos_iter_tolerance = 0;
40 +
41 +  oldPos = new double[3*nAtoms];
42 +  oldVel = new double[3*nAtoms];
43 +  oldJi = new double[3*nAtoms];
44 + #ifdef IS_MPI
45 +  Nparticles = mpiSim->getTotAtoms();
46 + #else
47 +  Nparticles = theInfo->n_atoms;
48 + #endif
49 +
50   }
51  
52 < void NPTi::moveA() {
53 <  
54 <  int i,j,k;
55 <  int atomIndex, aMatIndex;
52 > template<typename T> NPTi<T>::~NPTi() {
53 >  delete[] oldPos;
54 >  delete[] oldVel;
55 >  delete[] oldJi;
56 > }
57 >
58 > template<typename T> void NPTi<T>::moveA() {
59 >
60 >  //new version of NPTi
61 >  int i, j, k;
62    DirectionalAtom* dAtom;
63 <  double Tb[3];
64 <  double ji[3];
63 >  double Tb[3], ji[3];
64 >  double A[3][3], I[3][3];
65 >  double angle, mass;
66 >  double vel[3], pos[3], frc[3];
67 >
68    double rj[3];
69    double instaTemp, instaPress, instaVol;
70 <  double tt2, tb2;
71 <  double angle;
70 >  double tt2, tb2, scaleFactor;
71 >  double COM[3];
72  
46 –
73    tt2 = tauThermostat * tauThermostat;
74    tb2 = tauBarostat * tauBarostat;
75  
76    instaTemp = tStats->getTemperature();
77    instaPress = tStats->getPressure();
78    instaVol = tStats->getVolume();
53 –  
54 –   // first evolve chi a half step
79    
80 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
81 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
82 <                 (p_convert*NkBT*tb2));
59 <
80 >  tStats->getCOM(COM);
81 >  
82 >  //evolve velocity half step
83    for( i=0; i<nAtoms; i++ ){
61 –    atomIndex = i * 3;
62 –    aMatIndex = i * 9;
63 –    
64 –    // velocity half step
65 –    for( j=atomIndex; j<(atomIndex+3); j++ )
66 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
67 –                       - vel[j]*(chi+eta));
84  
85 <    // position whole step    
85 >    atoms[i]->getVel( vel );
86 >    atoms[i]->getFrc( frc );
87  
88 <    rj[0] = pos[atomIndex];
72 <    rj[1] = pos[atomIndex+1];
73 <    rj[2] = pos[atomIndex+2];
74 <    
75 <    info->wrapVector(rj);
88 >    mass = atoms[i]->getMass();
89  
90 <    pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]);
91 <    pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]);
92 <    pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]);
90 >    for (j=0; j < 3; j++) {
91 >      // velocity half step  (use chi from previous step here):
92 >      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
93 >  
94 >    }
95 >
96 >    atoms[i]->setVel( vel );
97    
98      if( atoms[i]->isDirectional() ){
99  
100        dAtom = (DirectionalAtom *)atoms[i];
101 <          
101 >
102        // get and convert the torque to body frame
103        
104 <      Tb[0] = dAtom->getTx();
88 <      Tb[1] = dAtom->getTy();
89 <      Tb[2] = dAtom->getTz();
90 <      
104 >      dAtom->getTrq( Tb );
105        dAtom->lab2Body( Tb );
106        
107        // get the angular momentum, and propagate a half step
108  
109 <      ji[0] = dAtom->getJx();
110 <      ji[1] = dAtom->getJy();
111 <      ji[2] = dAtom->getJz();
109 >      dAtom->getJ( ji );
110 >
111 >      for (j=0; j < 3; j++)
112 >        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
113        
99 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
100 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
101 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
102 –      
114        // use the angular velocities to propagate the rotation matrix a
115        // full time step
116 <      
116 >
117 >      dAtom->getA(A);
118 >      dAtom->getI(I);
119 >    
120        // rotate about the x-axis      
121 <      angle = dt2 * ji[0] / dAtom->getIxx();
122 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
123 <      
121 >      angle = dt2 * ji[0] / I[0][0];
122 >      this->rotate( 1, 2, angle, ji, A );
123 >
124        // rotate about the y-axis
125 <      angle = dt2 * ji[1] / dAtom->getIyy();
126 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
125 >      angle = dt2 * ji[1] / I[1][1];
126 >      this->rotate( 2, 0, angle, ji, A );
127        
128        // rotate about the z-axis
129 <      angle = dt * ji[2] / dAtom->getIzz();
130 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
129 >      angle = dt * ji[2] / I[2][2];
130 >      this->rotate( 0, 1, angle, ji, A);
131        
132        // rotate about the y-axis
133 <      angle = dt2 * ji[1] / dAtom->getIyy();
134 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
133 >      angle = dt2 * ji[1] / I[1][1];
134 >      this->rotate( 2, 0, angle, ji, A );
135        
136         // rotate about the x-axis
137 <      angle = dt2 * ji[0] / dAtom->getIxx();
138 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
137 >      angle = dt2 * ji[0] / I[0][0];
138 >      this->rotate( 1, 2, angle, ji, A );
139        
140 <      dAtom->setJx( ji[0] );
141 <      dAtom->setJy( ji[1] );
142 <      dAtom->setJz( ji[2] );
140 >      dAtom->setJ( ji );
141 >      dAtom->setA( A  );    
142 >    }    
143 >  }
144 >
145 >  // evolve chi and eta  half step
146 >  
147 >  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
148 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));
149 >
150 >  //calculate the integral of chidt
151 >  integralOfChidt += dt2*chi;
152 >
153 >  //save the old positions
154 >  for(i = 0; i < nAtoms; i++){
155 >    atoms[i]->getPos(pos);
156 >    for(j = 0; j < 3; j++)
157 >      oldPos[i*3 + j] = pos[j];
158 >  }
159 >  
160 >  //the first estimation of r(t+dt) is equal to  r(t)
161 >    
162 >  for(k = 0; k < 4; k ++){
163 >
164 >    for(i =0 ; i < nAtoms; i++){
165 >
166 >      atoms[i]->getVel(vel);
167 >      atoms[i]->getPos(pos);
168 >
169 >      for(j = 0; j < 3; j++)
170 >        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];    
171 >      
172 >      for(j = 0; j < 3; j++)
173 >        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
174 >
175 >      atoms[i]->setPos( pos );
176      }
177      
178 +    if (nConstrained){
179 +      constrainA();
180 +    }
181    }
182 +    
183 +
184    // Scale the box after all the positions have been moved:
185 +  
186 +  scaleFactor = exp(dt*eta);
187  
188 <  cerr << "eta = " << eta
189 <       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";
188 >  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
189 >    sprintf( painCave.errMsg,
190 >             "NPTi error: Attempting a Box scaling of more than 10 percent"
191 >             " check your tauBarostat, as it is probably too small!\n"
192 >             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
193 >             );
194 >    painCave.isFatal = 1;
195 >    simError();
196 >  } else {        
197 >    info->scaleBox(scaleFactor);      
198 >  }  
199  
137 –  info->scaleBox(exp(dt*eta));
138 –
200   }
201  
202 < void NPTi::moveB( void ){
203 <  int i,j,k;
204 <  int atomIndex;
202 > template<typename T> void NPTi<T>::moveB( void ){
203 >  
204 >  //new version of NPTi
205 >  int i, j, k;
206    DirectionalAtom* dAtom;
207 <  double Tb[3];
208 <  double ji[3];
209 <  double instaTemp, instaPress, instaVol;
207 >  double Tb[3], ji[3];
208 >  double vel[3], frc[3];
209 >  double mass;
210 >
211 >  double instTemp, instPress, instVol;
212    double tt2, tb2;
213 +  double oldChi, prevChi;
214 +  double oldEta, preEta;
215    
216    tt2 = tauThermostat * tauThermostat;
217    tb2 = tauBarostat * tauBarostat;
218  
153 –  instaTemp = tStats->getTemperature();
154 –  instaPress = tStats->getPressure();
155 –  instaVol = tStats->getVolume();
219  
220 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
221 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
222 <                 (p_convert*NkBT*tb2));
223 <  
220 >  // Set things up for the iteration:
221 >
222 >  oldChi = chi;
223 >  oldEta = eta;
224 >
225    for( i=0; i<nAtoms; i++ ){
226 <    atomIndex = i * 3;
227 <    
228 <    // velocity half step
229 <    for( j=atomIndex; j<(atomIndex+3); j++ )
230 <    for( j=atomIndex; j<(atomIndex+3); j++ )
231 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168 <                       - vel[j]*(chi+eta));
169 <    
226 >
227 >    atoms[i]->getVel( vel );
228 >
229 >    for (j=0; j < 3; j++)
230 >      oldVel[3*i + j]  = vel[j];
231 >
232      if( atoms[i]->isDirectional() ){
233 <      
233 >
234        dAtom = (DirectionalAtom *)atoms[i];
235 +
236 +      dAtom->getJ( ji );
237 +
238 +      for (j=0; j < 3; j++)
239 +        oldJi[3*i + j] = ji[j];
240 +
241 +    }
242 +  }
243 +
244 +  // do the iteration:
245 +
246 +  instVol = tStats->getVolume();
247 +  
248 +  for (k=0; k < 4; k++) {
249 +    
250 +    instTemp = tStats->getTemperature();
251 +    instPress = tStats->getPressure();
252 +
253 +    // evolve chi another half step using the temperature at t + dt/2
254 +
255 +    prevChi = chi;
256 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
257 +      (tauThermostat*tauThermostat);
258 +
259 +    preEta = eta;
260 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
261 +       (p_convert*NkBT*tb2));
262 +
263 +  
264 +    for( i=0; i<nAtoms; i++ ){
265 +
266 +      atoms[i]->getFrc( frc );
267 +      atoms[i]->getVel(vel);
268        
269 <      // get and convert the torque to body frame
269 >      mass = atoms[i]->getMass();
270        
271 <      Tb[0] = dAtom->getTx();
272 <      Tb[1] = dAtom->getTy();
273 <      Tb[2] = dAtom->getTz();
271 >      // velocity half step
272 >      for (j=0; j < 3; j++)
273 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
274        
275 <      dAtom->lab2Body( Tb );
275 >      atoms[i]->setVel( vel );
276        
277 <      // get the angular momentum, and complete the angular momentum
278 <      // half step
277 >      if( atoms[i]->isDirectional() ){
278 >
279 >        dAtom = (DirectionalAtom *)atoms[i];
280 >  
281 >        // get and convert the torque to body frame      
282 >  
283 >        dAtom->getTrq( Tb );
284 >        dAtom->lab2Body( Tb );      
285 >            
286 >        for (j=0; j < 3; j++)
287 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
288        
289 <      ji[0] = dAtom->getJx();
290 <      ji[1] = dAtom->getJy();
187 <      ji[2] = dAtom->getJz();
188 <      
189 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
190 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
191 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
192 <      
193 <      dAtom->setJx( ji[0] );
194 <      dAtom->setJy( ji[1] );
195 <      dAtom->setJz( ji[2] );
289 >          dAtom->setJ( ji );
290 >      }
291      }
292 +    
293 +    if (nConstrained){
294 +      constrainB();
295 +    }    
296 +    
297 +    if (fabs(prevChi - chi) <=
298 +        chiTolerance && fabs(preEta -eta) <= etaTolerance)
299 +      break;
300    }
301 +
302 +  //calculate integral of chida
303 +  integralOfChidt += dt2*chi;
304 +
305 +
306   }
307  
308 < int NPTi::readyCheck() {
308 > template<typename T> void NPTi<T>::resetIntegrator() {
309 >  chi = 0.0;
310 >  eta = 0.0;
311 > }
312 >
313 > template<typename T> int NPTi<T>::readyCheck() {
314 >
315 >  //check parent's readyCheck() first
316 >  if (T::readyCheck() == -1)
317 >    return -1;
318  
319    // First check to see if we have a target temperature.
320    // Not having one is fatal.
# Line 244 | Line 361 | int NPTi::readyCheck() {
361      return -1;
362    }    
363  
364 +  if (!have_chi_tolerance) {
365 +    sprintf( painCave.errMsg,
366 +             "NPTi warning: setting chi tolerance to 1e-6\n");
367 +    chiTolerance = 1e-6;
368 +    have_chi_tolerance = 1;
369 +    painCave.isFatal = 0;
370 +    simError();
371 +  }
372 +
373 +    if (!have_eta_tolerance) {
374 +    sprintf( painCave.errMsg,
375 +             "NPTi warning: setting eta tolerance to 1e-6\n");
376 +    etaTolerance = 1e-6;
377 +    have_eta_tolerance = 1;
378 +    painCave.isFatal = 0;
379 +    simError();
380 +  }
381    // We need NkBT a lot, so just set it here:
382  
383 <  NkBT = (double)info->ndf * kB * targetTemp;
383 >  NkBT = (double)Nparticles * kB * targetTemp;
384 >  fkBT = (double)info->ndf * kB * targetTemp;
385  
386    return 1;
387   }
388 +
389 + template<typename T> double NPTi<T>::getConservedQuantity(void){
390 +
391 +  double conservedQuantity;
392 +  double LkBT;
393 +  double fkBT;
394 +  double f1kBT;
395 +  double f2kBT;
396 +  double NkBT;
397 +  double Energy;
398 +  double thermostat_kinetic;
399 +  double thermostat_potential;
400 +  double barostat_kinetic;
401 +  double barostat_potential;
402 +  double tb2;
403 +  double eta2;  
404 +
405 +  LkBT = (double)(info->getNDF() + 4) * kB * targetTemp;  // 3N + 1
406 +  fkBT = (double)(info->getNDF()    ) * kB * targetTemp;  // 3N - 3
407 +  f1kBT = (double)(info->getNDF()+ 1) * kB * targetTemp;  // 3N - 3 + 1
408 +  NkBT = (double)(info->getNDF() + 3) * kB * targetTemp;  // 3N
409 +  f2kBT = (double)(info->getNDF()+ 2) * kB * targetTemp;  // 3N - 3 + 1
410 +
411 +  Energy = tStats->getTotalE();
412 +
413 +  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
414 +    (2.0 * eConvert);
415 +
416 +  thermostat_potential = fkBT* integralOfChidt / eConvert;
417 +
418 +
419 +  barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta /
420 +    (2.0 * eConvert);
421 +  
422 +  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
423 +    eConvert;
424 +
425 +  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
426 +    barostat_kinetic + barostat_potential;
427 +  
428 +  cout.width(8);
429 +  cout.precision(8);
430 +
431 +  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
432 +      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
433 +      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
434 +
435 +  return conservedQuantity;
436 + }

Diff Legend

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