ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/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.
branches/new-templateless/OOPSE/libmdtools/NPTxyz.cpp (file contents), Revision 851 by mmeineke, Wed Nov 5 19:18:17 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 )
26 > NPTxyz::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff):
27 >  NPT( theInfo, the_ff )
28   {
29 <  
29 >  GenericData* data;
30 >  double *etaArray;
31    int i,j;
32 <  
32 >
33    for(i = 0; i < 3; i++){
34      for (j = 0; j < 3; j++){
35 <      
35 >
36        eta[i][j] = 0.0;
37        oldEta[i][j] = 0.0;
38      }
39    }
40 +
41 +  // retrieve eta array from simInfo if it exists
42 +  data = info->getProperty(ETAVALUE_ID);
43 +  if(data != NULL){
44 +    
45 +    int test = data->getDarray(etaArray);
46 +    
47 +    if( test == 9 ){
48 +      
49 +      for(i = 0; i < 3; i++){
50 +        for (j = 0; j < 3; j++){
51 +          eta[i][j] = etaArray[3*i+j];
52 +          oldEta[i][j] = eta[i][j];
53 +        }
54 +      }    
55 +      delete[] etaArray;
56 +    }
57 +    else
58 +      std::cerr << "NPTxyz error: etaArray is not length 9 (actual = " << test
59 +                << ").\n"
60 +                << "            Simulation wil proceed with eta = 0;\n";
61 +  }
62   }
63  
64 < template<typename T> NPTxyz<T>::~NPTxyz() {
64 > NPTxyz::~NPTxyz() {
65  
66    // empty for now
67   }
68  
69 < template<typename T> void NPTxyz<T>::resetIntegrator() {
70 <  
69 > void NPTxyz::resetIntegrator() {
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 <  
77 <  T::resetIntegrator();
76 >
77 >  NPT::resetIntegrator();
78   }
79  
80 < template<typename T> void NPTxyz<T>::evolveEtaA() {
81 <  
80 > void NPTxyz::evolveEtaA() {
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 <  
99 > void NPTxyz::evolveEtaB() {
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>::evolveEtaB() {
116    }
117   }
118  
119 < template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
119 > void NPTxyz::getVelScaleA(double sc[3], double vel[3]) {
120    int i,j;
121    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    }
132 <  
132 >
133    info->matVecMul3( vScale, vel, sc );
134   }
135  
136 < template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
136 > void NPTxyz::getVelScaleB(double sc[3], int index ){
137    int i,j;
138    double myVel[3];
139    double vScale[3][3];
# Line 118 | Line 141 | template<typename T> void NPTxyz<T>::getVelScaleB(doub
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 NPTxyz<T>::getPosScale(double pos[3], double COM[3],
157 > void NPTxyz::getPosScale(double pos[3], double COM[3],
158                                                 int index, double sc[3]){
159    int j;
160    double rj[3];
# Line 142 | Line 165 | template<typename T> void NPTxyz<T>::getPosScale(doubl
165    info->matVecMul3( eta, rj, sc );
166   }
167  
168 < template<typename T> void NPTxyz<T>::scaleSimBox( void ){
168 > void NPTxyz::scaleSimBox( void ){
169  
170    int i,j,k;
171    double scaleMat[3][3];
172    double eta2ij, scaleFactor;
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        scaleMat[i][j] = 0.0;
190        if(i==j) scaleMat[i][j] = 1.0;
191      }
192    }
193 <  
193 >
194    for(i=0;i<3;i++){
195 <
195 >
196      // calculate the scaleFactors
197 <      
197 >
198      scaleFactor = exp(dt*eta[i][i]);
199 <    
199 >
200      scaleMat[i][i] = scaleFactor;
201  
202      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
# Line 182 | Line 205 | template<typename T> void NPTxyz<T>::scaleSimBox( void
205  
206   //   for(i=0; i<3; i++){
207   //     for(j=0; j<3; j++){
208 <      
208 >
209   //       // Calculate the matrix Product of the eta array (we only need
210   //       // the ij element right now):
211 <      
211 >
212   //       eta2ij = 0.0;
213   //       for(k=0; k<3; k++){
214   //         eta2ij += eta[i][k] * eta[k][j];
215   //       }
216 <      
216 >
217   //       scaleMat[i][j] = 0.0;
218   //       // identity matrix (see above):
219   //       if (i == j) scaleMat[i][j] = 1.0;
# Line 198 | Line 221 | template<typename T> void NPTxyz<T>::scaleSimBox( void
221   //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
222  
223   //       if (i != j)
224 < //         if (fabs(scaleMat[i][j]) > offDiagMax)
224 > //         if (fabs(scaleMat[i][j]) > offDiagMax)
225   //           offDiagMax = fabs(scaleMat[i][j]);
226   //     }
227  
228   //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
229   //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
230   //   }
231 <  
231 >
232    if ((bigScale > 1.1) || (smallScale < 0.9)) {
233      sprintf( painCave.errMsg,
234               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
# Line 225 | Line 248 | template<typename T> void NPTxyz<T>::scaleSimBox( void
248    }
249   }
250  
251 < template<typename T> bool NPTxyz<T>::etaConverged() {
251 > bool NPTxyz::etaConverged() {
252    int i;
253    double diffEta, sumEta;
254  
255    sumEta = 0;
256    for(i = 0; i < 3; i++)
257 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
258 <  
257 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
258 >
259    diffEta = sqrt( sumEta / 3.0 );
260 <  
260 >
261    return ( diffEta <= etaTolerance );
262   }
263  
264 < template<typename T> double NPTxyz<T>::getConservedQuantity(void){
265 <  
264 > double NPTxyz::getConservedQuantity(void){
265 >
266    double conservedQuantity;
267    double totalEnergy;
268    double thermostat_kinetic;
# Line 251 | Line 274 | template<typename T> double NPTxyz<T>::getConservedQua
274  
275    totalEnergy = tStats->getTotalE();
276  
277 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
277 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
278      (2.0 * eConvert);
279  
280    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 260 | Line 283 | template<typename T> double NPTxyz<T>::getConservedQua
283    info->matMul3(a, eta, b);
284    trEta = info->matTrace3(b);
285  
286 <  barostat_kinetic = NkBT * tb2 * trEta /
286 >  barostat_kinetic = NkBT * tb2 * trEta /
287      (2.0 * eConvert);
288 <  
289 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
288 >
289 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
290      eConvert;
291  
292    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
293      barostat_kinetic + barostat_potential;
294 <  
294 >
295   //   cout.width(8);
296   //   cout.precision(8);
297  
298 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
299 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
298 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
299 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
300   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
301  
302 <  return conservedQuantity;
303 <  
302 >  return conservedQuantity;
303 >
304   }
305 +
306 + char* NPTxyz::getAdditionalParameters(void){
307 +
308 +  sprintf(addParamBuffer,
309 +          "\t%G\t%G;"
310 +          "\t%G\t%G\t%G;"
311 +          "\t%G\t%G\t%G;"
312 +          "\t%G\t%G\t%G;",
313 +          chi, integralOfChidt,
314 +          eta[0][0], eta[0][1], eta[0][2],
315 +          eta[1][0], eta[1][1], eta[1][2],
316 +          eta[2][0], eta[2][1], eta[2][2]
317 +          );
318 +
319 +  return addParamBuffer;
320 + }

Diff Legend

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