57 |
|
painCave.isFatal = 1; |
58 |
|
simError(); |
59 |
|
} else { |
60 |
< |
surfaceTension= simParams->getSurfaceTension(); |
60 |
> |
surfaceTension= simParams->getSurfaceTension()* OOPSEConstant::surfaceTensorConvert * OOPSEConstant::energyConvert; |
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 - surfaceTension) / (NkBT*tb2); |
71 |
< |
eta(1,1) -= Axy * (sy - surfaceTension) / (NkBT*tb2); |
70 |
> |
eta(0,0) -= dt2* Axy * (sx - surfaceTension) / (NkBT*tb2); |
71 |
> |
eta(1,1) -= dt2* 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) = oldEta(0, 0) - Axy * (sx -surfaceTension) / (NkBT*tb2); |
84 |
< |
eta(1,1) = oldEta(1, 1) - Axy * (sy -surfaceTension) / (NkBT*tb2); |
83 |
> |
eta(0,0) = oldEta(0, 0) - dt2 * Axy * (sx -surfaceTension) / (NkBT*tb2); |
84 |
> |
eta(1,1) = oldEta(1, 1) - dt2 * Axy * (sy -surfaceTension) / (NkBT*tb2); |
85 |
|
eta(2,2) = oldEta(2, 2) + dt2 * instaVol * |
86 |
|
(press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
87 |
|
} |
175 |
|
double area = hmat(0,0) * hmat(1, 1); |
176 |
|
|
177 |
|
double conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
178 |
< |
barostat_kinetic + barostat_potential - surfaceTension * area; |
178 |
> |
barostat_kinetic + barostat_potential - surfaceTension * area/ OOPSEConstant::energyConvert; |
179 |
|
|
180 |
|
return conservedQuantity; |
181 |
|
|