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 852 by mmeineke, Thu Nov 6 18:20:47 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 66 | Line 81 | template<typename T> void NPT<T>::moveA() {
81    double COM[3];
82  
83    instaTemp = tStats->getTemperature();
84 +
85    tStats->getPressureTensor( press );
86    instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
87 +
88    instaVol = tStats->getVolume();
89 <  
89 >
90    tStats->getCOM(COM);
91 <  
91 >
92    //evolve velocity half step
93    for( i=0; i<nAtoms; i++ ){
94  
# Line 83 | Line 100 | template<typename T> void NPT<T>::moveA() {
100      getVelScaleA( sc, vel );
101  
102      for (j=0; j < 3; j++) {
103 <      
103 >
104        // velocity half step  (use chi from previous step here):
105        vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
106 <  
106 >
107      }
108  
109      atoms[i]->setVel( vel );
110 <  
110 >
111      if( atoms[i]->isDirectional() ){
112  
113        dAtom = (DirectionalAtom *)atoms[i];
114  
115        // get and convert the torque to body frame
116 <      
116 >
117        dAtom->getTrq( Tb );
118        dAtom->lab2Body( Tb );
119 <      
119 >
120        // get the angular momentum, and propagate a half step
121  
122        dAtom->getJ( ji );
123  
124 <      for (j=0; j < 3; j++)
124 >      for (j=0; j < 3; j++)
125          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
126 <      
126 >
127        this->rotationPropagation( dAtom, ji );
128 <      
128 >
129        dAtom->setJ( ji );
130 <    }    
130 >    }
131    }
132  
133    // evolve chi and eta  half step
134 <  
134 >
135    evolveChiA();
136    evolveEtaA();
137  
# Line 127 | Line 144 | template<typename T> void NPT<T>::moveA() {
144      for(j = 0; j < 3; j++)
145        oldPos[i*3 + j] = pos[j];
146    }
147 <  
148 <  //the first estimation of r(t+dt) is equal to  r(t)
149 <    
147 >
148 >  //the first estimation of r(t+dt) is equal to  r(t)
149 >
150    for(k = 0; k < 5; k ++){
151  
152      for(i =0 ; i < nAtoms; i++){
# Line 138 | Line 155 | template<typename T> void NPT<T>::moveA() {
155        atoms[i]->getPos(pos);
156  
157        this->getPosScale( pos, COM, i, sc );
158 <  
158 >
159        for(j = 0; j < 3; j++)
160          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
161  
162        atoms[i]->setPos( pos );
163      }
164 <    
164 >
165      if (nConstrained){
166        constrainA();
167      }
168    }
152 –    
169  
170 +
171    // Scale the box after all the positions have been moved:
172 <  
172 >
173    this->scaleSimBox();
174   }
175  
176 < template<typename T> void NPT<T>::moveB( void ){
177 <  
176 > void NPT::moveB( void ){
177 >
178    //new version of NPT
179    int i, j, k;
180    DirectionalAtom* dAtom;
181    double Tb[3], ji[3], sc[3];
182 <  double vel[3], frc[3];
182 >  double vel[3], frc[3], qVel[3];
183    double mass;
184 <  
184 >
185    // Set things up for the iteration:
186  
187    for( i=0; i<nAtoms; i++ ){
# Line 189 | Line 206 | template<typename T> void NPT<T>::moveB( void ){
206    // do the iteration:
207  
208    instaVol = tStats->getVolume();
209 <  
209 >
210    for (k=0; k < 4; k++) {
211 <    
211 >
212 >    atoms[0]->getVel(vel);
213 >
214      instaTemp = tStats->getTemperature();
215 +    
216      instaPress = tStats->getPressure();
217  
218      // evolve chi another half step using the temperature at t + dt/2
219  
220      this->evolveChiB();
221      this->evolveEtaB();
222 <    
222 >
223      for( i=0; i<nAtoms; i++ ){
224  
225        atoms[i]->getFrc( frc );
226        atoms[i]->getVel(vel);
227 <      
227 >
228        mass = atoms[i]->getMass();
229 <      
229 >
230        getVelScaleB( sc, i );
231  
232 +      for(j=0;j<3;j++)
233 +        qVel[j] = vel[j];
234 +
235        // velocity half step
236 <      for (j=0; j < 3; j++)
236 >      for (j=0; j < 3; j++)
237          vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
238 <      
238 >
239        atoms[i]->setVel( vel );
240 <      
240 >
241        if( atoms[i]->isDirectional() ){
242  
243          dAtom = (DirectionalAtom *)atoms[i];
244 <  
245 <        // get and convert the torque to body frame      
246 <  
244 >
245 >        // get and convert the torque to body frame
246 >
247          dAtom->getTrq( Tb );
248 <        dAtom->lab2Body( Tb );      
249 <            
250 <        for (j=0; j < 3; j++)
248 >        dAtom->lab2Body( Tb );
249 >
250 >        for (j=0; j < 3; j++)
251            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
252 <      
252 >
253            dAtom->setJ( ji );
254        }
255      }
256 <    
256 >
257      if (nConstrained){
258        constrainB();
259 <    }    
260 <    
259 >    }
260 >
261      if ( this->chiConverged() && this->etaConverged() ) break;
262    }
263  
# Line 244 | Line 267 | template<typename T> void NPT<T>::moveB( void ){
267  
268   }
269  
270 < template<typename T> void NPT<T>::resetIntegrator() {
270 > void NPT::resetIntegrator() {
271    chi = 0.0;
272 <  T::resetIntegrator();
272 >  Integrator::resetIntegrator();
273   }
274  
275 < template<typename T> void NPT<T>::evolveChiA() {
275 > void NPT::evolveChiA() {
276 >
277    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
278    oldChi = chi;
279   }
280  
281 < template<typename T> void NPT<T>::evolveChiB() {
282 <  
281 > void NPT::evolveChiB() {
282 >
283    prevChi = chi;
284    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
285   }
286  
287 < template<typename T> bool NPT<T>::chiConverged() {
288 <  
289 <  return ( fabs( prevChi - chi ) <= chiTolerance );  
287 > bool NPT::chiConverged() {
288 >
289 >  return ( fabs( prevChi - chi ) <= chiTolerance );
290   }
291  
292 < template<typename T> int NPT<T>::readyCheck() {
292 > int NPT::readyCheck() {
293  
294    //check parent's readyCheck() first
295 <  if (T::readyCheck() == -1)
295 >  if (Integrator::readyCheck() == -1)
296      return -1;
297 <
298 <  // First check to see if we have a target temperature.
299 <  // Not having one is fatal.
300 <  
297 >
298 >  // First check to see if we have a target temperature.
299 >  // Not having one is fatal.
300 >
301    if (!have_target_temp) {
302      sprintf( painCave.errMsg,
303               "NPT error: You can't use the NPT integrator\n"
# Line 293 | Line 317 | template<typename T> int NPT<T>::readyCheck() {
317      simError();
318      return -1;
319    }
320 <  
320 >
321    // We must set tauThermostat.
322 <  
322 >
323    if (!have_tau_thermostat) {
324      sprintf( painCave.errMsg,
325               "NPT error: If you use the NPT\n"
# Line 303 | Line 327 | template<typename T> int NPT<T>::readyCheck() {
327      painCave.isFatal = 1;
328      simError();
329      return -1;
330 <  }    
330 >  }
331  
332    // We must set tauBarostat.
333 <  
333 >
334    if (!have_tau_barostat) {
335      sprintf( painCave.errMsg,
336               "NPT error: If you use the NPT\n"
# Line 314 | Line 338 | template<typename T> int NPT<T>::readyCheck() {
338      painCave.isFatal = 1;
339      simError();
340      return -1;
341 <  }    
341 >  }
342  
343    if (!have_chi_tolerance) {
344      sprintf( painCave.errMsg,
# Line 323 | Line 347 | template<typename T> int NPT<T>::readyCheck() {
347      have_chi_tolerance = 1;
348      painCave.isFatal = 0;
349      simError();
350 <  }
350 >  }
351  
352    if (!have_eta_tolerance) {
353      sprintf( painCave.errMsg,
# Line 332 | Line 356 | template<typename T> int NPT<T>::readyCheck() {
356      have_eta_tolerance = 1;
357      painCave.isFatal = 0;
358      simError();
359 <  }
360 <  
359 >  }
360 >
361    // We need NkBT a lot, so just set it here: This is the RAW number
362    // of particles, so no subtraction or addition of constraints or
363    // orientational degrees of freedom:
364 <  
364 >
365    NkBT = (double)Nparticles * kB * targetTemp;
366 <  
366 >
367    // fkBT is used because the thermostat operates on more degrees of freedom
368    // than the barostat (when there are particles with orientational degrees
369    // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
370 <  
370 >
371    fkBT = (double)info->ndf * kB * targetTemp;
372  
373    tt2 = tauThermostat * tauThermostat;

Diff Legend

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