ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/NPTi.cpp
(Generate patch)

Comparing trunk/src/integrators/NPTi.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 49 | Line 49 | namespace oopse {
49  
50   namespace oopse {
51  
52 < // Basic isotropic thermostating and barostating via the Melchionna
53 < // modification of the Hoover algorithm:
54 < //
55 < //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
56 < //       Molec. Phys., 78, 533.
57 < //
58 < //           and
59 < //
60 < //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
52 >  // Basic isotropic thermostating and barostating via the Melchionna
53 >  // modification of the Hoover algorithm:
54 >  //
55 >  //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
56 >  //       Molec. Phys., 78, 533.
57 >  //
58 >  //           and
59 >  //
60 >  //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61  
62 < NPTi::NPTi ( SimInfo *info) : NPT(info){
62 >  NPTi::NPTi ( SimInfo *info) : NPT(info){
63  
64 < }
64 >  }
65  
66 < void NPTi::evolveEtaA() {
66 >  void NPTi::evolveEtaA() {
67      eta += dt2 * ( instaVol * (instaPress - targetPressure) /
68 <         (OOPSEConstant::pressureConvert*NkBT*tb2));
68 >                   (OOPSEConstant::pressureConvert*NkBT*tb2));
69      oldEta = eta;
70 < }
70 >  }
71  
72 < void NPTi::evolveEtaB() {
72 >  void NPTi::evolveEtaB() {
73  
74      prevEta = eta;
75      eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
76 <         (OOPSEConstant::pressureConvert*NkBT*tb2));
77 < }
76 >                           (OOPSEConstant::pressureConvert*NkBT*tb2));
77 >  }
78  
79 < void NPTi::calcVelScale() {
79 >  void NPTi::calcVelScale() {
80      vScale = chi + eta;
81 < }
81 >  }
82  
83 < void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
83 >  void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
84      sc = vel * vScale;
85 < }
85 >  }
86  
87 < void NPTi::getVelScaleB(Vector3d& sc, int index ){
87 >  void NPTi::getVelScaleB(Vector3d& sc, int index ){
88      sc = oldVel[index] * vScale;    
89 < }
89 >  }
90  
91  
92 < void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
93 <                           int index, Vector3d& sc){
92 >  void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
93 >                         int index, Vector3d& sc){
94      /**@todo*/
95 <    sc  = (oldPos[index] + pos)/2.0 -COM;
95 >    sc  = (oldPos[index] + pos)/(RealType)2.0 -COM;
96      sc *= eta;
97 < }
97 >  }
98  
99 < void NPTi::scaleSimBox(){
99 >  void NPTi::scaleSimBox(){
100  
101 <    double scaleFactor;
101 >    RealType scaleFactor;
102  
103      scaleFactor = exp(dt*eta);
104  
105      if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
106 <        sprintf( painCave.errMsg,
107 <             "NPTi error: Attempting a Box scaling of more than 10 percent"
108 <             " check your tauBarostat, as it is probably too small!\n"
109 <             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
110 <             );
111 <        painCave.isFatal = 1;
112 <        simError();
106 >      sprintf( painCave.errMsg,
107 >               "NPTi error: Attempting a Box scaling of more than 10 percent"
108 >               " check your tauBarostat, as it is probably too small!\n"
109 >               " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
110 >               );
111 >      painCave.isFatal = 1;
112 >      simError();
113      } else {
114 <        Mat3x3d hmat = currentSnapshot_->getHmat();
115 <        hmat *= scaleFactor;
116 <        currentSnapshot_->setHmat(hmat);
114 >      Mat3x3d hmat = currentSnapshot_->getHmat();
115 >      hmat *= scaleFactor;
116 >      currentSnapshot_->setHmat(hmat);
117      }
118  
119 < }
119 >  }
120  
121 < bool NPTi::etaConverged() {
121 >  bool NPTi::etaConverged() {
122  
123      return ( fabs(prevEta - eta) <= etaTolerance );
124 < }
124 >  }
125  
126 < double NPTi::calcConservedQuantity(){
126 >  RealType NPTi::calcConservedQuantity(){
127  
128      chi= currentSnapshot_->getChi();
129      integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
# Line 138 | Line 138 | double NPTi::calcConservedQuantity(){
138      // of freedom).  
139      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
140      
141 <    double conservedQuantity;
142 <    double Energy;
143 <    double thermostat_kinetic;
144 <    double thermostat_potential;
145 <    double barostat_kinetic;
146 <    double barostat_potential;
141 >    RealType conservedQuantity;
142 >    RealType Energy;
143 >    RealType thermostat_kinetic;
144 >    RealType thermostat_potential;
145 >    RealType barostat_kinetic;
146 >    RealType barostat_potential;
147  
148      Energy =thermo.getTotalE();
149  
# Line 155 | Line 155 | double NPTi::calcConservedQuantity(){
155      barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);
156  
157      barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
158 <        OOPSEConstant::energyConvert;
158 >      OOPSEConstant::energyConvert;
159  
160      conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
161 <        barostat_kinetic + barostat_potential;
161 >      barostat_kinetic + barostat_potential;
162      
163      return conservedQuantity;
164 < }
164 >  }
165  
166 < void NPTi::loadEta() {
166 >  void NPTi::loadEta() {
167      Mat3x3d etaMat = currentSnapshot_->getEta();
168      eta = etaMat(0,0);
169      //if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || !etaMat.isDiagonal()) {
# Line 172 | Line 172 | void NPTi::loadEta() {
172      //    painCave.isFatal = 1;
173      //    simError();
174      //}
175 < }
175 >  }
176  
177 < void NPTi::saveEta() {
177 >  void NPTi::saveEta() {
178      Mat3x3d etaMat(0.0);
179      etaMat(0, 0) = eta;
180      etaMat(1, 1) = eta;
181      etaMat(2, 2) = eta;
182      currentSnapshot_->setEta(etaMat);
183 +  }
184   }
184 –
185 – }

Diff Legend

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