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 851 by mmeineke, Wed Nov 5 19:18:17 2003 UTC

# Line 1 | Line 1
1 + #include <iostream>
2   #include <math.h>
3 +
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> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
30 <  T( theInfo, the_ff )
29 > NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
30 >  NPT( theInfo, the_ff )
31   {
32 +  GenericData* data;
33 +  double *etaArray;
34 +  int test;
35 +
36    eta = 0.0;
37    oldEta = 0.0;
38 +
39 +  // retrieve eta array from simInfo if it exists
40 +  data = info->getProperty(ETAVALUE_ID);
41 +  if(data != NULL){
42 +    
43 +    test = data->getDarray(etaArray);
44 +    
45 +    if( test == 9 ){
46 +      
47 +      eta = etaArray[0];
48 +      delete[] etaArray;
49 +    }
50 +    else
51 +      std::cerr << "NPTi error: etaArray is not length 9 (actual = " << test
52 +                << ").\n"
53 +                << "            Simulation wil proceed with eta = 0;\n";
54 +  }
55   }
56  
57 < template<typename T> NPTi<T>::~NPTi() {
57 > NPTi::~NPTi() {
58    //nothing for now
59   }
60  
61 < template<typename T> void NPTi<T>::resetIntegrator() {
61 > void NPTi::resetIntegrator() {
62    eta = 0.0;
63 <  T::resetIntegrator();
63 >  NPT::resetIntegrator();
64   }
65  
66 < template<typename T> void NPTi<T>::evolveEtaA() {
67 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
66 > void NPTi::evolveEtaA() {
67 >  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
68                   (p_convert*NkBT*tb2));
69    oldEta = eta;
70   }
71  
72 < template<typename T> void NPTi<T>::evolveEtaB() {
73 <  
72 > void NPTi::evolveEtaB() {
73 >
74    prevEta = eta;
75 <  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
75 >  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
76                   (p_convert*NkBT*tb2));
77   }
78  
79 < template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
79 > void NPTi::getVelScaleA(double sc[3], double vel[3]) {
80    int i;
81  
82    for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
83   }
84  
85 < template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
85 > void NPTi::getVelScaleB(double sc[3], int index ){
86    int i;
87  
88    for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * ( chi + eta );
89   }
90  
91  
92 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
92 > void NPTi::getPosScale(double pos[3], double COM[3],
93                                                 int index, double sc[3]){
94    int j;
95  
# Line 76 | Line 100 | template<typename T> void NPTi<T>::getPosScale(double
100      sc[j] *= eta;
101   }
102  
103 < template<typename T> void NPTi<T>::scaleSimBox( void ){
103 > void NPTi::scaleSimBox( void ){
104  
105    double scaleFactor;
106  
# Line 90 | Line 114 | template<typename T> void NPTi<T>::scaleSimBox( void )
114               );
115      painCave.isFatal = 1;
116      simError();
117 <  } else {        
118 <    info->scaleBox(scaleFactor);      
119 <  }  
117 >  } else {
118 >    info->scaleBox(scaleFactor);
119 >  }
120  
121   }
122  
123 < template<typename T> bool NPTi<T>::etaConverged() {
123 > bool NPTi::etaConverged() {
124  
125    return ( fabs(prevEta - eta) <= etaTolerance );
126   }
127  
128 < template<typename T> double NPTi<T>::getConservedQuantity(void){
128 > double NPTi::getConservedQuantity(void){
129  
130    double conservedQuantity;
131    double Energy;
# Line 109 | Line 133 | template<typename T> double NPTi<T>::getConservedQuant
133    double thermostat_potential;
134    double barostat_kinetic;
135    double barostat_potential;
136 <  
136 >
137    Energy = tStats->getTotalE();
138  
139 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
139 >  thermostat_kinetic = fkBT* tt2 * chi * chi /
140      (2.0 * eConvert);
141  
142    thermostat_potential = fkBT* integralOfChidt / eConvert;
143  
144  
145 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
145 >  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
146      (2.0 * eConvert);
147 <  
148 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
147 >
148 >  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
149      eConvert;
150  
151    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
152      barostat_kinetic + barostat_potential;
153 <  
153 >
154   //   cout.width(8);
155   //   cout.precision(8);
156  
157 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
158 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
157 > //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
158 > //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
159   //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
160 <  return conservedQuantity;
160 >  return conservedQuantity;
161   }
162 +
163 + char* NPTi::getAdditionalParameters(void){
164 +
165 +  sprintf(addParamBuffer,
166 +          "\t%G\t%G;"
167 +          "\t%G\t%0.0\t%0.0;"
168 +          "\t%0.0\t%G\t%0.0;"
169 +          "\t%0.0\t%0.0\t%G;",
170 +          chi, integralOfChidt,
171 +          eta, eta, eta
172 +          );
173 +
174 +  return addParamBuffer;
175 + }

Diff Legend

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