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

Comparing branches/new_design/OOPSE-3.0/src/integrators/NPTi.cpp (file contents):
Revision 1821 by tim, Tue Nov 23 23:12:23 2004 UTC vs.
Revision 1822 by tim, Thu Dec 2 02:08:29 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/NPT.hpp"
5 > #include "primitives/Molecule.hpp"
6 > #include "utils/OOPSEConstant.hpp"
7   #include "utils/simError.h"
8  
9 < #ifdef IS_MPI
13 < #include "brains/mpiSimulation.hpp"
14 < #endif
9 > namespace oopse {
10  
11   // Basic isotropic thermostating and barostating via the Melchionna
12   // modification of the Hoover algorithm:
# Line 23 | Line 18 | NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
18   //
19   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
20  
21 < NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
27 <  T( theInfo, the_ff )
28 < {
29 <  GenericData* data;
30 <  DoubleVectorGenericData * etaValue;
31 <  vector<double> etaArray;
21 > NPTi::NPTi ( SimInfo *info) : NPT(info){
22  
23 <  eta = 0.0;
34 <  oldEta = 0.0;
23 >    Globals* globals = info_->getGlobals();
24  
36 –  if( theInfo->useInitXSstate ){
37 –    // retrieve eta from simInfo if
38 –    data = info->getPropertyByName(ETAVALUE_ID);
39 –    if(data){
40 –      etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
41 –      
42 –      if(etaValue){
43 –        eta = (*etaValue)[0];
44 –        oldEta = eta;
45 –      }
46 –    }
47 –  }
25   }
26  
50 – NPTi::~NPTi() {
51 –  //nothing for now
52 – }
53 –
27   void NPTi::evolveEtaA() {
28 <  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
29 <                 (p_convert*NkBT*tb2));
30 <  oldEta = eta;
28 >    eta += dt2 * ( instaVol * (instaPress - targetPressure) /
29 >         (OOPSEConstant::pressureConvert*NkBT*tb2));
30 >    oldEta = eta;
31   }
32  
33   void NPTi::evolveEtaB() {
34  
35 <  prevEta = eta;
36 <  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
37 <                 (p_convert*NkBT*tb2));
35 >    prevEta = eta;
36 >    eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
37 >         (OOPSEConstant::pressureConvert*NkBT*tb2));
38   }
39  
40   void NPTi::calcVelScale() {
41 <  vScale = chi + eta;
41 >    vScale = chi + eta;
42   }
43  
44   void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
# Line 78 | Line 51 | void NPTi::getPosScale(const Vector3d& pos, const Vect
51  
52  
53   void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
54 <                                               int index, Vector3d& sc){
54 >                           int index, Vector3d& sc){
55      /**@todo*/
56      sc  = oldPos[index] + pos/2.0 -COM;
57      sc *= eta;
# Line 86 | Line 59 | void NPTi::scaleSimBox(){
59  
60   void NPTi::scaleSimBox(){
61  
62 <  double scaleFactor;
62 >    double scaleFactor;
63  
64 <  scaleFactor = exp(dt*eta);
64 >    scaleFactor = exp(dt*eta);
65  
66 <  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
67 <    sprintf( painCave.errMsg,
66 >    if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
67 >        sprintf( painCave.errMsg,
68               "NPTi error: Attempting a Box scaling of more than 10 percent"
69               " check your tauBarostat, as it is probably too small!\n"
70               " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
71               );
72 <    painCave.isFatal = 1;
73 <    simError();
74 <  } else {
75 <    info->scaleBox(scaleFactor);
76 <  }
72 >        painCave.isFatal = 1;
73 >        simError();
74 >    } else {
75 >        Mat3x3d hmat = currentSnapshot_->getHmat();
76 >        hmat *= scaleFactor;
77 >        currentSnapshot_->setHmat(hmat);
78 >    }
79  
80   }
81  
82   bool NPTi::etaConverged() {
83  
84 <  return ( fabs(prevEta - eta) <= etaTolerance );
84 >    return ( fabs(prevEta - eta) <= etaTolerance );
85   }
86  
87   double NPTi::calcConservedQuantity(){
88  
89 <  double conservedQuantity;
90 <  double Energy;
91 <  double thermostat_kinetic;
92 <  double thermostat_potential;
93 <  double barostat_kinetic;
94 <  double barostat_potential;
89 >    double conservedQuantity;
90 >    double Energy;
91 >    double thermostat_kinetic;
92 >    double thermostat_potential;
93 >    double barostat_kinetic;
94 >    double barostat_potential;
95  
96 <  Energy = tStats->getTotalE();
96 >    Energy =thermo.getTotalE();
97  
98 <  thermostat_kinetic = fkBT* tt2 * chi * chi /
124 <    (2.0 * eConvert);
98 >    thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * OOPSEConstant::energyConvert);
99  
100 <  thermostat_potential = fkBT* integralOfChidt / eConvert;
100 >    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
101  
102  
103 <  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
130 <    (2.0 * eConvert);
103 >    barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);
104  
105 <  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
106 <    eConvert;
105 >    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
106 >        OOPSEConstant::energyConvert;
107  
108 <  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
109 <    barostat_kinetic + barostat_potential;
108 >    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
109 >        barostat_kinetic + barostat_potential;
110  
111 <  return conservedQuantity;
111 >    return conservedQuantity;
112   }
113  
114 + void NPTi::loadEta() {
115 +    Mat3x3d etaMat = currentSnapshot_->getEta();
116 +    eta = etaMat(0,0);
117 +    if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || etaMat.isDiagonal()) {
118 +        sprintf( painCave.errMsg,
119 +                 "NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
120 +        painCave.isFatal = 1;
121 +        simError();
122 +    }
123 + }
124 +
125 + void NPTi::saveEta() {
126 +    Mat3x3d etaMat(0.0);
127 +    etaMat(0, 0) = eta;
128 +    etaMat(1, 1) = eta;
129 +    etaMat(2, 2) = eta;
130 +    currentSnapshot_->setEta(etaMat);
131 + }
132 +
133 + }

Diff Legend

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