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 1268 by tim, Fri Jun 11 17:16:21 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];
51 < #ifdef IS_MPI
52 <  Nparticles = mpiSim->getTotAtoms();
45 < #else
46 <  Nparticles = theInfo->n_atoms;
47 < #endif
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 +
69   }
70  
71   template<typename T> NPT<T>::~NPT() {
# Line 55 | Line 75 | template<typename T> void NPT<T>::moveA() {
75   }
76  
77   template<typename T> void NPT<T>::moveA() {
78 <
78 >
79    //new version of NPT
80    int i, j, k;
61 –  DirectionalAtom* dAtom;
81    double Tb[3], ji[3];
82    double mass;
83    double vel[3], pos[3], frc[3];
# Line 69 | Line 88 | template<typename T> void NPT<T>::moveA() {
88    tStats->getPressureTensor( press );
89    instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
90    instaVol = tStats->getVolume();
91 <  
91 >
92    tStats->getCOM(COM);
93 <  
93 >
94    //evolve velocity half step
76 –  for( i=0; i<nAtoms; i++ ){
95  
96 <    atoms[i]->getVel( vel );
97 <    atoms[i]->getFrc( frc );
96 >  calcVelScale();
97 >  
98 >  for( i=0; i<integrableObjects.size(); i++ ){
99  
100 <    mass = atoms[i]->getMass();
100 >    integrableObjects[i]->getVel( vel );
101 >    integrableObjects[i]->getFrc( frc );
102  
103 +    mass = integrableObjects[i]->getMass();
104 +
105      getVelScaleA( sc, vel );
106  
107      for (j=0; j < 3; j++) {
108 <      
108 >
109        // velocity half step  (use chi from previous step here):
110        vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
111 <  
111 >
112      }
113  
114 <    atoms[i]->setVel( vel );
93 <  
94 <    if( atoms[i]->isDirectional() ){
114 >    integrableObjects[i]->setVel( vel );
115  
116 <      dAtom = (DirectionalAtom *)atoms[i];
116 >    if( integrableObjects[i]->isDirectional() ){
117  
118        // get and convert the torque to body frame
119 <      
120 <      dAtom->getTrq( Tb );
121 <      dAtom->lab2Body( Tb );
122 <      
119 >
120 >      integrableObjects[i]->getTrq( Tb );
121 >      integrableObjects[i]->lab2Body( Tb );
122 >
123        // get the angular momentum, and propagate a half step
124  
125 <      dAtom->getJ( ji );
125 >      integrableObjects[i]->getJ( ji );
126  
127 <      for (j=0; j < 3; j++)
127 >      for (j=0; j < 3; j++)
128          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
129 <      
130 <      this->rotationPropagation( dAtom, ji );
131 <      
132 <      dAtom->setJ( ji );
133 <    }    
129 >
130 >      this->rotationPropagation( integrableObjects[i], ji );
131 >
132 >      integrableObjects[i]->setJ( ji );
133 >    }
134    }
135  
136    // evolve chi and eta  half step
137 <  
137 >
138    evolveChiA();
139    evolveEtaA();
140  
# Line 122 | Line 142 | template<typename T> void NPT<T>::moveA() {
142    integralOfChidt += dt2*chi;
143  
144    //save the old positions
145 <  for(i = 0; i < nAtoms; i++){
146 <    atoms[i]->getPos(pos);
145 >  for(i = 0; i < integrableObjects.size(); i++){
146 >    integrableObjects[i]->getPos(pos);
147      for(j = 0; j < 3; j++)
148        oldPos[i*3 + j] = pos[j];
149    }
150 <  
151 <  //the first estimation of r(t+dt) is equal to  r(t)
152 <    
150 >
151 >  //the first estimation of r(t+dt) is equal to  r(t)
152 >
153    for(k = 0; k < 5; k ++){
154  
155 <    for(i =0 ; i < nAtoms; i++){
155 >    for(i =0 ; i < integrableObjects.size(); i++){
156  
157 <      atoms[i]->getVel(vel);
158 <      atoms[i]->getPos(pos);
157 >      integrableObjects[i]->getVel(vel);
158 >      integrableObjects[i]->getPos(pos);
159  
160        this->getPosScale( pos, COM, i, sc );
161 <  
161 >
162        for(j = 0; j < 3; j++)
163          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
164  
165 <      atoms[i]->setPos( pos );
165 >      integrableObjects[i]->setPos( pos );
166      }
167 <    
168 <    if (nConstrained){
149 <      constrainA();
150 <    }
167 >
168 >    consFramework->doConstrainA();
169    }
152 –    
170  
171 +
172    // Scale the box after all the positions have been moved:
173 <  
173 >
174    this->scaleSimBox();
175   }
176  
177   template<typename T> void NPT<T>::moveB( void ){
178 <  
178 >
179    //new version of NPT
180    int i, j, k;
163 –  DirectionalAtom* dAtom;
181    double Tb[3], ji[3], sc[3];
182    double vel[3], frc[3];
183    double mass;
184 <  
184 >
185    // Set things up for the iteration:
186  
187 <  for( i=0; i<nAtoms; i++ ){
187 >  for( i=0; i<integrableObjects.size(); i++ ){
188  
189 <    atoms[i]->getVel( vel );
189 >    integrableObjects[i]->getVel( vel );
190  
191      for (j=0; j < 3; j++)
192        oldVel[3*i + j]  = vel[j];
193  
194 <    if( atoms[i]->isDirectional() ){
194 >    if( integrableObjects[i]->isDirectional() ){
195  
196 <      dAtom = (DirectionalAtom *)atoms[i];
196 >      integrableObjects[i]->getJ( ji );
197  
181 –      dAtom->getJ( ji );
182 –
198        for (j=0; j < 3; j++)
199          oldJi[3*i + j] = ji[j];
200  
# 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 <    
209 >
210      instaTemp = tStats->getTemperature();
211      instaPress = tStats->getPressure();
212  
# Line 199 | Line 214 | template<typename T> void NPT<T>::moveB( void ){
214  
215      this->evolveChiB();
216      this->evolveEtaB();
217 <    
203 <    for( i=0; i<nAtoms; i++ ){
217 >    this->calcVelScale();
218  
219 <      atoms[i]->getFrc( frc );
220 <      atoms[i]->getVel(vel);
221 <      
222 <      mass = atoms[i]->getMass();
223 <      
219 >    for( i=0; i<integrableObjects.size(); i++ ){
220 >
221 >      integrableObjects[i]->getFrc( frc );
222 >      integrableObjects[i]->getVel(vel);
223 >
224 >      mass = integrableObjects[i]->getMass();
225 >
226        getVelScaleB( sc, i );
227  
228        // velocity half step
229 <      for (j=0; j < 3; j++)
229 >      for (j=0; j < 3; j++)
230          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() ){
231  
232 <        dAtom = (DirectionalAtom *)atoms[i];
233 <  
234 <        // get and convert the torque to body frame      
235 <  
236 <        dAtom->getTrq( Tb );
237 <        dAtom->lab2Body( Tb );      
238 <            
239 <        for (j=0; j < 3; j++)
232 >      integrableObjects[i]->setVel( vel );
233 >
234 >      if( integrableObjects[i]->isDirectional() ){
235 >
236 >        // get and convert the torque to body frame
237 >
238 >        integrableObjects[i]->getTrq( Tb );
239 >        integrableObjects[i]->lab2Body( Tb );
240 >
241 >        for (j=0; j < 3; j++)
242            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
243 <      
244 <          dAtom->setJ( ji );
243 >
244 >          integrableObjects[i]->setJ( ji );
245        }
246      }
247 <    
248 <    if (nConstrained){
249 <      constrainB();
236 <    }    
237 <    
247 >
248 >    consFramework->doConstrainA();
249 >
250      if ( this->chiConverged() && this->etaConverged() ) break;
251    }
252  
# Line 255 | Line 267 | template<typename T> void NPT<T>::evolveChiB() {
267   }
268  
269   template<typename T> void NPT<T>::evolveChiB() {
270 <  
270 >
271    prevChi = chi;
272    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
273   }
274  
275   template<typename T> bool NPT<T>::chiConverged() {
276 <  
277 <  return ( fabs( prevChi - chi ) <= chiTolerance );  
276 >
277 >  return ( fabs( prevChi - chi ) <= chiTolerance );
278   }
279  
280   template<typename T> int NPT<T>::readyCheck() {
# Line 270 | Line 282 | template<typename T> int NPT<T>::readyCheck() {
282    //check parent's readyCheck() first
283    if (T::readyCheck() == -1)
284      return -1;
285 <
286 <  // First check to see if we have a target temperature.
287 <  // Not having one is fatal.
288 <  
285 >
286 >  // First check to see if we have a target temperature.
287 >  // Not having one is fatal.
288 >
289    if (!have_target_temp) {
290      sprintf( painCave.errMsg,
291               "NPT error: You can't use the NPT integrator\n"
# Line 293 | Line 305 | template<typename T> int NPT<T>::readyCheck() {
305      simError();
306      return -1;
307    }
308 <  
308 >
309    // We must set tauThermostat.
310 <  
310 >
311    if (!have_tau_thermostat) {
312      sprintf( painCave.errMsg,
313               "NPT error: If you use the NPT\n"
# Line 303 | Line 315 | template<typename T> int NPT<T>::readyCheck() {
315      painCave.isFatal = 1;
316      simError();
317      return -1;
318 <  }    
318 >  }
319  
320    // We must set tauBarostat.
321 <  
321 >
322    if (!have_tau_barostat) {
323      sprintf( painCave.errMsg,
324 <             "NPT error: If you use the NPT\n"
325 <             "   integrator, you must set tauBarostat.\n");
324 >             "If you use the NPT integrator, you must set tauBarostat.\n");
325 >    painCave.severity = OOPSE_ERROR;
326      painCave.isFatal = 1;
327      simError();
328      return -1;
329 <  }    
329 >  }
330  
331    if (!have_chi_tolerance) {
332      sprintf( painCave.errMsg,
333 <             "NPT warning: setting chi tolerance to 1e-6\n");
333 >             "Setting chi tolerance to 1e-6 in NPT integrator\n");
334      chiTolerance = 1e-6;
335      have_chi_tolerance = 1;
336 +    painCave.severity = OOPSE_INFO;
337      painCave.isFatal = 0;
338      simError();
339 <  }
339 >  }
340  
341    if (!have_eta_tolerance) {
342      sprintf( painCave.errMsg,
343 <             "NPT warning: setting eta tolerance to 1e-6\n");
343 >             "Setting eta tolerance to 1e-6 in NPT integrator");
344      etaTolerance = 1e-6;
345      have_eta_tolerance = 1;
346 +    painCave.severity = OOPSE_INFO;
347      painCave.isFatal = 0;
348      simError();
349 <  }
350 <  
349 >  }
350 >
351    // We need NkBT a lot, so just set it here: This is the RAW number
352 <  // of particles, so no subtraction or addition of constraints or
352 >  // of integrableObjects, so no subtraction or addition of constraints or
353    // orientational degrees of freedom:
354 <  
355 <  NkBT = (double)Nparticles * kB * targetTemp;
356 <  
354 >
355 >  NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
356 >
357    // fkBT is used because the thermostat operates on more degrees of freedom
358    // than the barostat (when there are particles with orientational degrees
359 <  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
346 <  
347 <  fkBT = (double)info->ndf * kB * targetTemp;
359 >  // of freedom).  
360  
361 +  fkBT = (double)(info->getNDF()) * kB * targetTemp;
362 +
363    tt2 = tauThermostat * tauThermostat;
364    tb2 = tauBarostat * tauBarostat;
365  

Diff Legend

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