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 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> 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 +  // retrieve eta from simInfo if
37 +  data = info->getProperty(ETAVALUE_ID);
38 +  if(data){
39 +    etaValue = dynamic_cast<DoubleArrayData*>(data);
40 +
41 +    if(etaValue){
42 +      etaArray = etaValue->getData();
43 +      eta = etaArray[0];
44 +      oldEta = eta;
45 +    }
46 +  }
47 +
48   }
49  
50   template<typename T> NPTi<T>::~NPTi() {
# Line 40 | Line 57 | template<typename T> void NPTi<T>::evolveEtaA() {
57   }
58  
59   template<typename T> void NPTi<T>::evolveEtaA() {
60 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
60 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
61                   (p_convert*NkBT*tb2));
62    oldEta = eta;
63   }
64  
65   template<typename T> void NPTi<T>::evolveEtaB() {
66 <  
66 >
67    prevEta = eta;
68 <  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
68 >  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
69                   (p_convert*NkBT*tb2));
70   }
71  
# Line 65 | Line 82 | template<typename T> void NPTi<T>::getPosScale(double
82   }
83  
84  
85 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
85 > template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
86                                                 int index, double sc[3]){
87    int j;
88  
# Line 90 | Line 107 | template<typename T> void NPTi<T>::scaleSimBox( void )
107               );
108      painCave.isFatal = 1;
109      simError();
110 <  } else {        
111 <    info->scaleBox(scaleFactor);      
112 <  }  
110 >  } else {
111 >    info->scaleBox(scaleFactor);
112 >  }
113  
114   }
115  
# Line 109 | Line 126 | template<typename T> double NPTi<T>::getConservedQuant
126    double thermostat_potential;
127    double barostat_kinetic;
128    double barostat_potential;
129 <  
129 >
130    Energy = tStats->getTotalE();
131  
132 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
132 >  thermostat_kinetic = fkBT* tt2 * chi * chi /
133      (2.0 * eConvert);
134  
135    thermostat_potential = fkBT* integralOfChidt / eConvert;
136  
137  
138 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
138 >  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
139      (2.0 * eConvert);
140 <  
141 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
140 >
141 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
142      eConvert;
143  
144    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
145      barostat_kinetic + barostat_potential;
146 <  
146 >
147   //   cout.width(8);
148   //   cout.precision(8);
149  
150 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
151 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
150 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
151 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
152   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
153 <  return conservedQuantity;
153 >  return conservedQuantity;
154   }
155 +
156 + template<typename T> string NPTi<T>::getAdditionalParameters(void){
157 +  string parameters;
158 +  const int BUFFERSIZE = 2000; // size of the read buffer
159 +  char buffer[BUFFERSIZE];
160 +
161 +  sprintf(buffer,"\t%g\t%g;", chi, integralOfChidt);
162 +  parameters += buffer;
163 +
164 +  sprintf(buffer,"\t%g\t0\t0;", eta);
165 +  parameters += buffer;
166 +
167 +  sprintf(buffer,"\t0\t%g\t0;", eta);
168 +  parameters += buffer;
169 +
170 +  sprintf(buffer,"\t0\t0\t%g;", eta);
171 +  parameters += buffer;
172 +
173 +  return parameters;
174 +
175 + }

Diff Legend

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