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

Comparing trunk/OOPSE/libmdtools/NPTf.cpp (file contents):
Revision 576 by gezelter, Tue Jul 8 21:10:16 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 1 | Line 1
1 + #include <math.h>
2 +
3 + #include "MatVec3.h"
4   #include "Atom.hpp"
5   #include "SRI.hpp"
6   #include "AbstractClasses.hpp"
# Line 6 | 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"
16 + #endif
17  
18 < // Basic isotropic thermostating and barostating via the Melchionna
18 > // Basic non-isotropic thermostating and barostating via the Melchionna
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 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
29 <  Integrator( theInfo, the_ff )
28 > template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
29 >  T( theInfo, the_ff )
30   {
31 <  int i;
32 <  chi = 0.0;
33 <  for(i = 0; i < 9; i++) eta[i] = 0.0;
34 <  have_tau_thermostat = 0;
29 <  have_tau_barostat = 0;
30 <  have_target_temp = 0;
31 <  have_target_pressure = 0;
32 < }
31 >  GenericData* data;
32 >  DoubleArrayData * etaValue;
33 >  vector<double> etaArray;
34 >  int i,j;
35  
36 < void NPTi::moveA() {
37 <  
36 <  int i,j,k;
37 <  int atomIndex, aMatIndex;
38 <  DirectionalAtom* dAtom;
39 <  double Tb[3];
40 <  double ji[3];
41 <  double rj[3];
42 <  double instaTemp, instaPress, instaVol;
43 <  double tt2, tb2;
44 <  double angle;
36 >  for(i = 0; i < 3; i++){
37 >    for (j = 0; j < 3; j++){
38  
39 <  tt2 = tauThermostat * tauThermostat;
40 <  tb2 = tauBarostat * tauBarostat;
39 >      eta[i][j] = 0.0;
40 >      oldEta[i][j] = 0.0;
41 >    }
42 >  }
43  
44 <  instaTemp = tStats->getTemperature();
45 <  instaPress = tStats->getPressure();
46 <  instaVol = tStats->getVolume();
47 <  
48 <  // first evolve chi a half step
49 <  
50 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
51 <  
52 <  for (i = 0; i < 9; i++) {
53 <    eta[i] += dt2 * ( instaVol * (sigma[i] - targetPressure*identMat[i]))
54 <      / (NkBT*tb2));
44 >
45 >  if( theInfo->useInitXSstate ){
46 >    // retrieve eta array from simInfo if it exists
47 >    data = info->getProperty(ETAVALUE_ID);
48 >    if(data){
49 >      etaValue = dynamic_cast<DoubleArrayData*>(data);
50 >      
51 >      if(etaValue){
52 >        etaArray = etaValue->getData();
53 >        
54 >        for(i = 0; i < 3; i++){
55 >          for (j = 0; j < 3; j++){
56 >            eta[i][j] = etaArray[3*i+j];
57 >            oldEta[i][j] = eta[i][j];
58 >          }
59 >        }
60 >      }
61 >    }
62 >  }
63 >
64   }
65  
66 <  for( i=0; i<nAtoms; i++ ){
63 <    atomIndex = i * 3;
64 <    aMatIndex = i * 9;
65 <    
66 <    // velocity half step
67 <    for( j=atomIndex; j<(atomIndex+3); j++ )
68 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
69 <                       - vel[j]*(chi+eta));
66 > template<typename T> NPTf<T>::~NPTf() {
67  
68 <    // position whole step    
68 >  // empty for now
69 > }
70  
71 <    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
74 <      rj[0] = pos[j];
75 <      rj[1] = pos[j+1];
76 <      rj[2] = pos[j+2];
71 > template<typename T> void NPTf<T>::resetIntegrator() {
72  
73 <      info->wrapVector(rj);
73 >  int i, j;
74  
75 <      pos[j] += dt * (vel[j] + eta*rj[0]);
76 <      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
77 <      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
75 >  for(i = 0; i < 3; i++)
76 >    for (j = 0; j < 3; j++)
77 >      eta[i][j] = 0.0;
78 >
79 >  T::resetIntegrator();
80 > }
81 >
82 > template<typename T> void NPTf<T>::evolveEtaA() {
83 >
84 >  int i, j;
85 >
86 >  for(i = 0; i < 3; i ++){
87 >    for(j = 0; j < 3; j++){
88 >      if( i == j)
89 >        eta[i][j] += dt2 *  instaVol *
90 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
91 >      else
92 >        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
93      }
94 +  }
95  
96 <    // Scale the box after all the positions have been moved:
96 >  for(i = 0; i < 3; i++)
97 >    for (j = 0; j < 3; j++)
98 >      oldEta[i][j] = eta[i][j];
99 > }
100  
101 <    info->scaleBox(exp(dt*eta));
88 <  
89 <    if( atoms[i]->isDirectional() ){
101 > template<typename T> void NPTf<T>::evolveEtaB() {
102  
103 <      dAtom = (DirectionalAtom *)atoms[i];
92 <          
93 <      // get and convert the torque to body frame
94 <      
95 <      Tb[0] = dAtom->getTx();
96 <      Tb[1] = dAtom->getTy();
97 <      Tb[2] = dAtom->getTz();
98 <      
99 <      dAtom->lab2Body( Tb );
100 <      
101 <      // get the angular momentum, and propagate a half step
103 >  int i,j;
104  
105 <      ji[0] = dAtom->getJx();
106 <      ji[1] = dAtom->getJy();
107 <      ji[2] = dAtom->getJz();
108 <      
109 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
110 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
111 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
112 <      
113 <      // use the angular velocities to propagate the rotation matrix a
114 <      // full time step
115 <      
116 <      // rotate about the x-axis      
115 <      angle = dt2 * ji[0] / dAtom->getIxx();
116 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
117 <      
118 <      // rotate about the y-axis
119 <      angle = dt2 * ji[1] / dAtom->getIyy();
120 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
121 <      
122 <      // rotate about the z-axis
123 <      angle = dt * ji[2] / dAtom->getIzz();
124 <      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
125 <      
126 <      // rotate about the y-axis
127 <      angle = dt2 * ji[1] / dAtom->getIyy();
128 <      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
129 <      
130 <       // rotate about the x-axis
131 <      angle = dt2 * ji[0] / dAtom->getIxx();
132 <      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
133 <      
134 <      dAtom->setJx( ji[0] );
135 <      dAtom->setJy( ji[1] );
136 <      dAtom->setJz( ji[2] );
105 >  for(i = 0; i < 3; i++)
106 >    for (j = 0; j < 3; j++)
107 >      prevEta[i][j] = eta[i][j];
108 >
109 >  for(i = 0; i < 3; i ++){
110 >    for(j = 0; j < 3; j++){
111 >      if( i == j) {
112 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
113 >          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
114 >      } else {
115 >        eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
116 >      }
117      }
138 –    
118    }
119   }
120  
121 < void NPTi::moveB( void ){
122 <  int i,j,k;
144 <  int atomIndex;
145 <  DirectionalAtom* dAtom;
146 <  double Tb[3];
147 <  double ji[3];
148 <  double instaTemp, instaPress, instaVol;
149 <  double tt2, tb2;
150 <  
151 <  tt2 = tauThermostat * tauThermostat;
152 <  tb2 = tauBarostat * tauBarostat;
121 > template<typename T> void NPTf<T>::calcVelScale(void){
122 >  int i,j;
123  
124 <  instaTemp = tStats->getTemperature();
125 <  instaPress = tStats->getPressure();
126 <  instaVol = tStats->getVolume();
124 >  for (i = 0; i < 3; i++ ) {
125 >    for (j = 0; j < 3; j++ ) {
126 >      vScale[i][j] = eta[i][j];
127  
128 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
129 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
130 <  
161 <  for( i=0; i<nAtoms; i++ ){
162 <    atomIndex = i * 3;
163 <    
164 <    // velocity half step
165 <    for( j=atomIndex; j<(atomIndex+3); j++ )
166 <    for( j=atomIndex; j<(atomIndex+3); j++ )
167 <      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168 <                       - vel[j]*(chi+eta));
169 <    
170 <    if( atoms[i]->isDirectional() ){
171 <      
172 <      dAtom = (DirectionalAtom *)atoms[i];
173 <      
174 <      // get and convert the torque to body frame
175 <      
176 <      Tb[0] = dAtom->getTx();
177 <      Tb[1] = dAtom->getTy();
178 <      Tb[2] = dAtom->getTz();
179 <      
180 <      dAtom->lab2Body( Tb );
181 <      
182 <      // get the angular momentum, and complete the angular momentum
183 <      // half step
184 <      
185 <      ji[0] = dAtom->getJx();
186 <      ji[1] = dAtom->getJy();
187 <      ji[2] = dAtom->getJz();
188 <      
189 <      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
190 <      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
191 <      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
192 <      
193 <      dAtom->setJx( ji[0] );
194 <      dAtom->setJy( ji[1] );
195 <      dAtom->setJz( ji[2] );
128 >      if (i == j) {
129 >        vScale[i][j] += chi;
130 >      }
131      }
132    }
133   }
134  
135 < int NPTi::readyCheck() {
135 > template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
136  
137 <  // First check to see if we have a target temperature.
138 <  // Not having one is fatal.
139 <  
140 <  if (!have_target_temp) {
141 <    sprintf( painCave.errMsg,
142 <             "NPTi error: You can't use the NPTi integrator\n"
143 <             "   without a targetTemp!\n"
144 <             );
145 <    painCave.isFatal = 1;
146 <    simError();
147 <    return -1;
137 >  matVecMul3( vScale, vel, sc );
138 > }
139 >
140 > template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
141 >  int j;
142 >  double myVel[3];
143 >  double vScale[3][3];
144 >
145 >  for (j = 0; j < 3; j++)
146 >    myVel[j] = oldVel[3*index + j];
147 >
148 >  matVecMul3( vScale, myVel, sc );
149 > }
150 >
151 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
152 >                                               int index, double sc[3]){
153 >  int j;
154 >  double rj[3];
155 >
156 >  for(j=0; j<3; j++)
157 >    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
158 >
159 >  matVecMul3( eta, rj, sc );
160 > }
161 >
162 > template<typename T> void NPTf<T>::scaleSimBox( void ){
163 >
164 >  int i,j,k;
165 >  double scaleMat[3][3];
166 >  double eta2ij;
167 >  double bigScale, smallScale, offDiagMax;
168 >  double hm[3][3], hmnew[3][3];
169 >
170 >
171 >
172 >  // Scale the box after all the positions have been moved:
173 >
174 >  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
175 >  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
176 >
177 >  bigScale = 1.0;
178 >  smallScale = 1.0;
179 >  offDiagMax = 0.0;
180 >
181 >  for(i=0; i<3; i++){
182 >    for(j=0; j<3; j++){
183 >
184 >      // Calculate the matrix Product of the eta array (we only need
185 >      // the ij element right now):
186 >
187 >      eta2ij = 0.0;
188 >      for(k=0; k<3; k++){
189 >        eta2ij += eta[i][k] * eta[k][j];
190 >      }
191 >
192 >      scaleMat[i][j] = 0.0;
193 >      // identity matrix (see above):
194 >      if (i == j) scaleMat[i][j] = 1.0;
195 >      // Taylor expansion for the exponential truncated at second order:
196 >      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
197 >
198 >      if (i != j)
199 >        if (fabs(scaleMat[i][j]) > offDiagMax)
200 >          offDiagMax = fabs(scaleMat[i][j]);
201 >    }
202 >
203 >    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
204 >    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
205    }
206  
207 <  if (!have_target_pressure) {
207 >  if ((bigScale > 1.01) || (smallScale < 0.99)) {
208      sprintf( painCave.errMsg,
209 <             "NPTi error: You can't use the NPTi integrator\n"
210 <             "   without a targetPressure!\n"
211 <             );
209 >             "NPTf error: Attempting a Box scaling of more than 1 percent.\n"
210 >             " Check your tauBarostat, as it is probably too small!\n\n"
211 >             " scaleMat = [%lf\t%lf\t%lf]\n"
212 >             "            [%lf\t%lf\t%lf]\n"
213 >             "            [%lf\t%lf\t%lf]\n",
214 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
215 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
216 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
217      painCave.isFatal = 1;
218      simError();
219 <    return -1;
223 <  }
224 <  
225 <  // We must set tauThermostat.
226 <  
227 <  if (!have_tau_thermostat) {
219 >  } else if (offDiagMax > 0.01) {
220      sprintf( painCave.errMsg,
221 <             "NPTi error: If you use the NPTi\n"
222 <             "   integrator, you must set tauThermostat.\n");
221 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
222 >             " Check your tauBarostat, as it is probably too small!\n\n"
223 >             " scaleMat = [%lf\t%lf\t%lf]\n"
224 >             "            [%lf\t%lf\t%lf]\n"
225 >             "            [%lf\t%lf\t%lf]\n",
226 >             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
227 >             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
228 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
229      painCave.isFatal = 1;
230      simError();
231 <    return -1;
232 <  }    
231 >  } else {
232 >    info->getBoxM(hm);
233 >    matMul3(hm, scaleMat, hmnew);
234 >    info->setBoxM(hmnew);
235 >  }
236 > }
237  
238 <  // We must set tauBarostat.
239 <  
240 <  if (!have_tau_barostat) {
239 <    sprintf( painCave.errMsg,
240 <             "NPTi error: If you use the NPTi\n"
241 <             "   integrator, you must set tauBarostat.\n");
242 <    painCave.isFatal = 1;
243 <    simError();
244 <    return -1;
245 <  }    
238 > template<typename T> bool NPTf<T>::etaConverged() {
239 >  int i;
240 >  double diffEta, sumEta;
241  
242 <  // We need NkBT a lot, so just set it here:
242 >  sumEta = 0;
243 >  for(i = 0; i < 3; i++)
244 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
245  
246 <  NkBT = (double)info->ndf * kB * targetTemp;
246 >  diffEta = sqrt( sumEta / 3.0 );
247  
248 <  return 1;
248 >  return ( diffEta <= etaTolerance );
249   }
250 +
251 + template<typename T> double NPTf<T>::getConservedQuantity(void){
252 +
253 +  double conservedQuantity;
254 +  double totalEnergy;
255 +  double thermostat_kinetic;
256 +  double thermostat_potential;
257 +  double barostat_kinetic;
258 +  double barostat_potential;
259 +  double trEta;
260 +  double a[3][3], b[3][3];
261 +
262 +  totalEnergy = tStats->getTotalE();
263 +
264 +  thermostat_kinetic = fkBT * tt2 * chi * chi /
265 +    (2.0 * eConvert);
266 +
267 +  thermostat_potential = fkBT* integralOfChidt / eConvert;
268 +
269 +  transposeMat3(eta, a);
270 +  matMul3(a, eta, b);
271 +  trEta = matTrace3(b);
272 +
273 +  barostat_kinetic = NkBT * tb2 * trEta /
274 +    (2.0 * eConvert);
275 +
276 +  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
277 +    eConvert;
278 +
279 +  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
280 +    barostat_kinetic + barostat_potential;
281 +
282 +  return conservedQuantity;
283 +
284 + }
285 +
286 + template<typename T> string NPTf<T>::getAdditionalParameters(void){
287 +  string parameters;
288 +  const int BUFFERSIZE = 2000; // size of the read buffer
289 +  char buffer[BUFFERSIZE];
290 +
291 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
292 +  parameters += buffer;
293 +
294 +  for(int i = 0; i < 3; i++){
295 +    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
296 +    parameters += buffer;
297 +  }
298 +
299 +  return parameters;
300 +
301 + }

Diff Legend

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