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 763 by tim, Mon Sep 15 16:52:02 2003 UTC vs.
Revision 767 by tim, Tue Sep 16 20:02:11 2003 UTC

# Line 56 | Line 56 | template<typename T> void NPTi<T>::moveA() {
56   }
57  
58   template<typename T> void NPTi<T>::moveA() {
59 –
60 –
61 – //   int i, j;
62 – //   DirectionalAtom* dAtom;
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, scaleFactor;
71 –
72 – //   tt2 = tauThermostat * tauThermostat;
73 – //   tb2 = tauBarostat * tauBarostat;
74 –
75 – //   instaTemp = tStats->getTemperature();
76 – //   instaPress = tStats->getPressure();
77 – //   instaVol = tStats->getVolume();
78 –  
79 – //    // first evolve chi a half step
80 –  
81 – //   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
82 – //   eta += dt2 * ( instaVol * (instaPress - targetPressure) /
83 – //               (p_convert*NkBT*tb2));
84 –
85 – //   integralOfChidt += dt2* chi;
86 –
87 – //   for( i=0; i<nAtoms; i++ ){
88 – //     atoms[i]->getVel( vel );
89 – //     atoms[i]->getPos( pos );
90 – //     atoms[i]->getFrc( frc );
91 –
92 – //     mass = atoms[i]->getMass();
93 –
94 – //     for (j=0; j < 3; j++) {
95 – //       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
96 – //       rj[j] = pos[j];
97 – //     }
98 –
99 – //     atoms[i]->setVel( vel );
100 –
101 – //     info->wrapVector(rj);
102 –
103 – //     for (j = 0; j < 3; j++)
104 – //       pos[j] += dt * (vel[j] + eta*rj[j]);
105 –
106 – //     atoms[i]->setPos( pos );
107 –
108 – //     if( atoms[i]->isDirectional() ){
109 –
110 – //       dAtom = (DirectionalAtom *)atoms[i];
111 –          
112 – //       // get and convert the torque to body frame
113 –      
114 – //       dAtom->getTrq( Tb );
115 – //       dAtom->lab2Body( Tb );
116 –      
117 – //       // get the angular momentum, and propagate a half step
118 –
119 – //       dAtom->getJ( ji );
120 –
121 – //       for (j=0; j < 3; j++)
122 – //         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
123 –      
124 – //       // use the angular velocities to propagate the rotation matrix a
125 – //       // full time step
126 –
127 – //       dAtom->getA(A);
128 – //       dAtom->getI(I);
129 –    
130 – //       // rotate about the x-axis      
131 – //       angle = dt2 * ji[0] / I[0][0];
132 – //       this->rotate( 1, 2, angle, ji, A );
133 –
134 – //       // rotate about the y-axis
135 – //       angle = dt2 * ji[1] / I[1][1];
136 – //       this->rotate( 2, 0, angle, ji, A );
137 –      
138 – //       // rotate about the z-axis
139 – //       angle = dt * ji[2] / I[2][2];
140 – //       this->rotate( 0, 1, angle, ji, A);
141 –      
142 – //       // rotate about the y-axis
143 – //       angle = dt2 * ji[1] / I[1][1];
144 – //       this->rotate( 2, 0, angle, ji, A );
145 –      
146 – //        // rotate about the x-axis
147 – //       angle = dt2 * ji[0] / I[0][0];
148 – //       this->rotate( 1, 2, angle, ji, A );
149 –      
150 – //       dAtom->setJ( ji );
151 – //       dAtom->setA( A  );    
152 – //     }                
153 –
154 – //   }
155 –
156 – //   // Scale the box after all the positions have been moved:
157 –  
158 – //   scaleFactor = exp(dt*eta);
159 –
160 – //   if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
161 – //     sprintf( painCave.errMsg,
162 – //              "NPTi error: Attempting a Box scaling of more than 10 percent"
163 – //              " check your tauBarostat, as it is probably too small!\n"
164 – //              " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
165 – //              );
166 – //     painCave.isFatal = 1;
167 – //     simError();
168 – //   } else {        
169 – //     info->scaleBox(exp(dt*eta));      
170 – //   }
171 –  
59  
60    //new version of NPTi
61    int i, j, k;
# Line 280 | Line 167 | template<typename T> void NPTi<T>::moveA() {
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];
284 <
170 >        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];    
171        
286 –      //wrapVector(r(t)) = r(t)-R0
287 –      //info->wrapVector(rj);
288 –      
172        for(j = 0; j < 3; j++)
173          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
174  
# Line 312 | Line 195 | template<typename T> void NPTi<T>::moveA() {
195      info->scaleBox(scaleFactor);      
196    }  
197  
315 –  //advance volume;
316 –  volume = volume * exp(dt*eta);
198   }
199  
200   template<typename T> void NPTi<T>::moveB( void ){
320 –
321 – /*
322 –  int i, j;
323 –  DirectionalAtom* dAtom;
324 –  double Tb[3], ji[3];
325 –  double vel[3], frc[3];
326 –  double mass;
327 –
328 –  double instaTemp, instaPress, instaVol;
329 –  double tt2, tb2;
330 –  
331 –  tt2 = tauThermostat * tauThermostat;
332 –  tb2 = tauBarostat * tauBarostat;
333 –
334 –  instaTemp = tStats->getTemperature();
335 –  instaPress = tStats->getPressure();
336 –  instaVol = tStats->getVolume();
337 –
338 –  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
339 –  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
340 –                 (p_convert*NkBT*tb2));
341 –  integralOfChidt += dt2*chi;
342 –  
343 –  for( i=0; i<nAtoms; i++ ){
344 –
345 –    atoms[i]->getVel( vel );
346 –    atoms[i]->getFrc( frc );
347 –
348 –    mass = atoms[i]->getMass();
349 –
350 –    // velocity half step
351 –    for (j=0; j < 3; j++)
352 –      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
353 –    
354 –    atoms[i]->setVel( vel );
355 –
356 –    if( atoms[i]->isDirectional() ){
357 –
358 –      dAtom = (DirectionalAtom *)atoms[i];
359 –
360 –      // get and convert the torque to body frame      
361 –
362 –      dAtom->getTrq( Tb );
363 –      dAtom->lab2Body( Tb );
364 –
365 –      // get the angular momentum, and propagate a half step
366 –
367 –      dAtom->getJ( ji );
368 –
369 –      for (j=0; j < 3; j++)
370 –        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);    
371 –
372 –      dAtom->setJ( ji );
373 –    }
374 –  }
375 –
376 – */
201    
202    //new version of NPTi
203    int i, j, k;
# Line 566 | Line 390 | template<typename T> double NPTi<T>::getConservedQuant
390    double PV;
391    double extra;
392  
569 –  static double pre_U;
570 –  static double pre_TS;
571 –  static double pre_PV;
572 –  static double pre_extra;
573 –  static int hackCount = 0;
574 –
575 –  double delta_U;
576 –  double delta_TS;
577 –  double delta_PV;
578 –  double delta_extra;
579 –
393    U = tStats->getTotalE();
394  
395    TS = fkBT *
# Line 587 | Line 400 | template<typename T> double NPTi<T>::getConservedQuant
400    tb2 = tauBarostat * tauBarostat;
401    eta2 = eta * eta;
402  
590 –  extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
591 –  /*
592 –  if(hackCount == 0){
593 –    pre_U = U;
594 –    pre_TS =TS;
595 –    pre_PV = PV;
596 –    pre_extra =extra;
597 –    hackCount ++;
598 –  }
403  
404 <  delta_U = U - pre_U;
405 <  delta_TS = TS - pre_TS;
602 <  delta_PV = PV - pre_PV;
603 <  delta_extra = extra - pre_extra;
604 < */
404 >  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;
405 >
406    cout.width(8);
407    cout.precision(8);
408  
# Line 615 | Line 416 | template<typename T> double NPTi<T>::getConservedQuant
416         << extra << "\t"
417         << U+TS+PV+extra << endl;
418  
618 – /*
619 –    pre_U = U;
620 –    pre_TS =TS;
621 –    pre_PV = PV;
622 –    pre_extra =extra;
623 –
624 –
625 –  cout << info->getTime() << "\t"
626 –       << U << "\t"
627 –       << U+TS << "\t"
628 –       << U+TS+PV << "\t"
629 –       << U+TS+PV+extra << endl;
630 – */
419    conservedQuantity = U+TS+PV+extra;
420    return conservedQuantity;
421   }

Diff Legend

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