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

Comparing trunk/src/integrators/NPrT.cpp (file contents):
Revision 536 by tim, Thu May 19 04:28:26 2005 UTC vs.
Revision 538 by tim, Thu May 19 21:31:23 2005 UTC

# Line 50 | Line 50 | namespace oopse {
50   namespace oopse {
51    NPrT::NPrT(SimInfo* info) : NPT(info) {
52      Globals* simParams = info_->getSimParams();
53 <    if (!simParams->haveTargetStress())
53 >    if (!simParams->haveSurfaceTension()) {
54        sprintf(painCave.errMsg,
55                "If you use the NPT integrator, you must set tauBarostat.\n");
56        painCave.severity = OOPSE_ERROR;
57        painCave.isFatal = 1;
58        simError();
59      } else {
60 <      targetStress= simParams->getTargetStress();
60 >      surfaceTension= simParams->getSurfaceTension();
61      }
62  
63    }
64    void NPrT::evolveEtaA() {
65 +    Mat3x3d hmat = currentSnapshot_->getHmat();
66 +    double hz = hmat(2, 2);
67 +    double Axy = hmat(0,0) * hmat(1, 1);
68      double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
69      double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert);
70 <    eta(0,0) -= Axy * (sx - targetStress) / (NkBT*tb2);
71 <    eta(1,1) -= Axy * (sy - targetStress) / (NkBT*tb2);
70 >    eta(0,0) -= Axy * (sx - surfaceTension) / (NkBT*tb2);
71 >    eta(1,1) -= Axy * (sy - surfaceTension) / (NkBT*tb2);
72      eta(2,2) += dt2 *  instaVol * (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
73      oldEta = eta;  
74    }
75  
76    void NPrT::evolveEtaB() {
77 <
77 >    Mat3x3d hmat = currentSnapshot_->getHmat();
78 >    double hz = hmat(2, 2);
79 >    double Axy = hmat(0,0) * hmat(1, 1);
80      prevEta = eta;
81      double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
82      double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert);
83 <    eta(0,0) -= Axy * (sx -targetStress) / (NkBT*tb2);
84 <    eta(1,1) -= Axy * (sy -targetStress) / (NkBT*tb2);
83 >    eta(0,0) = oldEta(0, 0) - Axy * (sx -surfaceTension) / (NkBT*tb2);
84 >    eta(1,1) = oldEta(1, 1) - Axy * (sy -surfaceTension) / (NkBT*tb2);
85      eta(2,2) = oldEta(2, 2) + dt2 *  instaVol *
86              (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
87    }
# Line 233 | Line 238 | namespace oopse {
238      // of freedom).  
239      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
240      
236 –    double conservedQuantity;
237 –    double totalEnergy;
238 –    double thermostat_kinetic;
239 –    double thermostat_potential;
240 –    double barostat_kinetic;
241 –    double barostat_potential;
242 –    double trEta;
241  
242 <    totalEnergy = thermo.getTotalE();
245 <
246 <    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
242 >    double totalEnergy = thermo.getTotalE();
243  
244 <    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
244 >    double thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
245  
246 +    double thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
247 +
248      SquareMatrix<double, 3> tmp = eta.transpose() * eta;
249 <    trEta = tmp.trace();
249 >    double trEta = tmp.trace();
250      
251 <    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
251 >    double barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
252  
253 <    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
253 >    double barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
254  
255 <    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
256 <      barostat_kinetic + barostat_potential;
255 >    Mat3x3d hmat = currentSnapshot_->getHmat();
256 >    double hz = hmat(2, 2);
257 >    double area = hmat(0,0) * hmat(1, 1);
258  
259 +    double conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
260 +      barostat_kinetic + barostat_potential - surfaceTension * area;
261 +
262      return conservedQuantity;
263  
264    }

Diff Legend

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