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 855 by mmeineke, Thu Nov 6 22:01:37 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  
# Line 65 | Line 83 | template<typename T> void NPTi<T>::getPosScale(double
83   }
84  
85  
86 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
86 > template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
87                                                 int index, double sc[3]){
88    int j;
89  
# Line 90 | Line 108 | template<typename T> void NPTi<T>::scaleSimBox( void )
108               );
109      painCave.isFatal = 1;
110      simError();
111 <  } else {        
112 <    info->scaleBox(scaleFactor);      
113 <  }  
111 >  } else {
112 >    info->scaleBox(scaleFactor);
113 >  }
114  
115   }
116  
# Line 109 | Line 127 | template<typename T> double NPTi<T>::getConservedQuant
127    double thermostat_potential;
128    double barostat_kinetic;
129    double barostat_potential;
130 <  
130 >
131    Energy = tStats->getTotalE();
132  
133 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
133 >  thermostat_kinetic = fkBT* tt2 * chi * chi /
134      (2.0 * eConvert);
135  
136    thermostat_potential = fkBT* integralOfChidt / eConvert;
137  
138  
139 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
139 >  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
140      (2.0 * eConvert);
141 <  
142 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
141 >
142 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
143      eConvert;
144  
145    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
146      barostat_kinetic + barostat_potential;
147 <  
147 >
148   //   cout.width(8);
149   //   cout.precision(8);
150  
151 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
152 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
151 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
152 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
153   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
154 <  return conservedQuantity;
154 >  return conservedQuantity;
155   }
156 +
157 + template<typename T> string NPTi<T>::getAdditionalParameters(void){
158 +  string parameters;
159 +  const int BUFFERSIZE = 2000; // size of the read buffer
160 +  char buffer[BUFFERSIZE];
161 +
162 +  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
163 +  parameters += buffer;
164 +
165 +  sprintf(buffer,"\t%G\t0\t0;", eta);
166 +  parameters += buffer;
167 +
168 +  sprintf(buffer,"\t0\t%G\t0;", eta);
169 +  parameters += buffer;
170 +
171 +  sprintf(buffer,"\t0\t0\t%G;", eta);
172 +  parameters += buffer;
173 +
174 +  return parameters;
175 +
176 + }

Diff Legend

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