| 1 |
+ |
#include "brains/SimInfo.hpp" |
| 2 |
+ |
#include "brains/Thermo.hpp" |
| 3 |
|
#include "integrators/NPTxyz.hpp" |
| 4 |
+ |
#include "primitives/Molecule.hpp" |
| 5 |
+ |
#include "utils/OOPSEConstant.hpp" |
| 6 |
+ |
#include "utils/simError.h" |
| 7 |
|
|
| 8 |
|
// Basic non-isotropic thermostating and barostating via the Melchionna |
| 9 |
|
// modification of the Hoover algorithm: |
| 60 |
|
for(i = 0; i < 3; i ++){ |
| 61 |
|
for(j = 0; j < 3; j++){ |
| 62 |
|
if( i == j) { |
| 63 |
< |
eta(i, j) += dt2 * instaVol *(press(i, j) - targetPressure/p_convert) / (NkBT*tb2); |
| 63 |
> |
eta(i, j) += dt2 * instaVol *(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
| 64 |
|
} else { |
| 65 |
|
eta(i, j) = 0.0; |
| 66 |
|
} |
| 87 |
|
for(j = 0; j < 3; j++){ |
| 88 |
|
if( i == j) { |
| 89 |
|
eta(i, j) = oldEta(i, j) + dt2 * instaVol * |
| 90 |
< |
(press(i, j) - targetPressure/p_convert) / (NkBT*tb2); |
| 90 |
> |
(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
| 91 |
|
} else { |
| 92 |
|
eta(i, j) = 0.0; |
| 93 |
|
} |
| 149 |
|
double barostat_kinetic; |
| 150 |
|
double barostat_potential; |
| 151 |
|
double trEta; |
| 147 |
– |
Mat3x3d a; |
| 148 |
– |
Mat3x3d b; |
| 152 |
|
|
| 153 |
< |
totalEnergy = tStats->getTotalE(); |
| 153 |
> |
totalEnergy = thermo.getTotalE(); |
| 154 |
|
|
| 155 |
< |
thermostat_kinetic = fkBT * tt2 * chi * chi / |
| 153 |
< |
(2.0 * eConvert); |
| 155 |
> |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
| 156 |
|
|
| 157 |
< |
thermostat_potential = fkBT* integralOfChidt / eConvert; |
| 157 |
> |
thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
| 158 |
|
|
| 159 |
< |
transposeMat3(eta, a); |
| 160 |
< |
matMul3(a, eta, b); |
| 159 |
< |
trEta = matTrace3(b); |
| 159 |
> |
SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
| 160 |
> |
trEta = tmp.trace(); |
| 161 |
|
|
| 162 |
< |
barostat_kinetic = NkBT * tb2 * trEta / |
| 162 |
< |
(2.0 * eConvert); |
| 162 |
> |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
| 163 |
|
|
| 164 |
< |
barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / |
| 165 |
< |
eConvert; |
| 164 |
> |
barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
| 165 |
|
|
| 166 |
|
conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
| 167 |
|
barostat_kinetic + barostat_potential; |
| 224 |
|
painCave.isFatal = 1; |
| 225 |
|
simError(); |
| 226 |
|
} else { |
| 227 |
< |
info->getBoxM(hm); |
| 228 |
< |
matMul3(hm, scaleMat, hmnew); |
| 229 |
< |
info->setBoxM(hmnew); |
| 227 |
> |
|
| 228 |
> |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
| 229 |
> |
hmat = hmat *scaleMat; |
| 230 |
> |
currentSnapshot_->setHmat(hmat); |
| 231 |
|
} |
| 232 |
|
} |
| 233 |
|
|
