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 836 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> NPTxyz<T>::NPTxyz ( 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   template<typename T> NPTxyz<T>::~NPTxyz() {
# Line 44 | Line 64 | template<typename T> void NPTxyz<T>::resetIntegrator()
64   }
65  
66   template<typename T> void NPTxyz<T>::resetIntegrator() {
67 <  
67 >
68    int i, j;
69 <  
69 >
70    for(i = 0; i < 3; i++)
71      for (j = 0; j < 3; j++)
72        eta[i][j] = 0.0;
73 <  
73 >
74    T::resetIntegrator();
75   }
76  
77   template<typename T> void NPTxyz<T>::evolveEtaA() {
78 <  
78 >
79    int i, j;
80 <  
80 >
81    for(i = 0; i < 3; i ++){
82      for(j = 0; j < 3; j++){
83        if( i == j)
84 <        eta[i][j] += dt2 *  instaVol *
84 >        eta[i][j] += dt2 *  instaVol *
85            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
86        else
87          eta[i][j] = 0.0;
88      }
89    }
90 <  
90 >
91    for(i = 0; i < 3; i++)
92      for (j = 0; j < 3; j++)
93        oldEta[i][j] = eta[i][j];
94   }
95  
96   template<typename T> void NPTxyz<T>::evolveEtaB() {
97 <  
97 >
98    int i,j;
99  
100    for(i = 0; i < 3; i++)
# Line 84 | Line 104 | template<typename T> void NPTxyz<T>::evolveEtaB() {
104    for(i = 0; i < 3; i ++){
105      for(j = 0; j < 3; j++){
106        if( i == j) {
107 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
107 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
108            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
109        } else {
110          eta[i][j] = 0.0;
# Line 100 | Line 120 | template<typename T> void NPTxyz<T>::getVelScaleA(doub
120    for (i = 0; i < 3; i++ ) {
121      for (j = 0; j < 3; j++ ) {
122        vScale[i][j] = eta[i][j];
123 <      
123 >
124        if (i == j) {
125 <        vScale[i][j] += chi;          
126 <      }              
125 >        vScale[i][j] += chi;
126 >      }
127      }
128    }
129 <  
129 >
130    info->matVecMul3( vScale, vel, sc );
131   }
132  
# Line 118 | Line 138 | template<typename T> void NPTxyz<T>::getVelScaleB(doub
138    for (i = 0; i < 3; i++ ) {
139      for (j = 0; j < 3; j++ ) {
140        vScale[i][j] = eta[i][j];
141 <      
141 >
142        if (i == j) {
143 <        vScale[i][j] += chi;          
144 <      }              
143 >        vScale[i][j] += chi;
144 >      }
145      }
146    }
147 <  
147 >
148    for (j = 0; j < 3; j++)
149      myVel[j] = oldVel[3*index + j];
150  
151    info->matVecMul3( vScale, myVel, sc );
152   }
153  
154 < template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
154 > template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
155                                                 int index, double sc[3]){
156    int j;
157    double rj[3];
# Line 149 | Line 169 | template<typename T> void NPTxyz<T>::scaleSimBox( void
169    double eta2ij, scaleFactor;
170    double bigScale, smallScale, offDiagMax;
171    double hm[3][3], hmnew[3][3];
152 –  
172  
173  
174 +
175    // Scale the box after all the positions have been moved:
176 <  
176 >
177    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
178    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
179 <  
179 >
180    bigScale = 1.0;
181    smallScale = 1.0;
182    offDiagMax = 0.0;
183 <  
183 >
184    for(i=0; i<3; i++){
185      for(j=0; j<3; j++){
186        scaleMat[i][j] = 0.0;
187        if(i==j) scaleMat[i][j] = 1.0;
188      }
189    }
190 <  
190 >
191    for(i=0;i<3;i++){
192 <
192 >
193      // calculate the scaleFactors
194 <      
194 >
195      scaleFactor = exp(dt*eta[i][i]);
196 <    
196 >
197      scaleMat[i][i] = scaleFactor;
198  
199      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 182 | Line 202 | template<typename T> void NPTxyz<T>::scaleSimBox( void
202  
203   //   for(i=0; i<3; i++){
204   //     for(j=0; j<3; j++){
205 <      
205 >
206   //       // Calculate the matrix Product of the eta array (we only need
207   //       // the ij element right now):
208 <      
208 >
209   //       eta2ij = 0.0;
210   //       for(k=0; k<3; k++){
211   //         eta2ij += eta[i][k] * eta[k][j];
212   //       }
213 <      
213 >
214   //       scaleMat[i][j] = 0.0;
215   //       // identity matrix (see above):
216   //       if (i == j) scaleMat[i][j] = 1.0;
# Line 198 | Line 218 | template<typename T> void NPTxyz<T>::scaleSimBox( void
218   //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
219  
220   //       if (i != j)
221 < //         if (fabs(scaleMat[i][j]) > offDiagMax)
221 > //         if (fabs(scaleMat[i][j]) > offDiagMax)
222   //           offDiagMax = fabs(scaleMat[i][j]);
223   //     }
224  
225   //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
226   //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
227   //   }
228 <  
228 >
229    if ((bigScale > 1.1) || (smallScale < 0.9)) {
230      sprintf( painCave.errMsg,
231               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
# Line 231 | Line 251 | template<typename T> bool NPTxyz<T>::etaConverged() {
251  
252    sumEta = 0;
253    for(i = 0; i < 3; i++)
254 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
255 <  
254 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
255 >
256    diffEta = sqrt( sumEta / 3.0 );
257 <  
257 >
258    return ( diffEta <= etaTolerance );
259   }
260  
261   template<typename T> double NPTxyz<T>::getConservedQuantity(void){
262 <  
262 >
263    double conservedQuantity;
264    double totalEnergy;
265    double thermostat_kinetic;
# Line 251 | Line 271 | template<typename T> double NPTxyz<T>::getConservedQua
271  
272    totalEnergy = tStats->getTotalE();
273  
274 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
274 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
275      (2.0 * eConvert);
276  
277    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 260 | Line 280 | template<typename T> double NPTxyz<T>::getConservedQua
280    info->matMul3(a, eta, b);
281    trEta = info->matTrace3(b);
282  
283 <  barostat_kinetic = NkBT * tb2 * trEta /
283 >  barostat_kinetic = NkBT * tb2 * trEta /
284      (2.0 * eConvert);
285 <  
286 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
285 >
286 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
287      eConvert;
288  
289    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
290      barostat_kinetic + barostat_potential;
291 <  
291 >
292   //   cout.width(8);
293   //   cout.precision(8);
294  
295 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
296 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
295 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
296 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
297   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
298  
299 <  return conservedQuantity;
300 <  
299 >  return conservedQuantity;
300 >
301   }
302 +
303 + template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
304 +  string parameters;
305 +  const int BUFFERSIZE = 2000; // size of the read buffer
306 +  char buffer[BUFFERSIZE];
307 +
308 +  sprintf(buffer,"\t%f\t%f;", chi, integralOfChidt);
309 +  parameters += buffer;
310 +
311 +  for(int i = 0; i < 3; i++){
312 +    sprintf(buffer,"\t%g\t%g\t%g;", eta[3*i], eta[3*i+1], eta[3*i+2]);
313 +    parameters += buffer;
314 +  }
315 +
316 +  return parameters;
317 +
318 + }

Diff Legend

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