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 847 by mmeineke, Fri Oct 31 18:28:52 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 39 | Line 41 | template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
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 <
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 <        }
56 <
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  
# Line 88 | Line 92 | template<typename T> void NPTf<T>::evolveEtaA() {
92          eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
93      }
94    }
95 <
95 >  
96    for(i = 0; i < 3; i++)
97      for (j = 0; j < 3; j++)
98        oldEta[i][j] = eta[i][j];
# Line 114 | 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;
119 –  double vScale[3][3];
123  
124    for (i = 0; i < 3; i++ ) {
125      for (j = 0; j < 3; j++ ) {
# Line 127 | Line 130 | template<typename T> void NPTf<T>::getVelScaleA(double
130        }
131      }
132    }
133 + }
134  
135 <  info->matVecMul3( vScale, vel, sc );
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];
137 –  double vScale[3][3];
143  
139 –  for (i = 0; i < 3; i++ ) {
140 –    for (j = 0; j < 3; j++ ) {
141 –      vScale[i][j] = eta[i][j];
142 –
143 –      if (i == j) {
144 –        vScale[i][j] += chi;
145 –      }
146 –    }
147 –  }
148 –
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],
# Line 160 | 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 198 | Line 193 | template<typename T> void NPTf<T>::scaleSimBox( void )
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)
# Line 208 | Line 204 | template<typename T> void NPTf<T>::scaleSimBox( void )
204      if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
205    }
206  
207 <  if ((bigScale > 1.1) || (smallScale < 0.9)) {
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 270 | Line 278 | template<typename T> double NPTf<T>::getConservedQuant
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 /
286      (2.0 * eConvert);
# Line 283 | Line 291 | template<typename T> double NPTf<T>::getConservedQuant
291    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
292      barostat_kinetic + barostat_potential;
293  
286 – //   cout.width(8);
287 – //   cout.precision(8);
288 –
289 – //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
290 – //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
291 – //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
292 –
294    return conservedQuantity;
295  
296   }
# Line 299 | Line 300 | template<typename T> string NPTf<T>::getAdditionalPara
300    const int BUFFERSIZE = 2000; // size of the read buffer
301    char buffer[BUFFERSIZE];
302  
303 <  sprintf(buffer,"\t%lf\t%lf;", chi, integralOfChidt);
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%lf\t%lf\t%lf;", eta[3*i], eta[3*i+1], eta[3*i+2]);
307 >    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
308      parameters += buffer;
309    }
310  

Diff Legend

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