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

Comparing trunk/OOPSE/libmdtools/NPTxyz.cpp (file contents):
Revision 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 1 | Line 1
1   #include <math.h>
2 + #include "MatVec3.h"
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 17 | Line 18
18   // modification of the Hoover algorithm:
19   //
20   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
21 < //       Molec. Phys., 78, 533.
21 > //       Molec. Phys., 78, 533.
22   //
23   //           and
24 < //
24 > //
25   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26  
27   template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff):
28    T( theInfo, the_ff )
29   {
30 <  
30 >  GenericData* data;
31 >  DoubleArrayData * etaValue;
32 >  vector<double> etaArray;
33    int i,j;
34 <  
34 >
35    for(i = 0; i < 3; i++){
36      for (j = 0; j < 3; j++){
37 <      
37 >
38        eta[i][j] = 0.0;
39        oldEta[i][j] = 0.0;
40      }
41    }
42 +
43 +
44 +  if( theInfo->useInitXSstate ){
45 +
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   template<typename T> NPTxyz<T>::~NPTxyz() {
# Line 44 | Line 68 | template<typename T> void NPTxyz<T>::resetIntegrator()
68   }
69  
70   template<typename T> void NPTxyz<T>::resetIntegrator() {
71 <  
71 >
72    int i, j;
73 <  
73 >
74    for(i = 0; i < 3; i++)
75      for (j = 0; j < 3; j++)
76        eta[i][j] = 0.0;
77 <  
77 >
78    T::resetIntegrator();
79   }
80  
81   template<typename T> void NPTxyz<T>::evolveEtaA() {
82 <  
82 >
83    int i, j;
84 <  
84 >
85    for(i = 0; i < 3; i ++){
86      for(j = 0; j < 3; j++){
87        if( i == j)
88 <        eta[i][j] += dt2 *  instaVol *
88 >        eta[i][j] += dt2 *  instaVol *
89            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
90        else
91          eta[i][j] = 0.0;
92      }
93    }
94 <  
94 >
95    for(i = 0; i < 3; i++)
96      for (j = 0; j < 3; j++)
97        oldEta[i][j] = eta[i][j];
98   }
99  
100   template<typename T> void NPTxyz<T>::evolveEtaB() {
101 <  
101 >
102    int i,j;
103  
104    for(i = 0; i < 3; i++)
# Line 84 | Line 108 | template<typename T> void NPTxyz<T>::evolveEtaB() {
108    for(i = 0; i < 3; i ++){
109      for(j = 0; j < 3; j++){
110        if( i == j) {
111 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
111 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
112            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
113        } else {
114          eta[i][j] = 0.0;
# Line 93 | Line 117 | template<typename T> void NPTxyz<T>::getVelScaleA(doub
117    }
118   }
119  
120 < template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
120 > template<typename T> void NPTxyz<T>::calcVelScale(void) {
121    int i,j;
98 –  double vScale[3][3];
122  
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    }
109 –  
110 –  info->matVecMul3( vScale, vel, sc );
132   }
133  
134 + template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
135 +  matVecMul3( vScale, vel, sc );
136 + }
137 +
138   template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
139 <  int i,j;
139 >  int j;
140    double myVel[3];
116 –  double vScale[3][3];
141  
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 –  
142    for (j = 0; j < 3; j++)
143      myVel[j] = oldVel[3*index + j];
144  
145 <  info->matVecMul3( vScale, myVel, sc );
145 >  matVecMul3( vScale, myVel, sc );
146   }
147  
148 < template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
148 > template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
149                                                 int index, double sc[3]){
150    int j;
151    double rj[3];
# Line 139 | Line 153 | template<typename T> void NPTxyz<T>::getPosScale(doubl
153    for(j=0; j<3; j++)
154      rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
155  
156 <  info->matVecMul3( eta, rj, sc );
156 >  matVecMul3( eta, rj, sc );
157   }
158  
159   template<typename T> void NPTxyz<T>::scaleSimBox( void ){
# Line 149 | Line 163 | template<typename T> void NPTxyz<T>::scaleSimBox( void
163    double eta2ij, scaleFactor;
164    double bigScale, smallScale, offDiagMax;
165    double hm[3][3], hmnew[3][3];
152 –  
166  
167  
168 +
169    // Scale the box after all the positions have been moved:
170 <  
170 >
171    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
172    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
173 <  
173 >
174    bigScale = 1.0;
175    smallScale = 1.0;
176    offDiagMax = 0.0;
177 <  
177 >
178    for(i=0; i<3; i++){
179      for(j=0; j<3; j++){
180        scaleMat[i][j] = 0.0;
181        if(i==j) scaleMat[i][j] = 1.0;
182      }
183    }
184 <  
184 >
185    for(i=0;i<3;i++){
186 <
186 >
187      // calculate the scaleFactors
188 <      
188 >
189      scaleFactor = exp(dt*eta[i][i]);
190 <    
190 >
191      scaleMat[i][i] = scaleFactor;
192  
193      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 182 | Line 196 | template<typename T> void NPTxyz<T>::scaleSimBox( void
196  
197   //   for(i=0; i<3; i++){
198   //     for(j=0; j<3; j++){
199 <      
199 >
200   //       // Calculate the matrix Product of the eta array (we only need
201   //       // the ij element right now):
202 <      
202 >
203   //       eta2ij = 0.0;
204   //       for(k=0; k<3; k++){
205   //         eta2ij += eta[i][k] * eta[k][j];
206   //       }
207 <      
207 >
208   //       scaleMat[i][j] = 0.0;
209   //       // identity matrix (see above):
210   //       if (i == j) scaleMat[i][j] = 1.0;
# Line 198 | Line 212 | template<typename T> void NPTxyz<T>::scaleSimBox( void
212   //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
213  
214   //       if (i != j)
215 < //         if (fabs(scaleMat[i][j]) > offDiagMax)
215 > //         if (fabs(scaleMat[i][j]) > offDiagMax)
216   //           offDiagMax = fabs(scaleMat[i][j]);
217   //     }
218  
219   //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
220   //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
221   //   }
222 <  
222 >
223    if ((bigScale > 1.1) || (smallScale < 0.9)) {
224      sprintf( painCave.errMsg,
225               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
# Line 220 | Line 234 | template<typename T> void NPTxyz<T>::scaleSimBox( void
234      simError();
235    } else {
236      info->getBoxM(hm);
237 <    info->matMul3(hm, scaleMat, hmnew);
237 >    matMul3(hm, scaleMat, hmnew);
238      info->setBoxM(hmnew);
239    }
240   }
# Line 231 | Line 245 | template<typename T> bool NPTxyz<T>::etaConverged() {
245  
246    sumEta = 0;
247    for(i = 0; i < 3; i++)
248 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
249 <  
248 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
249 >
250    diffEta = sqrt( sumEta / 3.0 );
251 <  
251 >
252    return ( diffEta <= etaTolerance );
253   }
254  
255   template<typename T> double NPTxyz<T>::getConservedQuantity(void){
256 <  
256 >
257    double conservedQuantity;
258    double totalEnergy;
259    double thermostat_kinetic;
# Line 251 | Line 265 | template<typename T> double NPTxyz<T>::getConservedQua
265  
266    totalEnergy = tStats->getTotalE();
267  
268 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
268 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
269      (2.0 * eConvert);
270  
271    thermostat_potential = fkBT* integralOfChidt / eConvert;
272  
273 <  info->transposeMat3(eta, a);
274 <  info->matMul3(a, eta, b);
275 <  trEta = info->matTrace3(b);
273 >  transposeMat3(eta, a);
274 >  matMul3(a, eta, b);
275 >  trEta = matTrace3(b);
276  
277 <  barostat_kinetic = NkBT * tb2 * trEta /
277 >  barostat_kinetic = NkBT * tb2 * trEta /
278      (2.0 * eConvert);
279 <  
280 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
279 >
280 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
281      eConvert;
282  
283    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
284      barostat_kinetic + barostat_potential;
285 <  
285 >
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 <<
289 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
290 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
291   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
292  
293 <  return conservedQuantity;
294 <  
293 >  return conservedQuantity;
294 >
295   }
296 +
297 + template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
298 +  string parameters;
299 +  const int BUFFERSIZE = 2000; // size of the read buffer
300 +  char buffer[BUFFERSIZE];
301 +
302 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
303 +  parameters += buffer;
304 +
305 +  for(int i = 0; i < 3; i++){
306 +    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
307 +    parameters += buffer;
308 +  }
309 +
310 +  return parameters;
311 +
312 + }

Diff Legend

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