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 837 by tim, Wed Oct 29 00:19:10 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 +    // retrieve eta array from simInfo if it exists
43 +    data = info->getProperty(ETAVALUE_ID);
44 +    if(data){
45 +      etaValue = dynamic_cast<DoubleArrayData*>(data);
46 +
47 +      if(etaValue){
48 +        etaArray = etaValue->getData();
49 +
50 +        for(i = 0; i < 3; i++){
51 +          for (j = 0; j < 3; j++){
52 +            eta[i][j] = etaArray[3*i+j];
53 +            oldEta[i][j] = eta[i][j];
54 +          }
55 +        }
56 +
57 +      }
58 +    }
59 +
60   }
61  
62   template<typename T> NPTf<T>::~NPTf() {
# Line 44 | Line 65 | template<typename T> void NPTf<T>::resetIntegrator() {
65   }
66  
67   template<typename T> void NPTf<T>::resetIntegrator() {
68 <  
68 >
69    int i, j;
70 <  
70 >
71    for(i = 0; i < 3; i++)
72      for (j = 0; j < 3; j++)
73        eta[i][j] = 0.0;
74 <  
74 >
75    T::resetIntegrator();
76   }
77  
78   template<typename T> void NPTf<T>::evolveEtaA() {
79 <  
79 >
80    int i, j;
81 <  
81 >
82    for(i = 0; i < 3; i ++){
83      for(j = 0; j < 3; j++){
84        if( i == j)
85 <        eta[i][j] += dt2 *  instaVol *
85 >        eta[i][j] += dt2 *  instaVol *
86            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
87        else
88          eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
89      }
90    }
91 <  
91 >
92    for(i = 0; i < 3; i++)
93      for (j = 0; j < 3; j++)
94        oldEta[i][j] = eta[i][j];
95   }
96  
97   template<typename T> void NPTf<T>::evolveEtaB() {
98 <  
98 >
99    int i,j;
100  
101    for(i = 0; i < 3; i++)
# Line 84 | Line 105 | template<typename T> void NPTf<T>::evolveEtaB() {
105    for(i = 0; i < 3; i ++){
106      for(j = 0; j < 3; j++){
107        if( i == j) {
108 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
108 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
109            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
110        } else {
111          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
# Line 100 | Line 121 | template<typename T> void NPTf<T>::getVelScaleA(double
121    for (i = 0; i < 3; i++ ) {
122      for (j = 0; j < 3; j++ ) {
123        vScale[i][j] = eta[i][j];
124 <      
124 >
125        if (i == j) {
126 <        vScale[i][j] += chi;          
127 <      }              
126 >        vScale[i][j] += chi;
127 >      }
128      }
129    }
130 <  
130 >
131    info->matVecMul3( vScale, vel, sc );
132   }
133  
# Line 118 | Line 139 | template<typename T> void NPTf<T>::getVelScaleB(double
139    for (i = 0; i < 3; i++ ) {
140      for (j = 0; j < 3; j++ ) {
141        vScale[i][j] = eta[i][j];
142 <      
142 >
143        if (i == j) {
144 <        vScale[i][j] += chi;          
145 <      }              
144 >        vScale[i][j] += chi;
145 >      }
146      }
147    }
148 <  
148 >
149    for (j = 0; j < 3; j++)
150      myVel[j] = oldVel[3*index + j];
151  
152    info->matVecMul3( vScale, myVel, sc );
153   }
154  
155 < template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
155 > template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
156                                                 int index, double sc[3]){
157    int j;
158    double rj[3];
# Line 149 | Line 170 | template<typename T> void NPTf<T>::scaleSimBox( void )
170    double eta2ij;
171    double bigScale, smallScale, offDiagMax;
172    double hm[3][3], hmnew[3][3];
152 –  
173  
174  
175 +
176    // Scale the box after all the positions have been moved:
177 <  
177 >
178    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
179    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
180 <  
180 >
181    bigScale = 1.0;
182    smallScale = 1.0;
183    offDiagMax = 0.0;
184 <  
184 >
185    for(i=0; i<3; i++){
186      for(j=0; j<3; j++){
187 <      
187 >
188        // Calculate the matrix Product of the eta array (we only need
189        // the ij element right now):
190 <      
190 >
191        eta2ij = 0.0;
192        for(k=0; k<3; k++){
193          eta2ij += eta[i][k] * eta[k][j];
194        }
195 <      
195 >
196        scaleMat[i][j] = 0.0;
197        // identity matrix (see above):
198        if (i == j) scaleMat[i][j] = 1.0;
# Line 179 | Line 200 | template<typename T> void NPTf<T>::scaleSimBox( void )
200        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
201  
202        if (i != j)
203 <        if (fabs(scaleMat[i][j]) > offDiagMax)
203 >        if (fabs(scaleMat[i][j]) > offDiagMax)
204            offDiagMax = fabs(scaleMat[i][j]);
205      }
206  
207      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
208      if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
209    }
210 <  
210 >
211    if ((bigScale > 1.1) || (smallScale < 0.9)) {
212      sprintf( painCave.errMsg,
213               "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
# Line 224 | Line 245 | template<typename T> bool NPTf<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 NPTf<T>::getConservedQuantity(void){
256 <  
256 >
257    double conservedQuantity;
258    double totalEnergy;
259    double thermostat_kinetic;
# Line 244 | Line 265 | template<typename T> double NPTf<T>::getConservedQuant
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;
# Line 253 | Line 274 | template<typename T> double NPTf<T>::getConservedQuant
274    info->matMul3(a, eta, b);
275    trEta = info->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 NPTf<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%f\t%f;", chi, integralOfChidt);
303 +  parameters += buffer;
304 +
305 +  for(int i = 0; i < 3; i++){
306 +    sprintf(buffer,"\t%g\t%g\t%g;", eta[3*i], eta[3*i+1], eta[3*i+2]);
307 +    parameters += buffer;
308 +  }
309 +
310 +  return parameters;
311 +
312 + }

Diff Legend

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