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 855 by mmeineke, Thu Nov 6 22:01:37 2003 UTC

# Line 7 | Line 7
7   #include "Thermo.hpp"
8   #include "ReadWrite.hpp"
9   #include "Integrator.hpp"
10 < #include "simError.h"
10 > #include "simError.h"
11  
12   #ifdef IS_MPI
13   #include "mpiSimulation.hpp"
# Line 17 | Line 17
17   // modification of the Hoover algorithm:
18   //
19   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
20 < //       Molec. Phys., 78, 533.
20 > //       Molec. Phys., 78, 533.
21   //
22   //           and
23 < //
23 > //
24   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
25  
26   template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
27    T( theInfo, the_ff )
28   {
29 <  
29 >  GenericData* data;
30 >  DoubleArrayData * etaValue;
31 >  vector<double> etaArray;
32    int i,j;
33 <  
33 >
34    for(i = 0; i < 3; i++){
35      for (j = 0; j < 3; j++){
36 <      
36 >
37        eta[i][j] = 0.0;
38        oldEta[i][j] = 0.0;
39      }
40    }
41 +
42 +
43 +  if( theInfo->useInitXSstate ){
44 +    // retrieve eta array from simInfo if it exists
45 +    data = info->getProperty(ETAVALUE_ID);
46 +    if(data){
47 +      etaValue = dynamic_cast<DoubleArrayData*>(data);
48 +      
49 +      if(etaValue){
50 +        etaArray = etaValue->getData();
51 +        
52 +        for(i = 0; i < 3; i++){
53 +          for (j = 0; j < 3; j++){
54 +            eta[i][j] = etaArray[3*i+j];
55 +            oldEta[i][j] = eta[i][j];
56 +          }
57 +        }
58 +      }
59 +    }
60 +  }
61 +
62   }
63  
64   template<typename T> NPTf<T>::~NPTf() {
# Line 44 | Line 67 | template<typename T> void NPTf<T>::resetIntegrator() {
67   }
68  
69   template<typename T> void NPTf<T>::resetIntegrator() {
70 <  
70 >
71    int i, j;
72 <  
72 >
73    for(i = 0; i < 3; i++)
74      for (j = 0; j < 3; j++)
75        eta[i][j] = 0.0;
76 <  
76 >
77    T::resetIntegrator();
78   }
79  
80   template<typename T> void NPTf<T>::evolveEtaA() {
81 <  
81 >
82    int i, j;
83 <  
83 >
84    for(i = 0; i < 3; i ++){
85      for(j = 0; j < 3; j++){
86        if( i == j)
87 <        eta[i][j] += dt2 *  instaVol *
87 >        eta[i][j] += dt2 *  instaVol *
88            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
89        else
90          eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
91      }
92    }
93 <  
93 >
94    for(i = 0; i < 3; i++)
95      for (j = 0; j < 3; j++)
96        oldEta[i][j] = eta[i][j];
97   }
98  
99   template<typename T> void NPTf<T>::evolveEtaB() {
100 <  
100 >
101    int i,j;
102  
103    for(i = 0; i < 3; i++)
# Line 84 | Line 107 | template<typename T> void NPTf<T>::evolveEtaB() {
107    for(i = 0; i < 3; i ++){
108      for(j = 0; j < 3; j++){
109        if( i == j) {
110 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
110 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
111            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
112        } else {
113          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
# Line 100 | Line 123 | template<typename T> void NPTf<T>::getVelScaleA(double
123    for (i = 0; i < 3; i++ ) {
124      for (j = 0; j < 3; j++ ) {
125        vScale[i][j] = eta[i][j];
126 <      
126 >
127        if (i == j) {
128 <        vScale[i][j] += chi;          
129 <      }              
128 >        vScale[i][j] += chi;
129 >      }
130      }
131    }
132 <  
132 >
133    info->matVecMul3( vScale, vel, sc );
134   }
135  
# Line 118 | Line 141 | template<typename T> void NPTf<T>::getVelScaleB(double
141    for (i = 0; i < 3; i++ ) {
142      for (j = 0; j < 3; j++ ) {
143        vScale[i][j] = eta[i][j];
144 <      
144 >
145        if (i == j) {
146 <        vScale[i][j] += chi;          
147 <      }              
146 >        vScale[i][j] += chi;
147 >      }
148      }
149    }
150 <  
150 >
151    for (j = 0; j < 3; j++)
152      myVel[j] = oldVel[3*index + j];
153  
154    info->matVecMul3( vScale, myVel, sc );
155   }
156  
157 < template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
157 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
158                                                 int index, double sc[3]){
159    int j;
160    double rj[3];
# Line 149 | Line 172 | template<typename T> void NPTf<T>::scaleSimBox( void )
172    double eta2ij;
173    double bigScale, smallScale, offDiagMax;
174    double hm[3][3], hmnew[3][3];
152 –  
175  
176  
177 +
178    // Scale the box after all the positions have been moved:
179 <  
179 >
180    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
181    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
182 <  
182 >
183    bigScale = 1.0;
184    smallScale = 1.0;
185    offDiagMax = 0.0;
186 <  
186 >
187    for(i=0; i<3; i++){
188      for(j=0; j<3; j++){
189 <      
189 >
190        // Calculate the matrix Product of the eta array (we only need
191        // the ij element right now):
192 <      
192 >
193        eta2ij = 0.0;
194        for(k=0; k<3; k++){
195          eta2ij += eta[i][k] * eta[k][j];
196        }
197 <      
197 >
198        scaleMat[i][j] = 0.0;
199        // identity matrix (see above):
200        if (i == j) scaleMat[i][j] = 1.0;
# Line 179 | Line 202 | template<typename T> void NPTf<T>::scaleSimBox( void )
202        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
203  
204        if (i != j)
205 <        if (fabs(scaleMat[i][j]) > offDiagMax)
205 >        if (fabs(scaleMat[i][j]) > offDiagMax)
206            offDiagMax = fabs(scaleMat[i][j]);
207      }
208  
209      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
210      if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
211    }
212 <  
213 <  if ((bigScale > 1.1) || (smallScale < 0.9)) {
212 >
213 >  if ((bigScale > 1.01) || (smallScale < 0.99)) {
214      sprintf( painCave.errMsg,
215 <             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
215 >             "NPTf error: Attempting a Box scaling of more than 1 percent.\n"
216               " Check your tauBarostat, as it is probably too small!\n\n"
217               " scaleMat = [%lf\t%lf\t%lf]\n"
218               "            [%lf\t%lf\t%lf]\n"
# Line 199 | Line 222 | template<typename T> void NPTf<T>::scaleSimBox( void )
222               scaleMat[2][0],scaleMat[2][1],scaleMat[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"
# Line 224 | Line 247 | template<typename T> bool NPTf<T>::etaConverged() {
247  
248    sumEta = 0;
249    for(i = 0; i < 3; i++)
250 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
251 <  
250 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
251 >
252    diffEta = sqrt( sumEta / 3.0 );
253 <  
253 >
254    return ( diffEta <= etaTolerance );
255   }
256  
257   template<typename T> double NPTf<T>::getConservedQuantity(void){
258 <  
258 >
259    double conservedQuantity;
260    double totalEnergy;
261    double thermostat_kinetic;
# Line 244 | Line 267 | template<typename T> double NPTf<T>::getConservedQuant
267  
268    totalEnergy = tStats->getTotalE();
269  
270 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
270 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
271      (2.0 * eConvert);
272  
273    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 253 | Line 276 | template<typename T> double NPTf<T>::getConservedQuant
276    info->matMul3(a, eta, b);
277    trEta = info->matTrace3(b);
278  
279 <  barostat_kinetic = NkBT * tb2 * trEta /
279 >  barostat_kinetic = NkBT * tb2 * trEta /
280      (2.0 * eConvert);
281 <  
282 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
281 >
282 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
283      eConvert;
284  
285    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
286      barostat_kinetic + barostat_potential;
264 –  
265 – //   cout.width(8);
266 – //   cout.precision(8);
287  
288 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270 < //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
288 >  return conservedQuantity;
289  
272 –  return conservedQuantity;
273 –  
290   }
291 +
292 + template<typename T> string NPTf<T>::getAdditionalParameters(void){
293 +  string parameters;
294 +  const int BUFFERSIZE = 2000; // size of the read buffer
295 +  char buffer[BUFFERSIZE];
296 +
297 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
298 +  parameters += buffer;
299 +
300 +  for(int i = 0; i < 3; i++){
301 +    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
302 +    parameters += buffer;
303 +  }
304 +
305 +  return parameters;
306 +
307 + }

Diff Legend

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