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

Comparing trunk/OOPSE/libmdtools/NPT.cpp (file contents):
Revision 780 by mmeineke, Mon Sep 22 21:23:25 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2 >
3   #include "Atom.hpp"
4   #include "SRI.hpp"
5   #include "AbstractClasses.hpp"
# Line 7 | Line 8
8   #include "Thermo.hpp"
9   #include "ReadWrite.hpp"
10   #include "Integrator.hpp"
11 < #include "simError.h"
11 > #include "simError.h"
12  
13   #ifdef IS_MPI
14   #include "mpiSimulation.hpp"
# Line 18 | Line 19
19   // modification of the Hoover algorithm:
20   //
21   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
22 < //       Molec. Phys., 78, 533.
22 > //       Molec. Phys., 78, 533.
23   //
24   //           and
25 < //
25 > //
26   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
27  
28   template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
29    T( theInfo, the_ff )
30   {
31 +  GenericData* data;
32 +  DoubleData * chiValue;
33 +  DoubleData * integralOfChidtValue;
34 +
35 +  chiValue = NULL;
36 +  integralOfChidtValue = NULL;
37 +
38    chi = 0.0;
39    integralOfChidt = 0.0;
40    have_tau_thermostat = 0;
# Line 37 | Line 45 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
45    have_eta_tolerance = 0;
46    have_pos_iter_tolerance = 0;
47  
48 <  oldPos = new double[3*nAtoms];
49 <  oldVel = new double[3*nAtoms];
50 <  oldJi = new double[3*nAtoms];
48 >  // retrieve chi and integralOfChidt from simInfo
49 >  data = info->getProperty(CHIVALUE_ID);
50 >  if(data){
51 >    chiValue = dynamic_cast<DoubleData*>(data);
52 >  }
53 >
54 >  data = info->getProperty(INTEGRALOFCHIDT_ID);
55 >  if(data){
56 >    integralOfChidtValue = dynamic_cast<DoubleData*>(data);
57 >  }
58 >
59 >  // chi and integralOfChidt should appear by pair
60 >  if(chiValue && integralOfChidtValue){
61 >    chi = chiValue->getData();
62 >    integralOfChidt = integralOfChidtValue->getData();
63 >  }
64 >
65 >  oldPos = new double[3*integrableObjects.size()];
66 >  oldVel = new double[3*integrableObjects.size()];
67 >  oldJi = new double[3*integrableObjects.size()];
68   #ifdef IS_MPI
69    Nparticles = mpiSim->getTotAtoms();
70   #else
# Line 55 | Line 80 | template<typename T> void NPT<T>::moveA() {
80   }
81  
82   template<typename T> void NPT<T>::moveA() {
83 <
83 >
84    //new version of NPT
85    int i, j, k;
61 –  DirectionalAtom* dAtom;
86    double Tb[3], ji[3];
87    double mass;
88    double vel[3], pos[3], frc[3];
# Line 69 | Line 93 | template<typename T> void NPT<T>::moveA() {
93    tStats->getPressureTensor( press );
94    instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
95    instaVol = tStats->getVolume();
96 <  
96 >
97    tStats->getCOM(COM);
98 <  
98 >
99    //evolve velocity half step
76 –  for( i=0; i<nAtoms; i++ ){
100  
101 <    atoms[i]->getVel( vel );
102 <    atoms[i]->getFrc( frc );
101 >  calcVelScale();
102 >  
103 >  for( i=0; i<integrableObjects.size(); i++ ){
104  
105 <    mass = atoms[i]->getMass();
105 >    integrableObjects[i]->getVel( vel );
106 >    integrableObjects[i]->getFrc( frc );
107  
108 +    mass = integrableObjects[i]->getMass();
109 +
110      getVelScaleA( sc, vel );
111  
112      for (j=0; j < 3; j++) {
113 <      
113 >
114        // velocity half step  (use chi from previous step here):
115        vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
116 <  
116 >
117      }
118  
119 <    atoms[i]->setVel( vel );
93 <  
94 <    if( atoms[i]->isDirectional() ){
119 >    integrableObjects[i]->setVel( vel );
120  
121 <      dAtom = (DirectionalAtom *)atoms[i];
121 >    if( integrableObjects[i]->isDirectional() ){
122  
123        // get and convert the torque to body frame
124 <      
125 <      dAtom->getTrq( Tb );
126 <      dAtom->lab2Body( Tb );
127 <      
124 >
125 >      integrableObjects[i]->getTrq( Tb );
126 >      integrableObjects[i]->lab2Body( Tb );
127 >
128        // get the angular momentum, and propagate a half step
129  
130 <      dAtom->getJ( ji );
130 >      integrableObjects[i]->getJ( ji );
131  
132 <      for (j=0; j < 3; j++)
132 >      for (j=0; j < 3; j++)
133          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
134 <      
135 <      this->rotationPropagation( dAtom, ji );
136 <      
137 <      dAtom->setJ( ji );
138 <    }    
134 >
135 >      this->rotationPropagation( integrableObjects[i], ji );
136 >
137 >      integrableObjects[i]->setJ( ji );
138 >    }
139    }
140  
141    // evolve chi and eta  half step
142 <  
142 >
143    evolveChiA();
144    evolveEtaA();
145  
# Line 122 | Line 147 | template<typename T> void NPT<T>::moveA() {
147    integralOfChidt += dt2*chi;
148  
149    //save the old positions
150 <  for(i = 0; i < nAtoms; i++){
151 <    atoms[i]->getPos(pos);
150 >  for(i = 0; i < integrableObjects.size(); i++){
151 >    integrableObjects[i]->getPos(pos);
152      for(j = 0; j < 3; j++)
153        oldPos[i*3 + j] = pos[j];
154    }
155 <  
156 <  //the first estimation of r(t+dt) is equal to  r(t)
157 <    
155 >
156 >  //the first estimation of r(t+dt) is equal to  r(t)
157 >
158    for(k = 0; k < 5; k ++){
159  
160 <    for(i =0 ; i < nAtoms; i++){
160 >    for(i =0 ; i < integrableObjects.size(); i++){
161  
162 <      atoms[i]->getVel(vel);
163 <      atoms[i]->getPos(pos);
162 >      integrableObjects[i]->getVel(vel);
163 >      integrableObjects[i]->getPos(pos);
164  
165        this->getPosScale( pos, COM, i, sc );
166 <  
166 >
167        for(j = 0; j < 3; j++)
168          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
169  
170 <      atoms[i]->setPos( pos );
170 >      integrableObjects[i]->setPos( pos );
171      }
172 <    
172 >
173      if (nConstrained){
174        constrainA();
175      }
176    }
152 –    
177  
178 +
179    // Scale the box after all the positions have been moved:
180 <  
180 >
181    this->scaleSimBox();
182   }
183  
184   template<typename T> void NPT<T>::moveB( void ){
185 <  
185 >
186    //new version of NPT
187    int i, j, k;
163 –  DirectionalAtom* dAtom;
188    double Tb[3], ji[3], sc[3];
189    double vel[3], frc[3];
190    double mass;
191 <  
191 >
192    // Set things up for the iteration:
193  
194 <  for( i=0; i<nAtoms; i++ ){
194 >  for( i=0; i<integrableObjects.size(); i++ ){
195  
196 <    atoms[i]->getVel( vel );
196 >    integrableObjects[i]->getVel( vel );
197  
198      for (j=0; j < 3; j++)
199        oldVel[3*i + j]  = vel[j];
200  
201 <    if( atoms[i]->isDirectional() ){
201 >    if( integrableObjects[i]->isDirectional() ){
202  
203 <      dAtom = (DirectionalAtom *)atoms[i];
203 >      integrableObjects[i]->getJ( ji );
204  
181 –      dAtom->getJ( ji );
182 –
205        for (j=0; j < 3; j++)
206          oldJi[3*i + j] = ji[j];
207  
# Line 189 | Line 211 | template<typename T> void NPT<T>::moveB( void ){
211    // do the iteration:
212  
213    instaVol = tStats->getVolume();
214 <  
214 >
215    for (k=0; k < 4; k++) {
216 <    
216 >
217      instaTemp = tStats->getTemperature();
218      instaPress = tStats->getPressure();
219  
# Line 199 | Line 221 | template<typename T> void NPT<T>::moveB( void ){
221  
222      this->evolveChiB();
223      this->evolveEtaB();
224 <    
203 <    for( i=0; i<nAtoms; i++ ){
224 >    this->calcVelScale();
225  
226 <      atoms[i]->getFrc( frc );
227 <      atoms[i]->getVel(vel);
228 <      
229 <      mass = atoms[i]->getMass();
230 <      
226 >    for( i=0; i<integrableObjects.size(); i++ ){
227 >
228 >      integrableObjects[i]->getFrc( frc );
229 >      integrableObjects[i]->getVel(vel);
230 >
231 >      mass = integrableObjects[i]->getMass();
232 >
233        getVelScaleB( sc, i );
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]);
215 –      
216 –      atoms[i]->setVel( vel );
217 –      
218 –      if( atoms[i]->isDirectional() ){
238  
239 <        dAtom = (DirectionalAtom *)atoms[i];
240 <  
241 <        // get and convert the torque to body frame      
242 <  
243 <        dAtom->getTrq( Tb );
244 <        dAtom->lab2Body( Tb );      
245 <            
246 <        for (j=0; j < 3; j++)
239 >      integrableObjects[i]->setVel( vel );
240 >
241 >      if( integrableObjects[i]->isDirectional() ){
242 >
243 >        // get and convert the torque to body frame
244 >
245 >        integrableObjects[i]->getTrq( Tb );
246 >        integrableObjects[i]->lab2Body( Tb );
247 >
248 >        for (j=0; j < 3; j++)
249            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
250 <      
251 <          dAtom->setJ( ji );
250 >
251 >          integrableObjects[i]->setJ( ji );
252        }
253      }
254 <    
254 >
255      if (nConstrained){
256        constrainB();
257 <    }    
258 <    
257 >    }
258 >
259      if ( this->chiConverged() && this->etaConverged() ) break;
260    }
261  
# Line 255 | Line 276 | template<typename T> void NPT<T>::evolveChiB() {
276   }
277  
278   template<typename T> void NPT<T>::evolveChiB() {
279 <  
279 >
280    prevChi = chi;
281    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
282   }
283  
284   template<typename T> bool NPT<T>::chiConverged() {
285 <  
286 <  return ( fabs( prevChi - chi ) <= chiTolerance );  
285 >
286 >  return ( fabs( prevChi - chi ) <= chiTolerance );
287   }
288  
289   template<typename T> int NPT<T>::readyCheck() {
# Line 270 | Line 291 | template<typename T> int NPT<T>::readyCheck() {
291    //check parent's readyCheck() first
292    if (T::readyCheck() == -1)
293      return -1;
294 <
295 <  // First check to see if we have a target temperature.
296 <  // Not having one is fatal.
297 <  
294 >
295 >  // First check to see if we have a target temperature.
296 >  // Not having one is fatal.
297 >
298    if (!have_target_temp) {
299      sprintf( painCave.errMsg,
300               "NPT error: You can't use the NPT integrator\n"
# Line 293 | Line 314 | template<typename T> int NPT<T>::readyCheck() {
314      simError();
315      return -1;
316    }
317 <  
317 >
318    // We must set tauThermostat.
319 <  
319 >
320    if (!have_tau_thermostat) {
321      sprintf( painCave.errMsg,
322               "NPT error: If you use the NPT\n"
# Line 303 | Line 324 | template<typename T> int NPT<T>::readyCheck() {
324      painCave.isFatal = 1;
325      simError();
326      return -1;
327 <  }    
327 >  }
328  
329    // We must set tauBarostat.
330 <  
330 >
331    if (!have_tau_barostat) {
332      sprintf( painCave.errMsg,
333               "NPT error: If you use the NPT\n"
# Line 314 | Line 335 | template<typename T> int NPT<T>::readyCheck() {
335      painCave.isFatal = 1;
336      simError();
337      return -1;
338 <  }    
338 >  }
339  
340    if (!have_chi_tolerance) {
341      sprintf( painCave.errMsg,
# Line 323 | Line 344 | template<typename T> int NPT<T>::readyCheck() {
344      have_chi_tolerance = 1;
345      painCave.isFatal = 0;
346      simError();
347 <  }
347 >  }
348  
349    if (!have_eta_tolerance) {
350      sprintf( painCave.errMsg,
# Line 332 | Line 353 | template<typename T> int NPT<T>::readyCheck() {
353      have_eta_tolerance = 1;
354      painCave.isFatal = 0;
355      simError();
356 <  }
357 <  
356 >  }
357 >
358    // We need NkBT a lot, so just set it here: This is the RAW number
359    // of particles, so no subtraction or addition of constraints or
360    // orientational degrees of freedom:
361 <  
361 >
362    NkBT = (double)Nparticles * kB * targetTemp;
363 <  
363 >
364    // fkBT is used because the thermostat operates on more degrees of freedom
365    // than the barostat (when there are particles with orientational degrees
366    // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
367 <  
367 >
368    fkBT = (double)info->ndf * kB * targetTemp;
369  
370    tt2 = tauThermostat * tauThermostat;

Diff Legend

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