ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/integrators/NPTi.cpp
(Generate patch)

Comparing:
trunk/OOPSE-2.0/src/integrators/NPTi.cpp (file contents), Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/new_design/OOPSE-2.0/src/integrators/NPTi.cpp (file contents), Revision 1837 by tim, Thu Dec 2 22:15:31 2004 UTC

# Line 1 | Line 1
1 < #include <math.h>
2 < #include "primitives/Atom.hpp"
3 < #include "primitives/SRI.hpp"
4 < #include "primitives/AbstractClasses.hpp"
1 > #include "NPTi.hpp"
2   #include "brains/SimInfo.hpp"
6 – #include "UseTheForce/ForceFields.hpp"
3   #include "brains/Thermo.hpp"
4 < #include "io/ReadWrite.hpp"
5 < #include "integrators/Integrator.hpp"
4 > #include "integrators/IntegratorCreator.hpp"
5 > #include "integrators/NPT.hpp"
6 > #include "primitives/Molecule.hpp"
7 > #include "utils/OOPSEConstant.hpp"
8   #include "utils/simError.h"
9  
10 < #ifdef IS_MPI
13 < #include "brains/mpiSimulation.hpp"
14 < #endif
10 > namespace oopse {
11  
12 + static IntegratorBuilder<NPTi>* NPTiCreator = new IntegratorBuilder<NPTi>("NPTi");
13 +
14   // Basic isotropic thermostating and barostating via the Melchionna
15   // modification of the Hoover algorithm:
16   //
# Line 23 | Line 21 | template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo,
21   //
22   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
23  
24 < 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;
24 > NPTi::NPTi ( SimInfo *info) : NPT(info){
25  
26 <  eta = 0.0;
34 <  oldEta = 0.0;
26 >    Globals* globals = info_->getGlobals();
27  
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 –  }
28   }
29  
30 < template<typename T> NPTi<T>::~NPTi() {
31 <  //nothing for now
30 > void NPTi::evolveEtaA() {
31 >    eta += dt2 * ( instaVol * (instaPress - targetPressure) /
32 >         (OOPSEConstant::pressureConvert*NkBT*tb2));
33 >    oldEta = eta;
34   }
35  
36 < template<typename T> void NPTi<T>::resetIntegrator() {
56 <  eta = 0.0;
57 <  T::resetIntegrator();
58 < }
36 > void NPTi::evolveEtaB() {
37  
38 < template<typename T> void NPTi<T>::evolveEtaA() {
39 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
40 <                 (p_convert*NkBT*tb2));
63 <  oldEta = eta;
38 >    prevEta = eta;
39 >    eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
40 >         (OOPSEConstant::pressureConvert*NkBT*tb2));
41   }
42  
43 < template<typename T> void NPTi<T>::evolveEtaB() {
44 <
68 <  prevEta = eta;
69 <  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
70 <                 (p_convert*NkBT*tb2));
43 > void NPTi::calcVelScale() {
44 >    vScale = chi + eta;
45   }
46  
47 < template<typename T> void NPTi<T>::calcVelScale(void) {
48 <  vScale = chi + eta;
47 > void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
48 >    sc = vel * vScale;
49   }
50  
51 < template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
52 <  int i;
79 <
80 <  for(i=0; i<3; i++) sc[i] = vel[i] * vScale;
51 > void NPTi::getVelScaleB(Vector3d& sc, int index ){
52 >    sc = oldVel[index] * vScale;    
53   }
54  
83 – template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
84 –  int i;
55  
56 <  for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * vScale;
56 > void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
57 >                           int index, Vector3d& sc){
58 >    /**@todo*/
59 >    sc  = oldPos[index] + pos/2.0 -COM;
60 >    sc *= eta;
61   }
62  
63 + void NPTi::scaleSimBox(){
64  
65 < template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
91 <                                               int index, double sc[3]){
92 <  int j;
65 >    double scaleFactor;
66  
67 <  for(j=0; j<3; j++)
95 <    sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
67 >    scaleFactor = exp(dt*eta);
68  
69 <  for(j=0; j<3; j++)
70 <    sc[j] *= eta;
99 < }
100 <
101 < template<typename T> void NPTi<T>::scaleSimBox( void ){
102 <
103 <  double scaleFactor;
104 <
105 <  scaleFactor = exp(dt*eta);
106 <
107 <  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
108 <    sprintf( painCave.errMsg,
69 >    if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
70 >        sprintf( painCave.errMsg,
71               "NPTi error: Attempting a Box scaling of more than 10 percent"
72               " check your tauBarostat, as it is probably too small!\n"
73               " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
74               );
75 <    painCave.isFatal = 1;
76 <    simError();
77 <  } else {
78 <    info->scaleBox(scaleFactor);
79 <  }
75 >        painCave.isFatal = 1;
76 >        simError();
77 >    } else {
78 >        Mat3x3d hmat = currentSnapshot_->getHmat();
79 >        hmat *= scaleFactor;
80 >        currentSnapshot_->setHmat(hmat);
81 >    }
82  
83   }
84  
85 < template<typename T> bool NPTi<T>::etaConverged() {
85 > bool NPTi::etaConverged() {
86  
87 <  return ( fabs(prevEta - eta) <= etaTolerance );
87 >    return ( fabs(prevEta - eta) <= etaTolerance );
88   }
89  
90 < template<typename T> double NPTi<T>::getConservedQuantity(void){
90 > double NPTi::calcConservedQuantity(){
91  
92 <  double conservedQuantity;
93 <  double Energy;
94 <  double thermostat_kinetic;
95 <  double thermostat_potential;
96 <  double barostat_kinetic;
97 <  double barostat_potential;
92 >    double conservedQuantity;
93 >    double Energy;
94 >    double thermostat_kinetic;
95 >    double thermostat_potential;
96 >    double barostat_kinetic;
97 >    double barostat_potential;
98  
99 <  Energy = tStats->getTotalE();
99 >    Energy =thermo.getTotalE();
100  
101 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
138 <    (2.0 * eConvert);
101 >    thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * OOPSEConstant::energyConvert);
102  
103 <  thermostat_potential = fkBT* integralOfChidt / eConvert;
103 >    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
104  
105  
106 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
144 <    (2.0 * eConvert);
106 >    barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);
107  
108 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
109 <    eConvert;
108 >    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
109 >        OOPSEConstant::energyConvert;
110  
111 <  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
112 <    barostat_kinetic + barostat_potential;
111 >    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
112 >        barostat_kinetic + barostat_potential;
113  
114 < //   cout.width(8);
115 < //   cout.precision(8);
114 >    return conservedQuantity;
115 > }
116  
117 < //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
118 < //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
119 < //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
120 <  return conservedQuantity;
117 > void NPTi::loadEta() {
118 >    Mat3x3d etaMat = currentSnapshot_->getEta();
119 >    eta = etaMat(0,0);
120 >    if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || etaMat.isDiagonal()) {
121 >        sprintf( painCave.errMsg,
122 >                 "NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
123 >        painCave.isFatal = 1;
124 >        simError();
125 >    }
126   }
127  
128 < template<typename T> string NPTi<T>::getAdditionalParameters(void){
129 <  string parameters;
130 <  const int BUFFERSIZE = 2000; // size of the read buffer
131 <  char buffer[BUFFERSIZE];
128 > void NPTi::saveEta() {
129 >    Mat3x3d etaMat(0.0);
130 >    etaMat(0, 0) = eta;
131 >    etaMat(1, 1) = eta;
132 >    etaMat(2, 2) = eta;
133 >    currentSnapshot_->setEta(etaMat);
134 > }
135  
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 –
136   }

Diff Legend

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