ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/OOPSE/libmdtools/NPTf.cpp
(Generate patch)

Comparing:
trunk/OOPSE/libmdtools/NPTf.cpp (file contents), Revision 787 by mmeineke, Thu Sep 25 19:27:15 2003 UTC vs.
branches/new-templateless/OOPSE/libmdtools/NPTf.cpp (file contents), Revision 851 by mmeineke, Wed Nov 5 19:18:17 2003 UTC

# Line 1 | Line 1
1 < #include <cmath>
1 > #include <stdlib.h>
2 > #include <math.h>
3 > #include <string.h>
4 >
5   #include "Atom.hpp"
6   #include "SRI.hpp"
7   #include "AbstractClasses.hpp"
# Line 7 | Line 10
10   #include "Thermo.hpp"
11   #include "ReadWrite.hpp"
12   #include "Integrator.hpp"
13 < #include "simError.h"
13 > #include "simError.h"
14  
15   #ifdef IS_MPI
16   #include "mpiSimulation.hpp"
# Line 17 | Line 20
20   // modification of the Hoover algorithm:
21   //
22   //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
23 < //       Molec. Phys., 78, 533.
23 > //       Molec. Phys., 78, 533.
24   //
25   //           and
26 < //
26 > //
27   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28  
29 < template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
30 <  T( theInfo, the_ff )
29 > NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
30 >  NPT( theInfo, the_ff )
31   {
32 <  
32 >  GenericData* data;
33 >  double *etaArray;
34    int i,j;
35 <  
35 >
36    for(i = 0; i < 3; i++){
37      for (j = 0; j < 3; j++){
38 <      
38 >
39        eta[i][j] = 0.0;
40        oldEta[i][j] = 0.0;
41      }
42    }
43 +
44 +  // retrieve eta array from simInfo if it exists
45 +  data = info->getProperty(ETAVALUE_ID);
46 +  if(data != NULL){
47 +    
48 +    int test = data->getDarray(etaArray);
49 +    
50 +    if( test == 9 ){
51 +      
52 +      for(i = 0; i < 3; i++){
53 +        for (j = 0; j < 3; j++){
54 +          eta[i][j] = etaArray[3*i+j];
55 +          oldEta[i][j] = eta[i][j];
56 +        }
57 +      }    
58 +      delete[] etaArray;
59 +    }
60 +    else
61 +      std::cerr << "NPTf error: etaArray is not length 9 (actual = " << test
62 +                << ").\n"
63 +                << "            Simulation wil proceed with eta = 0;\n";
64 +  }
65   }
66  
67 < template<typename T> NPTf<T>::~NPTf() {
67 > NPTf::~NPTf() {
68  
69    // empty for now
70   }
71  
72 < template<typename T> void NPTf<T>::resetIntegrator() {
73 <  
72 > void NPTf::resetIntegrator() {
73 >
74    int i, j;
75 <  
75 >
76    for(i = 0; i < 3; i++)
77      for (j = 0; j < 3; j++)
78        eta[i][j] = 0.0;
79 <  
80 <  T::resetIntegrator();
79 >
80 >  NPT::resetIntegrator();
81   }
82  
83 < template<typename T> void NPTf<T>::evolveEtaA() {
84 <  
83 > void NPTf::evolveEtaA() {
84 >
85    int i, j;
86 <  
86 >
87    for(i = 0; i < 3; i ++){
88      for(j = 0; j < 3; j++){
89        if( i == j)
90 <        eta[i][j] += dt2 *  instaVol *
90 >        eta[i][j] += dt2 *  instaVol *
91            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
92        else
93          eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
94      }
95    }
96 <  
96 >
97    for(i = 0; i < 3; i++)
98      for (j = 0; j < 3; j++)
99        oldEta[i][j] = eta[i][j];
100   }
101  
102 < template<typename T> void NPTf<T>::evolveEtaB() {
103 <  
102 > void NPTf::evolveEtaB() {
103 >
104    int i,j;
105  
106    for(i = 0; i < 3; i++)
# Line 84 | Line 110 | template<typename T> void NPTf<T>::evolveEtaB() {
110    for(i = 0; i < 3; i ++){
111      for(j = 0; j < 3; j++){
112        if( i == j) {
113 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
113 >        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
114            (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
115        } else {
116          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
# Line 93 | Line 119 | template<typename T> void NPTf<T>::evolveEtaB() {
119    }
120   }
121  
122 < template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
122 > void NPTf::getVelScaleA(double sc[3], double vel[3]) {
123    int i,j;
124    double vScale[3][3];
125  
126    for (i = 0; i < 3; i++ ) {
127      for (j = 0; j < 3; j++ ) {
128        vScale[i][j] = eta[i][j];
129 <      
129 >
130        if (i == j) {
131 <        vScale[i][j] += chi;          
132 <      }              
131 >        vScale[i][j] += chi;
132 >      }
133      }
134    }
135 <  
135 >
136    info->matVecMul3( vScale, vel, sc );
137   }
138  
139 < template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
139 > void NPTf::getVelScaleB(double sc[3], int index ){
140    int i,j;
141    double myVel[3];
142    double vScale[3][3];
143  
144 + //   std::cerr << "velScaleB chi = " << chi << "\n";
145 +
146    for (i = 0; i < 3; i++ ) {
147      for (j = 0; j < 3; j++ ) {
148        vScale[i][j] = eta[i][j];
149 <      
149 >
150        if (i == j) {
151 <        vScale[i][j] += chi;          
152 <      }              
151 >        vScale[i][j] += chi;
152 >      }
153      }
154    }
155 <  
155 >
156    for (j = 0; j < 3; j++)
157      myVel[j] = oldVel[3*index + j];
158  
159 + //   std::cerr << "velScaleB = \n"
160 + //          << "[ " << vScale[0][0] << " , " << vScale[0][1] << " , " << vScale[0][2] << "]\n"
161 + //          << "[ " << vScale[1][0] << " , " << vScale[1][1] << " , " << vScale[1][2] << "]\n"
162 + //          << "[ " << vScale[2][0] << " , " << vScale[2][1] << " , " << vScale[2][2] << "]\n\n";
163 +
164 +
165 + //  std::cerr << "myVel " << index << " in => "
166 + //          << myVel[0] << ", " << myVel[1] << ", " << myVel[2] << "\n";
167 +
168    info->matVecMul3( vScale, myVel, sc );
169 +
170 + //  std::cerr << "sc " << index << " out => "
171 + //          << sc[0] << ", " << sc[1] << ", " << sc[2] << "\n";
172   }
173  
174 < template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
174 > void NPTf::getPosScale(double pos[3], double COM[3],
175                                                 int index, double sc[3]){
176    int j;
177    double rj[3];
# Line 142 | Line 182 | template<typename T> void NPTf<T>::getPosScale(double
182    info->matVecMul3( eta, rj, sc );
183   }
184  
185 < template<typename T> void NPTf<T>::scaleSimBox( void ){
185 > void NPTf::scaleSimBox( void ){
186  
187    int i,j,k;
188    double scaleMat[3][3];
189    double eta2ij;
190    double bigScale, smallScale, offDiagMax;
191    double hm[3][3], hmnew[3][3];
152 –  
192  
193  
194 +
195    // Scale the box after all the positions have been moved:
196 <  
196 >
197    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
198    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
199 <  
199 >
200    bigScale = 1.0;
201    smallScale = 1.0;
202    offDiagMax = 0.0;
203 <  
203 >
204    for(i=0; i<3; i++){
205      for(j=0; j<3; j++){
206 <      
206 >
207        // Calculate the matrix Product of the eta array (we only need
208        // the ij element right now):
209 <      
209 >
210        eta2ij = 0.0;
211        for(k=0; k<3; k++){
212          eta2ij += eta[i][k] * eta[k][j];
213        }
214 <      
214 >
215        scaleMat[i][j] = 0.0;
216        // identity matrix (see above):
217        if (i == j) scaleMat[i][j] = 1.0;
# Line 179 | Line 219 | template<typename T> void NPTf<T>::scaleSimBox( void )
219        scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
220  
221        if (i != j)
222 <        if (fabs(scaleMat[i][j]) > offDiagMax)
222 >        if (fabs(scaleMat[i][j]) > offDiagMax)
223            offDiagMax = fabs(scaleMat[i][j]);
224      }
225  
226      if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
227      if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
228    }
229 <  
229 >
230    if ((bigScale > 1.1) || (smallScale < 0.9)) {
231      sprintf( painCave.errMsg,
232               "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
# Line 218 | Line 258 | template<typename T> void NPTf<T>::scaleSimBox( void )
258    }
259   }
260  
261 < template<typename T> bool NPTf<T>::etaConverged() {
261 > bool NPTf::etaConverged() {
262    int i;
263    double diffEta, sumEta;
264  
265    sumEta = 0;
266    for(i = 0; i < 3; i++)
267 <    sumEta += pow(prevEta[i][i] - eta[i][i], 2);    
268 <  
267 >    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
268 >
269    diffEta = sqrt( sumEta / 3.0 );
270 <  
270 >
271    return ( diffEta <= etaTolerance );
272   }
273  
274 < template<typename T> double NPTf<T>::getConservedQuantity(void){
275 <  
274 > double NPTf::getConservedQuantity(void){
275 >
276    double conservedQuantity;
277    double totalEnergy;
278    double thermostat_kinetic;
# Line 244 | Line 284 | template<typename T> double NPTf<T>::getConservedQuant
284  
285    totalEnergy = tStats->getTotalE();
286  
287 <  thermostat_kinetic = fkBT * tt2 * chi * chi /
287 >  thermostat_kinetic = fkBT * tt2 * chi * chi /
288      (2.0 * eConvert);
289  
290    thermostat_potential = fkBT* integralOfChidt / eConvert;
# Line 253 | Line 293 | template<typename T> double NPTf<T>::getConservedQuant
293    info->matMul3(a, eta, b);
294    trEta = info->matTrace3(b);
295  
296 <  barostat_kinetic = NkBT * tb2 * trEta /
296 >  barostat_kinetic = NkBT * tb2 * trEta /
297      (2.0 * eConvert);
298 <  
299 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
298 >
299 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
300      eConvert;
301  
302    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
303      barostat_kinetic + barostat_potential;
304 <  
304 >
305   //   cout.width(8);
306   //   cout.precision(8);
307  
308 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
309 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
308 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
309 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
310   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
311  
312 <  return conservedQuantity;
313 <  
312 >  return conservedQuantity;
313 >
314   }
315 +
316 + char* NPTf::getAdditionalParameters(void){
317 +
318 +  sprintf(addParamBuffer,
319 +          "\t%G\t%G;"
320 +          "\t%G\t%G\t%G;"
321 +          "\t%G\t%G\t%G;"
322 +          "\t%G\t%G\t%G;",
323 +          chi, integralOfChidt,
324 +          eta[0][0], eta[0][1], eta[0][2],
325 +          eta[1][0], eta[1][1], eta[1][2],
326 +          eta[2][0], eta[2][1], eta[2][2]
327 +          );
328 +
329 +  return addParamBuffer;
330 + }

Diff Legend

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