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

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 574 by gezelter, Tue Jul 8 20:56:10 2003 UTC vs.
Revision 768 by mmeineke, Wed Sep 17 14:22:15 2003 UTC

# Line 1 | Line 1
1 + #include <cmath>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# Line 8 | 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 19 | Line 24
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 48 | Line 76 | void NPTi::moveA() {
76    instaTemp = tStats->getTemperature();
77    instaPress = tStats->getPressure();
78    instaVol = tStats->getVolume();
51 –  
52 –  // 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++ ){
58 –    atomIndex = i * 3;
59 –    aMatIndex = i * 9;
60 –    
61 –    // velocity half step
62 –    for( j=atomIndex; j<(atomIndex+3); j++ )
63 –      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
64 –                       - vel[j]*(chi+eta));
84  
85 <    // position whole step    
85 >    atoms[i]->getVel( vel );
86 >    atoms[i]->getFrc( frc );
87  
88 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
69 <      rj[0] = pos[j];
70 <      rj[1] = pos[j+1];
71 <      rj[2] = pos[j+2];
88 >    mass = atoms[i]->getMass();
89  
90 <      info->wrapVector(rj);
91 <
92 <      pos[j] += dt * (vel[j] + eta*rj[0]);
93 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
77 <      pos[j+2] += dt * (vel[j+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 <    // Scale the box after all the positions have been moved:
81 <
82 <    info->scaleBox(dt*eta);
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();
91 <      Tb[1] = dAtom->getTy();
92 <      Tb[2] = dAtom->getTz();
93 <      
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        
102 –      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
103 –      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
104 –      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
105 –      
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 +  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  
149 –  instaTemp = tStats->getTemperature();
150 –  instaPress = tStats->getPressure();
151 –  instaVol = tStats->getVolume();
219  
220 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
221 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
222 <  
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
163 <                       - vel[j]*(chi+eta));
164 <    
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();
182 <      ji[2] = dAtom->getJz();
183 <      
184 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
185 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
186 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
187 <      
188 <      dAtom->setJx( ji[0] );
189 <      dAtom->setJy( ji[1] );
190 <      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 239 | 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 tb2;
393 +  double eta2;  
394 +  double E_NPT;
395 +  double U;
396 +  double TS;
397 +  double PV;
398 +  double extra;
399 +
400 +  U = tStats->getTotalE();
401 +
402 +  TS = fkBT *
403 +    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
404 +
405 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
406 +
407 +  tb2 = tauBarostat * tauBarostat;
408 +  eta2 = eta * eta;
409 +
410 +
411 +  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;
412 +
413 +  cout.width(8);
414 +  cout.precision(8);
415 +
416 +  
417 + //   cout << info->getTime() << "\t"
418 + //        << chi << "\t"
419 + //        << eta << "\t"
420 + //        << U << "\t"
421 + //        << TS << "\t"
422 + //        << PV << "\t"
423 + //        << extra << "\t"
424 + //        << U+TS+PV+extra << endl;
425 +
426 +  conservedQuantity = U+TS+PV+extra;
427 +  return conservedQuantity;
428 + }

Diff Legend

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