ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/OOPSE/libmdtools/NPT.cpp
(Generate patch)

Comparing:
trunk/OOPSE/libmdtools/NPT.cpp (file contents), Revision 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
branches/new-templateless/OOPSE/libmdtools/NPT.cpp (file contents), Revision 851 by mmeineke, Wed Nov 5 19:18:17 2003 UTC

# Line 1 | Line 1
1 + #include <stdlib.h>
2   #include <math.h>
3 +
4   #include "Atom.hpp"
5   #include "SRI.hpp"
6   #include "AbstractClasses.hpp"
# Line 7 | Line 9
9   #include "Thermo.hpp"
10   #include "ReadWrite.hpp"
11   #include "Integrator.hpp"
12 < #include "simError.h"
12 > #include "simError.h"
13  
14   #ifdef IS_MPI
15   #include "mpiSimulation.hpp"
# Line 18 | Line 20
20   // modification of the Hoover algorithm:
21   //
22   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
23 < //       Molec. Phys., 78, 533.
23 > //       Molec. Phys., 78, 533.
24   //
25   //           and
26 < //
26 > //
27   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28  
29 < template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
30 <  T( theInfo, the_ff )
29 > NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff):
30 >  Integrator( theInfo, the_ff )
31   {
32 +  GenericData* data;
33 +  
34    chi = 0.0;
35    integralOfChidt = 0.0;
36    have_tau_thermostat = 0;
# Line 37 | Line 41 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
41    have_eta_tolerance = 0;
42    have_pos_iter_tolerance = 0;
43  
44 +  // retrieve chi and integralOfChidt from simInfo
45 +  data = info->getProperty(CHIVALUE_ID);
46 +  if(data != NULL ){
47 +    chi = data->getDval();
48 +  }
49 +
50 +  data = info->getProperty(INTEGRALOFCHIDT_ID);
51 +  if(data != NULL ){
52 +    integralOfChidt = data->getDval();
53 +  }
54 +
55    oldPos = new double[3*nAtoms];
56    oldVel = new double[3*nAtoms];
57    oldJi = new double[3*nAtoms];
# Line 48 | Line 63 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
63  
64   }
65  
66 < template<typename T> NPT<T>::~NPT() {
66 > NPT::~NPT() {
67    delete[] oldPos;
68    delete[] oldVel;
69    delete[] oldJi;
70   }
71  
72 < template<typename T> void NPT<T>::moveA() {
73 <
72 > void NPT::moveA() {
73 >
74    //new version of NPT
75    int i, j, k;
76    DirectionalAtom* dAtom;
# Line 69 | Line 84 | template<typename T> void NPT<T>::moveA() {
84    tStats->getPressureTensor( press );
85    instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
86    instaVol = tStats->getVolume();
87 <  
87 >
88    tStats->getCOM(COM);
89 <  
89 >
90    //evolve velocity half step
91    for( i=0; i<nAtoms; i++ ){
92  
# Line 83 | Line 98 | template<typename T> void NPT<T>::moveA() {
98      getVelScaleA( sc, vel );
99  
100      for (j=0; j < 3; j++) {
101 <      
101 >
102        // velocity half step  (use chi from previous step here):
103        vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
104 <  
104 >
105      }
106  
107      atoms[i]->setVel( vel );
108 <  
108 >
109      if( atoms[i]->isDirectional() ){
110  
111        dAtom = (DirectionalAtom *)atoms[i];
112  
113        // get and convert the torque to body frame
114 <      
114 >
115        dAtom->getTrq( Tb );
116        dAtom->lab2Body( Tb );
117 <      
117 >
118        // get the angular momentum, and propagate a half step
119  
120        dAtom->getJ( ji );
121  
122 <      for (j=0; j < 3; j++)
122 >      for (j=0; j < 3; j++)
123          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
124 <      
124 >
125        this->rotationPropagation( dAtom, ji );
126 <      
126 >
127        dAtom->setJ( ji );
128 <    }    
128 >    }
129    }
130  
131    // evolve chi and eta  half step
132 <  
132 >
133    evolveChiA();
134    evolveEtaA();
135  
# Line 127 | Line 142 | template<typename T> void NPT<T>::moveA() {
142      for(j = 0; j < 3; j++)
143        oldPos[i*3 + j] = pos[j];
144    }
145 <  
146 <  //the first estimation of r(t+dt) is equal to  r(t)
147 <    
145 >
146 >  //the first estimation of r(t+dt) is equal to  r(t)
147 >
148    for(k = 0; k < 5; k ++){
149  
150      for(i =0 ; i < nAtoms; i++){
# Line 138 | Line 153 | template<typename T> void NPT<T>::moveA() {
153        atoms[i]->getPos(pos);
154  
155        this->getPosScale( pos, COM, i, sc );
156 <  
156 >
157        for(j = 0; j < 3; j++)
158          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
159  
160        atoms[i]->setPos( pos );
161      }
162 <    
162 >
163      if (nConstrained){
164        constrainA();
165      }
166    }
152 –    
167  
168 +
169    // Scale the box after all the positions have been moved:
170 <  
170 >
171    this->scaleSimBox();
172   }
173  
174 < template<typename T> void NPT<T>::moveB( void ){
175 <  
174 > void NPT::moveB( void ){
175 >
176    //new version of NPT
177    int i, j, k;
178    DirectionalAtom* dAtom;
179    double Tb[3], ji[3], sc[3];
180 <  double vel[3], frc[3];
180 >  double vel[3], frc[3], qVel[3];
181    double mass;
182 <  
182 >
183    // Set things up for the iteration:
184  
185    for( i=0; i<nAtoms; i++ ){
# Line 189 | Line 204 | template<typename T> void NPT<T>::moveB( void ){
204    // do the iteration:
205  
206    instaVol = tStats->getVolume();
207 <  
207 >
208    for (k=0; k < 4; k++) {
209 +
210 +    atoms[0]->getVel(vel);
211 +
212 + //    std::cerr << "vel atom0 = "
213 + //            << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";
214      
215      instaTemp = tStats->getTemperature();
216 +    
217 +    std::cerr << "instaTemp MoveB = " << instaTemp << "\n";
218 +    
219      instaPress = tStats->getPressure();
220  
221      // evolve chi another half step using the temperature at t + dt/2
222  
223      this->evolveChiB();
224      this->evolveEtaB();
225 <    
225 >
226      for( i=0; i<nAtoms; i++ ){
227  
228        atoms[i]->getFrc( frc );
229        atoms[i]->getVel(vel);
230 <      
230 >
231        mass = atoms[i]->getMass();
232 <      
232 >
233        getVelScaleB( sc, i );
234  
235 +      for(j=0;j<3;j++)
236 +        qVel[j] = vel[j];
237 +
238        // velocity half step
239 <      for (j=0; j < 3; j++)
239 >      for (j=0; j < 3; j++)
240          vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
241 +
242 + //       std::cerr << " vel of " << i << " => "
243 + //              << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";
244 +
245 +      double delVel = 0.0;
246 +      for(j=0;j<3;j++)
247 +        delVel += pow( (qVel[j] - vel[j]), 2 );
248        
249 <      atoms[i]->setVel( vel );
249 >      delVel = sqrt(delVel);
250        
251 +      std::cerr << "delVel[" << i << "] = " << delVel << "\n";
252 +
253 +      atoms[i]->setVel( vel );
254 +
255        if( atoms[i]->isDirectional() ){
256  
257          dAtom = (DirectionalAtom *)atoms[i];
258 <  
259 <        // get and convert the torque to body frame      
260 <  
261 <        dAtom->getTrq( Tb );
262 <        dAtom->lab2Body( Tb );      
263 <            
264 <        for (j=0; j < 3; j++)
258 >
259 >        // get and convert the torque to body frame
260 >
261 >        dAtom->getTrq( Tb );
262 >        dAtom->lab2Body( Tb );
263 >
264 >        for (j=0; j < 3; j++)
265            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
266 <      
266 >
267            dAtom->setJ( ji );
268        }
269      }
270 <    
270 >
271      if (nConstrained){
272        constrainB();
273 <    }    
274 <    
273 >    }
274 >
275      if ( this->chiConverged() && this->etaConverged() ) break;
276    }
277  
# Line 244 | Line 281 | template<typename T> void NPT<T>::moveB( void ){
281  
282   }
283  
284 < template<typename T> void NPT<T>::resetIntegrator() {
284 > void NPT::resetIntegrator() {
285    chi = 0.0;
286 <  T::resetIntegrator();
286 >  Integrator::resetIntegrator();
287   }
288  
289 < template<typename T> void NPT<T>::evolveChiA() {
289 > void NPT::evolveChiA() {
290 >
291 >  std::cerr << "\n"
292 >            << "evolveChiA before:\n"
293 >            << "chi = " << chi << "\n"
294 >            << "oldChi = " << oldChi << "\n"
295 >            << "\n"
296 >            << "instaTemp = " << instaTemp << "\n"
297 >            << "targetTemp = " << targetTemp << "\n"
298 >            << "dt2 = " << dt2 << "\n"
299 >            << "tt2 = " << tt2 << "\n";
300 >
301    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
302    oldChi = chi;
303 +
304 +  std::cerr << "\n"
305 +            << "evolveChiA after:\n"
306 +            << "chi = " << chi << "\n"
307 +            << "oldChi = " << oldChi << "\n"
308 +            << "\n"
309 +            << "instaTemp = " << instaTemp << "\n"
310 +            << "targetTemp = " << targetTemp << "\n"
311 +            << "dt2 = " << dt2 << "\n"
312 +            << "tt2 = " << tt2 << "\n";
313 +
314   }
315  
316 < template<typename T> void NPT<T>::evolveChiB() {
317 <  
316 > void NPT::evolveChiB() {
317 >
318 >  std::cerr << "\n"
319 >            << "evolveChiB before:\n"
320 >            << "chi = " << chi << "\n"
321 >            << "oldChi = " << oldChi << "\n"
322 >            << "\n"
323 >            << "instaTemp = " << instaTemp << "\n"
324 >            << "targetTemp = " << targetTemp << "\n"
325 >            << "dt2 = " << dt2 << "\n"
326 >            << "tt2 = " << tt2 << "\n";
327 >
328    prevChi = chi;
329    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
330 +
331 +  std::cerr << "\n"
332 +            << "evolveChiB after:\n"
333 +            << "chi = " << chi << "\n"
334 +            << "oldChi = " << oldChi << "\n"
335 +            << "\n"
336 +            << "instaTemp = " << instaTemp << "\n"
337 +            << "targetTemp = " << targetTemp << "\n"
338 +            << "dt2 = " << dt2 << "\n"
339 +            << "tt2 = " << tt2 << "\n";
340   }
341  
342 < template<typename T> bool NPT<T>::chiConverged() {
343 <  
344 <  return ( fabs( prevChi - chi ) <= chiTolerance );  
342 > bool NPT::chiConverged() {
343 >
344 >  return ( fabs( prevChi - chi ) <= chiTolerance );
345   }
346  
347 < template<typename T> int NPT<T>::readyCheck() {
347 > int NPT::readyCheck() {
348  
349    //check parent's readyCheck() first
350 <  if (T::readyCheck() == -1)
350 >  if (Integrator::readyCheck() == -1)
351      return -1;
352 <
353 <  // First check to see if we have a target temperature.
354 <  // Not having one is fatal.
355 <  
352 >
353 >  // First check to see if we have a target temperature.
354 >  // Not having one is fatal.
355 >
356    if (!have_target_temp) {
357      sprintf( painCave.errMsg,
358               "NPT error: You can't use the NPT integrator\n"
# Line 293 | Line 372 | template<typename T> int NPT<T>::readyCheck() {
372      simError();
373      return -1;
374    }
375 <  
375 >
376    // We must set tauThermostat.
377 <  
377 >
378    if (!have_tau_thermostat) {
379      sprintf( painCave.errMsg,
380               "NPT error: If you use the NPT\n"
# Line 303 | Line 382 | template<typename T> int NPT<T>::readyCheck() {
382      painCave.isFatal = 1;
383      simError();
384      return -1;
385 <  }    
385 >  }
386  
387    // We must set tauBarostat.
388 <  
388 >
389    if (!have_tau_barostat) {
390      sprintf( painCave.errMsg,
391               "NPT error: If you use the NPT\n"
# Line 314 | Line 393 | template<typename T> int NPT<T>::readyCheck() {
393      painCave.isFatal = 1;
394      simError();
395      return -1;
396 <  }    
396 >  }
397  
398    if (!have_chi_tolerance) {
399      sprintf( painCave.errMsg,
# Line 323 | Line 402 | template<typename T> int NPT<T>::readyCheck() {
402      have_chi_tolerance = 1;
403      painCave.isFatal = 0;
404      simError();
405 <  }
405 >  }
406  
407    if (!have_eta_tolerance) {
408      sprintf( painCave.errMsg,
# Line 332 | Line 411 | template<typename T> int NPT<T>::readyCheck() {
411      have_eta_tolerance = 1;
412      painCave.isFatal = 0;
413      simError();
414 <  }
415 <  
414 >  }
415 >
416    // We need NkBT a lot, so just set it here: This is the RAW number
417    // of particles, so no subtraction or addition of constraints or
418    // orientational degrees of freedom:
419 <  
419 >
420    NkBT = (double)Nparticles * kB * targetTemp;
421 <  
421 >
422    // fkBT is used because the thermostat operates on more degrees of freedom
423    // than the barostat (when there are particles with orientational degrees
424    // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
425 <  
425 >
426    fkBT = (double)info->ndf * kB * targetTemp;
427  
428    tt2 = tauThermostat * tauThermostat;

Diff Legend

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