ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPTi.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPTi.cpp (file contents):
Revision 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 857 by mmeineke, Fri Nov 7 17:09:48 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> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
27    T( theInfo, the_ff )
28   {
29 +  GenericData* data;
30 +  DoubleArrayData * etaValue;
31 +  vector<double> etaArray;
32 +
33    eta = 0.0;
34    oldEta = 0.0;
35 +
36 +  if( theInfo->useInitXSstate ){
37 +    // retrieve eta from simInfo if
38 +    data = info->getProperty(ETAVALUE_ID);
39 +    if(data){
40 +      etaValue = dynamic_cast<DoubleArrayData*>(data);
41 +      
42 +      if(etaValue){
43 +        etaArray = etaValue->getData();
44 +        eta = etaArray[0];
45 +        oldEta = eta;
46 +      }
47 +    }
48 +  }
49   }
50  
51   template<typename T> NPTi<T>::~NPTi() {
# Line 40 | Line 58 | template<typename T> void NPTi<T>::evolveEtaA() {
58   }
59  
60   template<typename T> void NPTi<T>::evolveEtaA() {
61 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
61 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
62                   (p_convert*NkBT*tb2));
63    oldEta = eta;
64   }
65  
66   template<typename T> void NPTi<T>::evolveEtaB() {
67 <  
67 >
68    prevEta = eta;
69 <  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
69 >  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
70                   (p_convert*NkBT*tb2));
71   }
72  
73 + template<typename T> void NPTi<T>::calcVelScale(void) {
74 +  vScale = chi + eta;
75 + }
76 +
77   template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
78    int i;
79  
80 <  for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
80 >  for(i=0; i<3; i++) sc[i] = vel[i] * vScale;
81   }
82  
83   template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
84    int i;
85  
86 <  for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * ( chi + eta );
86 >  for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * vScale;
87   }
88  
89  
90 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
90 > template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
91                                                 int index, double sc[3]){
92    int j;
93  
# Line 90 | Line 112 | template<typename T> void NPTi<T>::scaleSimBox( void )
112               );
113      painCave.isFatal = 1;
114      simError();
115 <  } else {        
116 <    info->scaleBox(scaleFactor);      
117 <  }  
115 >  } else {
116 >    info->scaleBox(scaleFactor);
117 >  }
118  
119   }
120  
# Line 109 | Line 131 | template<typename T> double NPTi<T>::getConservedQuant
131    double thermostat_potential;
132    double barostat_kinetic;
133    double barostat_potential;
134 <  
134 >
135    Energy = tStats->getTotalE();
136  
137 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
137 >  thermostat_kinetic = fkBT* tt2 * chi * chi /
138      (2.0 * eConvert);
139  
140    thermostat_potential = fkBT* integralOfChidt / eConvert;
141  
142  
143 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
143 >  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
144      (2.0 * eConvert);
145 <  
146 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
145 >
146 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
147      eConvert;
148  
149    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
150      barostat_kinetic + barostat_potential;
151 <  
151 >
152   //   cout.width(8);
153   //   cout.precision(8);
154  
155 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
156 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
155 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
156 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
157   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
158 <  return conservedQuantity;
158 >  return conservedQuantity;
159   }
160 +
161 + template<typename T> string NPTi<T>::getAdditionalParameters(void){
162 +  string parameters;
163 +  const int BUFFERSIZE = 2000; // size of the read buffer
164 +  char buffer[BUFFERSIZE];
165 +
166 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
167 +  parameters += buffer;
168 +
169 +  sprintf(buffer,"\t%G\t0\t0;", eta);
170 +  parameters += buffer;
171 +
172 +  sprintf(buffer,"\t0\t%G\t0;", eta);
173 +  parameters += buffer;
174 +
175 +  sprintf(buffer,"\t0\t0\t%G;", eta);
176 +  parameters += buffer;
177 +
178 +  return parameters;
179 +
180 + }

Diff Legend

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