ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/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.
branches/new-templateless/OOPSE/libmdtools/NPTi.cpp (file contents), Revision 849 by mmeineke, Fri Oct 31 21:06:47 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 )
26 > NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
27 >  Integrator( 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() {
50 > NPTi::~NPTi() {
51    //nothing for now
52   }
53  
54 < template<typename T> void NPTi<T>::resetIntegrator() {
54 > void NPTi::resetIntegrator() {
55    eta = 0.0;
56 <  T::resetIntegrator();
56 >  Integrator::resetIntegrator();
57   }
58  
59 < template<typename T> void NPTi<T>::evolveEtaA() {
60 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
59 > void NPTi::evolveEtaA() {
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 <  
65 > void NPTi::evolveEtaB() {
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  
72 < template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
72 > void NPTi::getVelScaleA(double sc[3], double vel[3]) {
73    int i;
74  
75    for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
76   }
77  
78 < template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
78 > void NPTi::getVelScaleB(double sc[3], int index ){
79    int i;
80  
81    for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * ( chi + eta );
82   }
83  
84  
85 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
85 > void NPTi::getPosScale(double pos[3], double COM[3],
86                                                 int index, double sc[3]){
87    int j;
88  
# Line 76 | Line 93 | template<typename T> void NPTi<T>::getPosScale(double
93      sc[j] *= eta;
94   }
95  
96 < template<typename T> void NPTi<T>::scaleSimBox( void ){
96 > void NPTi::scaleSimBox( void ){
97  
98    double scaleFactor;
99  
# 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  
116 < template<typename T> bool NPTi<T>::etaConverged() {
116 > bool NPTi::etaConverged() {
117  
118    return ( fabs(prevEta - eta) <= etaTolerance );
119   }
120  
121 < template<typename T> double NPTi<T>::getConservedQuantity(void){
121 > double NPTi::getConservedQuantity(void){
122  
123    double conservedQuantity;
124    double Energy;
# 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 + string NPTi::getAdditionalParameters(void){
157 +  string parameters;
158 +  const int BUFFERSIZE = 2000; // size of the read buffer
159 +  char buffer[BUFFERSIZE];
160 +
161 +  sprintf(buffer,"\t%lf\t%lf;", chi, integralOfChidt);
162 +  parameters += buffer;
163 +
164 +  sprintf(buffer,"\t%lf\t0\t0;", eta);
165 +  parameters += buffer;
166 +
167 +  sprintf(buffer,"\t0\t%lf\t0;", eta);
168 +  parameters += buffer;
169 +
170 +  sprintf(buffer,"\t0\t0\t%lf;", eta);
171 +  parameters += buffer;
172 +
173 +  return parameters;
174 +
175 + }

Diff Legend

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