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 578 by gezelter, Wed Jul 9 02:15:29 2003 UTC vs.
Revision 770 by gezelter, Fri Sep 19 14:55:41 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  
73    tt2 = tauThermostat * tauThermostat;
74    tb2 = tauBarostat * tauBarostat;
# Line 49 | Line 76 | void NPTi::moveA() {
76    instaTemp = tStats->getTemperature();
77    instaPress = tStats->getPressure();
78    instaVol = tStats->getVolume();
52 –  
53 –  // first evolve chi a half step
79    
80 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
81 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
82 <
80 >  tStats->getCOM(COM);
81 >  
82 >  //evolve velocity half step
83    for( i=0; i<nAtoms; i++ ){
59 –    atomIndex = i * 3;
60 –    aMatIndex = i * 9;
61 –    
62 –    // velocity half step
63 –    for( j=atomIndex; j<(atomIndex+3); j++ )
64 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
65 –                       - 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];
70 <    rj[1] = pos[atomIndex+1];
71 <    rj[2] = pos[atomIndex+2];
72 <    
73 <    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();
86 <      Tb[1] = dAtom->getTy();
87 <      Tb[2] = dAtom->getTz();
88 <      
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        
97 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
98 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
99 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
100 –      
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 <  info->scaleBox(exp(dt*eta));
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  
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  
219 <  instaTemp = tStats->getTemperature();
149 <  instaPress = tStats->getPressure();
150 <  instaVol = tStats->getVolume();
219 >  // Set things up for the iteration:
220  
221 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
222 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
223 <  
221 >  oldChi = chi;
222 >  oldEta = eta;
223 >
224    for( i=0; i<nAtoms; i++ ){
225 <    atomIndex = i * 3;
226 <    
227 <    // velocity half step
228 <    for( j=atomIndex; j<(atomIndex+3); j++ )
229 <    for( j=atomIndex; j<(atomIndex+3); j++ )
230 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
162 <                       - vel[j]*(chi+eta));
163 <    
225 >
226 >    atoms[i]->getVel( vel );
227 >
228 >    for (j=0; j < 3; j++)
229 >      oldVel[3*i + j]  = vel[j];
230 >
231      if( atoms[i]->isDirectional() ){
232 <      
232 >
233        dAtom = (DirectionalAtom *)atoms[i];
234 +
235 +      dAtom->getJ( ji );
236 +
237 +      for (j=0; j < 3; j++)
238 +        oldJi[3*i + j] = ji[j];
239 +
240 +    }
241 +  }
242 +
243 +  // do the iteration:
244 +
245 +  instVol = tStats->getVolume();
246 +  
247 +  for (k=0; k < 4; k++) {
248 +    
249 +    instTemp = tStats->getTemperature();
250 +    instPress = tStats->getPressure();
251 +
252 +    // evolve chi another half step using the temperature at t + dt/2
253 +
254 +    prevChi = chi;
255 +    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
256 +      (tauThermostat*tauThermostat);
257 +
258 +    preEta = eta;
259 +    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
260 +       (p_convert*NkBT*tb2));
261 +
262 +  
263 +    for( i=0; i<nAtoms; i++ ){
264 +
265 +      atoms[i]->getFrc( frc );
266 +      atoms[i]->getVel(vel);
267        
268 <      // get and convert the torque to body frame
268 >      mass = atoms[i]->getMass();
269        
270 <      Tb[0] = dAtom->getTx();
271 <      Tb[1] = dAtom->getTy();
272 <      Tb[2] = dAtom->getTz();
270 >      // velocity half step
271 >      for (j=0; j < 3; j++)
272 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
273        
274 <      dAtom->lab2Body( Tb );
274 >      atoms[i]->setVel( vel );
275        
276 <      // get the angular momentum, and complete the angular momentum
277 <      // half step
276 >      if( atoms[i]->isDirectional() ){
277 >
278 >        dAtom = (DirectionalAtom *)atoms[i];
279 >  
280 >        // get and convert the torque to body frame      
281 >  
282 >        dAtom->getTrq( Tb );
283 >        dAtom->lab2Body( Tb );      
284 >            
285 >        for (j=0; j < 3; j++)
286 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
287        
288 <      ji[0] = dAtom->getJx();
289 <      ji[1] = dAtom->getJy();
181 <      ji[2] = dAtom->getJz();
182 <      
183 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
184 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
185 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
186 <      
187 <      dAtom->setJx( ji[0] );
188 <      dAtom->setJy( ji[1] );
189 <      dAtom->setJz( ji[2] );
288 >          dAtom->setJ( ji );
289 >      }
290      }
291 +    
292 +    if (nConstrained){
293 +      constrainB();
294 +    }    
295 +    
296 +    if (fabs(prevChi - chi) <=
297 +        chiTolerance && fabs(preEta -eta) <= etaTolerance)
298 +      break;
299    }
300 +
301 +  //calculate integral of chida
302 +  integralOfChidt += dt2*chi;
303 +
304 +
305   }
306  
307 < int NPTi::readyCheck() {
307 > template<typename T> void NPTi<T>::resetIntegrator() {
308 >  chi = 0.0;
309 >  eta = 0.0;
310 > }
311 >
312 > template<typename T> int NPTi<T>::readyCheck() {
313 >
314 >  //check parent's readyCheck() first
315 >  if (T::readyCheck() == -1)
316 >    return -1;
317  
318    // First check to see if we have a target temperature.
319    // Not having one is fatal.
# Line 238 | Line 360 | int NPTi::readyCheck() {
360      return -1;
361    }    
362  
363 <  // We need NkBT a lot, so just set it here:
363 >  if (!have_chi_tolerance) {
364 >    sprintf( painCave.errMsg,
365 >             "NPTi warning: setting chi tolerance to 1e-6\n");
366 >    chiTolerance = 1e-6;
367 >    have_chi_tolerance = 1;
368 >    painCave.isFatal = 0;
369 >    simError();
370 >  }
371  
372 <  NkBT = (double)info->ndf * kB * targetTemp;
372 >  if (!have_eta_tolerance) {
373 >    sprintf( painCave.errMsg,
374 >             "NPTi warning: setting eta tolerance to 1e-6\n");
375 >    etaTolerance = 1e-6;
376 >    have_eta_tolerance = 1;
377 >    painCave.isFatal = 0;
378 >    simError();
379 >  }
380 >  
381 >  
382 >  // We need NkBT a lot, so just set it here: This is the RAW number
383 >  // of particles, so no subtraction or addition of constraints or
384 >  // orientational degrees of freedom:
385 >  
386 >  NkBT = (double)Nparticles * kB * targetTemp;
387 >  
388 >  // fkBT is used because the thermostat operates on more degrees of freedom
389 >  // than the barostat (when there are particles with orientational degrees
390 >  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
391 >  
392 >  fkBT = (double)info->ndf * kB * targetTemp;
393  
394    return 1;
395   }
396 +
397 + template<typename T> double NPTi<T>::getConservedQuantity(void){
398 +
399 +  double conservedQuantity;
400 +  double Three_NkBT;
401 +  double Energy;
402 +  double thermostat_kinetic;
403 +  double thermostat_potential;
404 +  double barostat_kinetic;
405 +  double barostat_potential;
406 +  double tb2;
407 +  double eta2;
408 +
409 +  Energy = tStats->getTotalE();
410 +
411 +  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
412 +    (2.0 * eConvert);
413 +
414 +  thermostat_potential = fkBT* integralOfChidt / eConvert;
415 +
416 +
417 +  barostat_kinetic = 3.0 * NkBT * tauBarostat * tauBarostat * eta * eta /
418 +    (2.0 * eConvert);
419 +  
420 +  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
421 +    eConvert;
422 +
423 +  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
424 +    barostat_kinetic + barostat_potential;
425 +  
426 +  cout.width(8);
427 +  cout.precision(8);
428 +
429 +  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
430 +      "\t" << thermostat_potential << "\t" << barostat_kinetic <<
431 +      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
432 +
433 +  return conservedQuantity;
434 + }

Diff Legend

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