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 812 by mmeineke, Wed Oct 22 21:17:32 2003 UTC vs.
Revision 857 by mmeineke, Fri Nov 7 17:09:48 2003 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <math.h>
2   #include "Atom.hpp"
3   #include "SRI.hpp"
4   #include "AbstractClasses.hpp"
# 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 +
43 +  if( theInfo->useInitXSstate ){
44 +
45 +    // retrieve eta array from simInfo if it exists
46 +    data = info->getProperty(ETAVALUE_ID);
47 +    if(data){
48 +      etaValue = dynamic_cast<DoubleArrayData*>(data);
49 +      
50 +      if(etaValue){
51 +        etaArray = etaValue->getData();
52 +        
53 +        for(i = 0; i < 3; i++){
54 +          for (j = 0; j < 3; j++){
55 +            eta[i][j] = etaArray[3*i+j];
56 +            oldEta[i][j] = eta[i][j];
57 +          }
58 +        }
59 +      }
60 +    }
61 +  }
62   }
63  
64   template<typename T> NPTxyz<T>::~NPTxyz() {
# Line 44 | Line 67 | template<typename T> void NPTxyz<T>::resetIntegrator()
67   }
68  
69   template<typename T> void NPTxyz<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 NPTxyz<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] = 0.0;
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 NPTxyz<T>::evolveEtaB() {
100 <  
100 >
101    int i,j;
102  
103    for(i = 0; i < 3; i++)
# Line 84 | Line 107 | template<typename T> void NPTxyz<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] = 0.0;
# Line 93 | Line 116 | template<typename T> void NPTxyz<T>::getVelScaleA(doub
116    }
117   }
118  
119 < template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
119 > template<typename T> void NPTxyz<T>::calcVelScale(void) {
120    int i,j;
98 –  double vScale[3][3];
121  
122    for (i = 0; i < 3; i++ ) {
123      for (j = 0; j < 3; j++ ) {
124        vScale[i][j] = eta[i][j];
125 <      
125 >
126        if (i == j) {
127 <        vScale[i][j] += chi;          
128 <      }              
127 >        vScale[i][j] += chi;
128 >      }
129      }
130    }
131 <  
131 > }
132 >
133 > template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
134    info->matVecMul3( vScale, vel, sc );
135   }
136  
137   template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
138 <  int i,j;
138 >  int j;
139    double myVel[3];
116 –  double vScale[3][3];
140  
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 –  
141    for (j = 0; j < 3; j++)
142      myVel[j] = oldVel[3*index + j];
143  
144    info->matVecMul3( vScale, myVel, sc );
145   }
146  
147 < template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
147 > template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
148                                                 int index, double sc[3]){
149    int j;
150    double rj[3];
# Line 149 | Line 162 | template<typename T> void NPTxyz<T>::scaleSimBox( void
162    double eta2ij, scaleFactor;
163    double bigScale, smallScale, offDiagMax;
164    double hm[3][3], hmnew[3][3];
152 –  
165  
166  
167 +
168    // Scale the box after all the positions have been moved:
169 <  
169 >
170    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
171    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
172 <  
172 >
173    bigScale = 1.0;
174    smallScale = 1.0;
175    offDiagMax = 0.0;
176 <  
176 >
177    for(i=0; i<3; i++){
178      for(j=0; j<3; j++){
179        scaleMat[i][j] = 0.0;
180        if(i==j) scaleMat[i][j] = 1.0;
181      }
182    }
183 <  
183 >
184    for(i=0;i<3;i++){
185 <
185 >
186      // calculate the scaleFactors
187 <      
187 >
188      scaleFactor = exp(dt*eta[i][i]);
189 <    
189 >
190      scaleMat[i][i] = scaleFactor;
191  
192      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 182 | Line 195 | template<typename T> void NPTxyz<T>::scaleSimBox( void
195  
196   //   for(i=0; i<3; i++){
197   //     for(j=0; j<3; j++){
198 <      
198 >
199   //       // Calculate the matrix Product of the eta array (we only need
200   //       // the ij element right now):
201 <      
201 >
202   //       eta2ij = 0.0;
203   //       for(k=0; k<3; k++){
204   //         eta2ij += eta[i][k] * eta[k][j];
205   //       }
206 <      
206 >
207   //       scaleMat[i][j] = 0.0;
208   //       // identity matrix (see above):
209   //       if (i == j) scaleMat[i][j] = 1.0;
# Line 198 | Line 211 | template<typename T> void NPTxyz<T>::scaleSimBox( void
211   //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
212  
213   //       if (i != j)
214 < //         if (fabs(scaleMat[i][j]) > offDiagMax)
214 > //         if (fabs(scaleMat[i][j]) > offDiagMax)
215   //           offDiagMax = fabs(scaleMat[i][j]);
216   //     }
217  
218   //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
219   //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
220   //   }
221 <  
221 >
222    if ((bigScale > 1.1) || (smallScale < 0.9)) {
223      sprintf( painCave.errMsg,
224               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
# Line 231 | Line 244 | template<typename T> bool NPTxyz<T>::etaConverged() {
244  
245    sumEta = 0;
246    for(i = 0; i < 3; i++)
247 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
248 <  
247 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
248 >
249    diffEta = sqrt( sumEta / 3.0 );
250 <  
250 >
251    return ( diffEta <= etaTolerance );
252   }
253  
254   template<typename T> double NPTxyz<T>::getConservedQuantity(void){
255 <  
255 >
256    double conservedQuantity;
257    double totalEnergy;
258    double thermostat_kinetic;
# Line 251 | Line 264 | template<typename T> double NPTxyz<T>::getConservedQua
264  
265    totalEnergy = tStats->getTotalE();
266  
267 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
267 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
268      (2.0 * eConvert);
269  
270    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 260 | Line 273 | template<typename T> double NPTxyz<T>::getConservedQua
273    info->matMul3(a, eta, b);
274    trEta = info->matTrace3(b);
275  
276 <  barostat_kinetic = NkBT * tb2 * trEta /
276 >  barostat_kinetic = NkBT * tb2 * trEta /
277      (2.0 * eConvert);
278 <  
279 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
278 >
279 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
280      eConvert;
281  
282    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
283      barostat_kinetic + barostat_potential;
284 <  
284 >
285   //   cout.width(8);
286   //   cout.precision(8);
287  
288 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
289 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
288 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
289 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
290   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
291  
292 <  return conservedQuantity;
293 <  
292 >  return conservedQuantity;
293 >
294   }
295 +
296 + template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
297 +  string parameters;
298 +  const int BUFFERSIZE = 2000; // size of the read buffer
299 +  char buffer[BUFFERSIZE];
300 +
301 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
302 +  parameters += buffer;
303 +
304 +  for(int i = 0; i < 3; i++){
305 +    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
306 +    parameters += buffer;
307 +  }
308 +
309 +  return parameters;
310 +
311 + }

Diff Legend

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