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 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 1253 by gezelter, Tue Jun 8 16:49:46 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 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 17 | 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> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
29    T( theInfo, the_ff )
30   {
31 <  
31 >  GenericData* data;
32 >  DoubleArrayData * etaValue;
33 >  vector<double> etaArray;
34    int i,j;
35 <  
35 >
36    for(i = 0; i < 3; i++){
37      for (j = 0; j < 3; j++){
38 <      
38 >
39        eta[i][j] = 0.0;
40        oldEta[i][j] = 0.0;
41      }
42    }
43 +
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   template<typename T> NPTf<T>::~NPTf() {
# Line 44 | Line 69 | template<typename T> void NPTf<T>::resetIntegrator() {
69   }
70  
71   template<typename T> void NPTf<T>::resetIntegrator() {
72 <  
72 >
73    int i, j;
74 <  
74 >
75    for(i = 0; i < 3; i++)
76      for (j = 0; j < 3; j++)
77        eta[i][j] = 0.0;
78 <  
78 >
79    T::resetIntegrator();
80   }
81  
82   template<typename T> void NPTf<T>::evolveEtaA() {
83 <  
83 >
84    int i, j;
85 <  
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 *
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);
# Line 74 | Line 99 | template<typename T> void NPTf<T>::evolveEtaB() {
99   }
100  
101   template<typename T> void NPTf<T>::evolveEtaB() {
102 <  
102 >
103    int i,j;
104  
105    for(i = 0; i < 3; i++)
# Line 84 | Line 109 | template<typename T> void NPTf<T>::evolveEtaB() {
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 *
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);
# Line 93 | Line 118 | template<typename T> void NPTf<T>::evolveEtaB() {
118    }
119   }
120  
121 < template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
121 > template<typename T> void NPTf<T>::calcVelScale(void){
122    int i,j;
98 –  double vScale[3][3];
123  
124    for (i = 0; i < 3; i++ ) {
125      for (j = 0; j < 3; j++ ) {
126        vScale[i][j] = eta[i][j];
127 <      
127 >
128        if (i == j) {
129 <        vScale[i][j] += chi;          
130 <      }              
129 >        vScale[i][j] += chi;
130 >      }
131      }
132    }
109 –  
110 –  info->matVecMul3( vScale, vel, sc );
133   }
134  
135 + template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
136 +
137 +  matVecMul3( vScale, vel, sc );
138 + }
139 +
140   template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
141 <  int i,j;
141 >  int j;
142    double myVel[3];
116 –  double vScale[3][3];
143  
118 –  for (i = 0; i < 3; i++ ) {
119 –    for (j = 0; j < 3; j++ ) {
120 –      vScale[i][j] = eta[i][j];
121 –      
122 –      if (i == j) {
123 –        vScale[i][j] += chi;          
124 –      }              
125 –    }
126 –  }
127 –  
144    for (j = 0; j < 3; j++)
145      myVel[j] = oldVel[3*index + j];
146 <
147 <  info->matVecMul3( vScale, myVel, sc );
146 >  
147 >  matVecMul3( vScale, myVel, sc );
148   }
149  
150 < template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
150 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
151                                                 int index, double sc[3]){
152    int j;
153    double rj[3];
# Line 139 | Line 155 | template<typename T> void NPTf<T>::getPosScale(double
155    for(j=0; j<3; j++)
156      rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
157  
158 <  info->matVecMul3( eta, rj, sc );
158 >  matVecMul3( eta, rj, sc );
159   }
160  
161   template<typename T> void NPTf<T>::scaleSimBox( void ){
# Line 149 | Line 165 | template<typename T> void NPTf<T>::scaleSimBox( void )
165    double eta2ij;
166    double bigScale, smallScale, offDiagMax;
167    double hm[3][3], hmnew[3][3];
152 –  
168  
169  
170 +
171    // Scale the box after all the positions have been moved:
172 <  
172 >
173    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
174    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
175 <  
175 >
176    bigScale = 1.0;
177    smallScale = 1.0;
178    offDiagMax = 0.0;
179 <  
179 >
180    for(i=0; i<3; i++){
181      for(j=0; j<3; j++){
182 <      
182 >
183        // Calculate the matrix Product of the eta array (we only need
184        // the ij element right now):
185 <      
185 >
186        eta2ij = 0.0;
187        for(k=0; k<3; k++){
188          eta2ij += eta[i][k] * eta[k][j];
189        }
190 <      
190 >
191        scaleMat[i][j] = 0.0;
192        // identity matrix (see above):
193        if (i == j) scaleMat[i][j] = 1.0;
194        // Taylor expansion for the exponential truncated at second order:
195        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
196 +      
197  
198        if (i != j)
199 <        if (fabs(scaleMat[i][j]) > offDiagMax)
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 ((bigScale > 1.1) || (smallScale < 0.9)) {
206 >
207 >  if ((bigScale > 1.01) || (smallScale < 0.99)) {
208      sprintf( painCave.errMsg,
209 <             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
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 +             "      eta = [%lf\t%lf\t%lf]\n"
215 +             "            [%lf\t%lf\t%lf]\n"
216               "            [%lf\t%lf\t%lf]\n",
217               scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
218               scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
219 <             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
219 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2],
220 >             eta[0][0],eta[0][1],eta[0][2],
221 >             eta[1][0],eta[1][1],eta[1][2],
222 >             eta[2][0],eta[2][1],eta[2][2]);
223      painCave.isFatal = 1;
224      simError();
225 <  } else if (offDiagMax > 0.1) {
225 >  } else if (offDiagMax > 0.01) {
226      sprintf( painCave.errMsg,
227 <             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
227 >             "NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
228               " Check your tauBarostat, as it is probably too small!\n\n"
229               " scaleMat = [%lf\t%lf\t%lf]\n"
230               "            [%lf\t%lf\t%lf]\n"
231 +             "            [%lf\t%lf\t%lf]\n"
232 +             "      eta = [%lf\t%lf\t%lf]\n"
233 +             "            [%lf\t%lf\t%lf]\n"
234               "            [%lf\t%lf\t%lf]\n",
235               scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
236               scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
237 <             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
237 >             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2],
238 >             eta[0][0],eta[0][1],eta[0][2],
239 >             eta[1][0],eta[1][1],eta[1][2],
240 >             eta[2][0],eta[2][1],eta[2][2]);
241      painCave.isFatal = 1;
242      simError();
243    } else {
244      info->getBoxM(hm);
245 <    info->matMul3(hm, scaleMat, hmnew);
245 >    matMul3(hm, scaleMat, hmnew);
246      info->setBoxM(hmnew);
247    }
248   }
# Line 224 | Line 253 | template<typename T> bool NPTf<T>::etaConverged() {
253  
254    sumEta = 0;
255    for(i = 0; i < 3; i++)
256 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
257 <  
256 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
257 >
258    diffEta = sqrt( sumEta / 3.0 );
259 <  
259 >
260    return ( diffEta <= etaTolerance );
261   }
262  
263   template<typename T> double NPTf<T>::getConservedQuantity(void){
264 <  
264 >
265    double conservedQuantity;
266    double totalEnergy;
267    double thermostat_kinetic;
# Line 244 | Line 273 | template<typename T> double NPTf<T>::getConservedQuant
273  
274    totalEnergy = tStats->getTotalE();
275  
276 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
276 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
277      (2.0 * eConvert);
278  
279    thermostat_potential = fkBT* integralOfChidt / eConvert;
280  
281 <  info->transposeMat3(eta, a);
282 <  info->matMul3(a, eta, b);
283 <  trEta = info->matTrace3(b);
281 >  transposeMat3(eta, a);
282 >  matMul3(a, eta, b);
283 >  trEta = matTrace3(b);
284  
285 <  barostat_kinetic = NkBT * tb2 * trEta /
285 >  barostat_kinetic = NkBT * tb2 * trEta /
286      (2.0 * eConvert);
287 <  
288 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
287 >
288 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
289      eConvert;
290  
291    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
292      barostat_kinetic + barostat_potential;
264 –  
265 – //   cout.width(8);
266 – //   cout.precision(8);
293  
294 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270 < //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
294 >  return conservedQuantity;
295  
272 –  return conservedQuantity;
273 –  
296   }
297 +
298 + template<typename T> string NPTf<T>::getAdditionalParameters(void){
299 +  string parameters;
300 +  const int BUFFERSIZE = 2000; // size of the read buffer
301 +  char buffer[BUFFERSIZE];
302 +
303 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
304 +  parameters += buffer;
305 +
306 +  for(int i = 0; i < 3; i++){
307 +    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
308 +    parameters += buffer;
309 +  }
310 +
311 +  return parameters;
312 +
313 + }

Diff Legend

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