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 |
|
} |
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 |
|
} |