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 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 25 | Line 25 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
25   {
26    int i, j;
27    chi = 0.0;
28 +  integralOfChidt = 0.0;
29  
30    for(i = 0; i < 3; i++)
31      for (j = 0; j < 3; j++)
# Line 34 | Line 35 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
35    have_tau_barostat = 0;
36    have_target_temp = 0;
37    have_target_pressure = 0;
38 +
39 +  have_chi_tolerance = 0;
40 +  have_eta_tolerance = 0;
41 +  have_pos_iter_tolerance = 0;
42 +
43 +  oldPos = new double[3*nAtoms];
44 +  oldVel = new double[3*nAtoms];
45 +  oldJi = new double[3*nAtoms];
46 + #ifdef IS_MPI
47 +  Nparticles = mpiSim->getTotAtoms();
48 + #else
49 +  Nparticles = theInfo->n_atoms;
50 + #endif
51   }
52  
53 + template<typename T> NPTf<T>::~NPTf() {
54 +  delete[] oldPos;
55 +  delete[] oldVel;
56 +  delete[] oldJi;
57 + }
58 +
59   template<typename T> void NPTf<T>::moveA() {
60    
61    int i, j, k;
# Line 52 | Line 72 | template<typename T> void NPTf<T>::moveA() {
72    double eta2ij;
73    double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
74    double bigScale, smallScale, offDiagMax;
75 +  double COM[3];
76  
77    tt2 = tauThermostat * tauThermostat;
78    tb2 = tauBarostat * tauBarostat;
# Line 59 | Line 80 | template<typename T> void NPTf<T>::moveA() {
80    instaTemp = tStats->getTemperature();
81    tStats->getPressureTensor(press);
82    instaVol = tStats->getVolume();
62 –  
63 –  // first evolve chi a half step
83    
84 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
84 >  tStats->getCOM(COM);
85  
86 +  //calculate scale factor of veloity
87    for (i = 0; i < 3; i++ ) {
88      for (j = 0; j < 3; j++ ) {
89 +      vScale[i][j] = eta[i][j];
90 +      
91        if (i == j) {
92 <        
93 <        eta[i][j] += dt2 * instaVol *
72 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
73 <        
74 <        vScale[i][j] = eta[i][j] + chi;
75 <        
76 <      } else {
77 <        
78 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
79 <
80 <        vScale[i][j] = eta[i][j];
81 <        
82 <      }
92 >        vScale[i][j] += chi;          
93 >      }              
94      }
95    }
96 <
96 >  
97 >  //evolve velocity half step
98    for( i=0; i<nAtoms; i++ ){
99  
100      atoms[i]->getVel( vel );
89 –    atoms[i]->getPos( pos );
101      atoms[i]->getFrc( frc );
102  
103      mass = atoms[i]->getMass();
104      
94 –    // velocity half step
95 –        
105      info->matVecMul3( vScale, vel, sc );
106 <    
107 <    for (j = 0; j < 3; j++) {
106 >
107 >    for (j=0; j < 3; j++) {
108 >      // velocity half step  (use chi from previous step here):
109        vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
110 <      rj[j] = pos[j];
110 >  
111      }
112  
113      atoms[i]->setVel( vel );
104 –
105 –    // position whole step    
106 –
107 –    info->wrapVector(rj);
108 –
109 –    info->matVecMul3( eta, rj, sc );
110 –
111 –    for (j = 0; j < 3; j++ )
112 –      pos[j] += dt * (vel[j] + sc[j]);
113 –
114 –    atoms[i]->setPos( pos );
114    
115      if( atoms[i]->isDirectional() ){
116  
117        dAtom = (DirectionalAtom *)atoms[i];
118 <          
118 >
119        // get and convert the torque to body frame
120        
121        dAtom->getTrq( Tb );
# Line 157 | Line 156 | template<typename T> void NPTf<T>::moveA() {
156        
157        dAtom->setJ( ji );
158        dAtom->setA( A  );    
159 <    }                    
159 >    }    
160    }
161 +
162 +  // advance chi half step
163 +  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
164 +
165 +  //calculate the integral of chidt
166 +  integralOfChidt += dt2*chi;
167 +
168 +  //advance eta half step
169 +  for(i = 0; i < 3; i ++)
170 +    for(j = 0; j < 3; j++){
171 +      if( i == j)
172 +        eta[i][j] += dt2 *  instaVol *
173 +          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
174 +      else
175 +        eta[i][j] += dt2 * instaVol * press[i][j] / ( NkBT*tb2);
176 +    }
177 +    
178 +  //save the old positions
179 +  for(i = 0; i < nAtoms; i++){
180 +    atoms[i]->getPos(pos);
181 +    for(j = 0; j < 3; j++)
182 +      oldPos[i*3 + j] = pos[j];
183 +  }
184    
185 +  //the first estimation of r(t+dt) is equal to  r(t)
186 +    
187 +  for(k = 0; k < 4; k ++){
188 +
189 +    for(i =0 ; i < nAtoms; i++){
190 +
191 +      atoms[i]->getVel(vel);
192 +      atoms[i]->getPos(pos);
193 +
194 +      for(j = 0; j < 3; j++)
195 +        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
196 +      
197 +      info->matVecMul3( eta, rj, sc );
198 +      
199 +      for(j = 0; j < 3; j++)
200 +        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
201 +
202 +      atoms[i]->setPos( pos );
203 +
204 +    }
205 +
206 +  }  
207 +
208 +
209    // Scale the box after all the positions have been moved:
210    
211    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
# Line 189 | Line 235 | template<typename T> void NPTf<T>::moveA() {
235        if (i != j)
236          if (fabs(scaleMat[i][j]) > offDiagMax)
237            offDiagMax = fabs(scaleMat[i][j]);
192 –      
238      }
239  
240      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 230 | Line 275 | template<typename T> void NPTf<T>::moveB( void ){
275  
276   template<typename T> void NPTf<T>::moveB( void ){
277  
278 <  int i, j;
278 >  int i, j, k;
279    DirectionalAtom* dAtom;
280    double Tb[3], ji[3];
281    double vel[3], frc[3];
# Line 240 | Line 285 | template<typename T> void NPTf<T>::moveB( void ){
285    double tt2, tb2;
286    double sc[3];
287    double press[3][3], vScale[3][3];
288 +  double oldChi, prevChi;
289 +  double oldEta[3][3], preEta[3][3], diffEta;
290    
291    tt2 = tauThermostat * tauThermostat;
292    tb2 = tauBarostat * tauBarostat;
293  
247 –  instaTemp = tStats->getTemperature();
248 –  tStats->getPressureTensor(press);
249 –  instaVol = tStats->getVolume();
250 –  
251 –  // first evolve chi a half step
252 –  
253 –  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
254 –  
255 –  for (i = 0; i < 3; i++ ) {
256 –    for (j = 0; j < 3; j++ ) {
257 –      if (i == j) {
294  
295 <        eta[i][j] += dt2 * instaVol *
260 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
295 >  // Set things up for the iteration:
296  
297 <        vScale[i][j] = eta[i][j] + chi;
298 <        
299 <      } else {
300 <        
301 <        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
297 >  oldChi = chi;
298 >  
299 >  for(i = 0; i < 3; i++)
300 >    for(j = 0; j < 3; j++)
301 >      oldEta[i][j] = eta[i][j];
302  
268 –        vScale[i][j] = eta[i][j];
269 –        
270 –      }
271 –    }
272 –  }
273 –
303    for( i=0; i<nAtoms; i++ ){
304  
305      atoms[i]->getVel( vel );
277 –    atoms[i]->getFrc( frc );
306  
307 <    mass = atoms[i]->getMass();
308 <    
281 <    // velocity half step
282 <        
283 <    info->matVecMul3( vScale, vel, sc );
284 <    
285 <    for (j = 0; j < 3; j++) {
286 <      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
287 <    }
307 >    for (j=0; j < 3; j++)
308 >      oldVel[3*i + j]  = vel[j];
309  
289 –    atoms[i]->setVel( vel );
290 –    
310      if( atoms[i]->isDirectional() ){
311  
312        dAtom = (DirectionalAtom *)atoms[i];
313 <          
314 <      // get and convert the torque to body frame
313 >
314 >      dAtom->getJ( ji );
315 >
316 >      for (j=0; j < 3; j++)
317 >        oldJi[3*i + j] = ji[j];
318 >
319 >    }
320 >  }
321 >
322 >  // do the iteration:
323 >
324 >  instaVol = tStats->getVolume();
325 >  
326 >  for (k=0; k < 4; k++) {
327 >    
328 >    instaTemp = tStats->getTemperature();
329 >    tStats->getPressureTensor(press);
330 >
331 >    // evolve chi another half step using the temperature at t + dt/2
332 >
333 >    prevChi = chi;
334 >    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
335 >    
336 >    for(i = 0; i < 3; i++)
337 >      for(j = 0; j < 3; j++)
338 >        preEta[i][j] = eta[i][j];
339 >
340 >    //advance eta half step and calculate scale factor for velocity
341 >    for(i = 0; i < 3; i ++)
342 >      for(j = 0; j < 3; j++){
343 >        if( i == j){
344 >          eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
345 >            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
346 >          vScale[i][j] = eta[i][j] + chi;
347 >        }
348 >        else
349 >        {
350 >          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
351 >          vScale[i][j] = eta[i][j];
352 >        }
353 >    }      
354 >
355 >    //advance velocity half step
356 >    for( i=0; i<nAtoms; i++ ){
357 >
358 >      atoms[i]->getFrc( frc );
359 >      atoms[i]->getVel(vel);
360        
361 <      dAtom->getTrq( Tb );
298 <      dAtom->lab2Body( Tb );
361 >      mass = atoms[i]->getMass();
362        
363 <      // get the angular momentum, and propagate a half step
363 >      info->matVecMul3( vScale, vel, sc );
364 >
365 >      for (j=0; j < 3; j++) {
366 >        // velocity half step  (use chi from previous step here):
367 >        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
368 >      }
369        
370 <      dAtom->getJ( ji );
370 >      atoms[i]->setVel( vel );
371        
372 <      for (j=0; j < 3; j++)
373 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
372 >      if( atoms[i]->isDirectional() ){
373 >
374 >        dAtom = (DirectionalAtom *)atoms[i];
375 >  
376 >        // get and convert the torque to body frame      
377 >  
378 >        dAtom->getTrq( Tb );
379 >        dAtom->lab2Body( Tb );      
380 >            
381 >        for (j=0; j < 3; j++)
382 >          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
383        
384 <      dAtom->setJ( ji );
384 >          dAtom->setJ( ji );
385 >      }
386 >    }
387  
388 <    }                    
388 >    
389 >    diffEta = 0;
390 >    for(i = 0; i < 3; i++)
391 >      diffEta += pow(preEta[i][i] - eta[i][i], 2);    
392 >    
393 >    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
394 >      break;
395    }
396 +
397 +  //calculate integral of chida
398 +  integralOfChidt += dt2*chi;
399 +
400 +  
401   }
402  
403 + template<typename T> void NPTf<T>::resetIntegrator() {
404 +  int i,j;
405 +  
406 +  chi = 0.0;
407 +
408 +  for(i = 0; i < 3; i++)
409 +    for (j = 0; j < 3; j++)
410 +      eta[i][j] = 0.0;
411 +
412 + }
413 +
414   template<typename T> int NPTf<T>::readyCheck() {
415  
416    //check parent's readyCheck() first
# Line 363 | Line 464 | template<typename T> int NPTf<T>::readyCheck() {
464  
465    // We need NkBT a lot, so just set it here:
466  
467 <  NkBT = (double)info->ndf * kB * targetTemp;
467 >  NkBT = (double)Nparticles * kB * targetTemp;
468 >  fkBT = (double)info->ndf * kB * targetTemp;
469  
470    return 1;
471   }
472 +
473 + template<typename T> double NPTf<T>::getConservedQuantity(void){
474 +
475 +  double conservedQuantity;
476 +  double tb2;
477 +  double trEta;  
478 +  double U;
479 +  double thermo;
480 +  double integral;
481 +  double baro;
482 +  double PV;
483 +
484 +  U = tStats->getTotalE();
485 +  thermo = (fkBT * tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
486 +
487 +  tb2 = tauBarostat * tauBarostat;
488 +  trEta = info->matTrace3(eta);
489 +  baro = ((double)info->ndfTrans * kB * targetTemp * tb2 * trEta * trEta / 2.0) / eConvert;
490 +
491 +  integral = ((double)(info->ndf + 1) * kB * targetTemp * integralOfChidt) /eConvert;
492 +
493 +  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
494 +
495 +
496 +  cout.width(8);
497 +  cout.precision(8);
498 +  
499 +  cout << info->getTime() << "\t"
500 +       << chi << "\t"
501 +       << trEta << "\t"
502 +       << U << "\t"
503 +       << thermo << "\t"
504 +       << baro << "\t"
505 +       << integral << "\t"
506 +       << PV << "\t"
507 +       << U+thermo+integral+PV+baro << endl;
508 +
509 +  conservedQuantity = U+thermo+integral+PV+baro;
510 +  return conservedQuantity;
511 +  
512 + }

Diff Legend

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