| 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 |
|
} |
| 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 |
|
} |
| 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 |
|
} |
| 238 |
|
// of freedom). |
| 239 |
|
fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp; |
| 240 |
|
|
| 241 |
– |
double conservedQuantity; |
| 242 |
– |
double totalEnergy; |
| 243 |
– |
double thermostat_kinetic; |
| 244 |
– |
double thermostat_potential; |
| 245 |
– |
double barostat_kinetic; |
| 246 |
– |
double barostat_potential; |
| 247 |
– |
double trEta; |
| 241 |
|
|
| 242 |
< |
totalEnergy = thermo.getTotalE(); |
| 242 |
> |
double totalEnergy = thermo.getTotalE(); |
| 243 |
|
|
| 244 |
< |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
| 252 |
< |
|
| 253 |
< |
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 |
|
} |
