1 |
– |
#include <math.h> |
2 |
– |
|
3 |
– |
#include "math/MatVec3.h" |
4 |
– |
#include "primitives/Atom.hpp" |
5 |
– |
#include "primitives/SRI.hpp" |
6 |
– |
#include "primitives/AbstractClasses.hpp" |
1 |
|
#include "brains/SimInfo.hpp" |
8 |
– |
#include "UseTheForce/ForceFields.hpp" |
2 |
|
#include "brains/Thermo.hpp" |
3 |
< |
#include "io/ReadWrite.hpp" |
4 |
< |
#include "integrators/Integrator.hpp" |
3 |
> |
#include "integrators/NPTf.hpp" |
4 |
> |
#include "primitives/Molecule.hpp" |
5 |
> |
#include "utils/OOPSEConstant.hpp" |
6 |
|
#include "utils/simError.h" |
7 |
|
|
14 |
– |
#ifdef IS_MPI |
15 |
– |
#include "brains/mpiSimulation.hpp" |
16 |
– |
#endif |
17 |
– |
|
8 |
|
// Basic non-isotropic thermostating and barostating via the Melchionna |
9 |
|
// modification of the Hoover algorithm: |
10 |
|
// |
16 |
|
// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. |
17 |
|
|
18 |
|
NPTf::NPTf (SimInfo* info): NPT(info){ |
29 |
– |
GenericData* data; |
30 |
– |
DoubleVectorGenericData * etaValue; |
31 |
– |
int i,j; |
32 |
– |
|
33 |
– |
for(i = 0; i < 3; i++){ |
34 |
– |
for (j = 0; j < 3; j++){ |
35 |
– |
|
36 |
– |
eta(i, j) = 0.0; |
37 |
– |
oldEta(i, j) = 0.0; |
38 |
– |
} |
39 |
– |
} |
40 |
– |
|
41 |
– |
|
42 |
– |
if( theInfo->useInitXSstate ){ |
43 |
– |
// retrieve eta array from simInfo if it exists |
44 |
– |
data = info->getPropertyByName(ETAVALUE_ID); |
45 |
– |
if(data){ |
46 |
– |
etaValue = dynamic_cast<DoubleVectorGenericData*>(data); |
47 |
– |
|
48 |
– |
if(etaValue){ |
49 |
– |
|
50 |
– |
for(i = 0; i < 3; i++){ |
51 |
– |
for (j = 0; j < 3; j++){ |
52 |
– |
eta(i, j) = (*etaValue)[3*i+j]; |
53 |
– |
oldEta(i, j) = eta(i, j); |
54 |
– |
} |
55 |
– |
} |
56 |
– |
} |
57 |
– |
} |
58 |
– |
} |
19 |
|
|
20 |
|
} |
21 |
|
|
62 |
– |
NPTf::~NPTf() { |
22 |
|
|
64 |
– |
// empty for now |
65 |
– |
} |
66 |
– |
|
23 |
|
void NPTf::evolveEtaA() { |
24 |
|
|
25 |
|
int i, j; |
27 |
|
for(i = 0; i < 3; i ++){ |
28 |
|
for(j = 0; j < 3; j++){ |
29 |
|
if( i == j) { |
30 |
< |
eta(i, j) += dt2 * instaVol * (press(i, j) - targetPressure/p_convert) / (NkBT*tb2); |
30 |
> |
eta(i, j) += dt2 * instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
31 |
|
} else { |
32 |
|
eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2); |
33 |
|
} |
57 |
|
for(j = 0; j < 3; j++){ |
58 |
|
if( i == j) { |
59 |
|
eta(i, j) = oldEta(i, j) + dt2 * instaVol * |
60 |
< |
(press(i, j) - targetPressure/p_convert) / (NkBT*tb2); |
60 |
> |
(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
61 |
|
} else { |
62 |
|
eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2); |
63 |
|
} |
80 |
|
} |
81 |
|
} |
82 |
|
|
83 |
< |
void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel);{ |
83 |
> |
void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){ |
84 |
|
sc = vScale * vel; |
85 |
|
} |
86 |
|
|
88 |
|
sc = vScale * oldVel[index]; |
89 |
|
} |
90 |
|
|
91 |
< |
void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, |
136 |
< |
int index, Vector3d& sc) { |
91 |
> |
void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) { |
92 |
|
|
93 |
|
/**@todo */ |
94 |
< |
Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM; |
94 |
> |
Vector3d rj = (oldPos[index] + pos)/2.0 -COM; |
95 |
|
sc = eta * rj; |
96 |
|
} |
97 |
|
|
181 |
|
painCave.isFatal = 1; |
182 |
|
simError(); |
183 |
|
} else { |
184 |
< |
info->getBoxM(hm); |
185 |
< |
matMul3(hm, scaleMat, hmnew); |
186 |
< |
info->setBoxM(hmnew); |
184 |
> |
|
185 |
> |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
186 |
> |
hmat = hmat *scaleMat; |
187 |
> |
currentSnapshot_->setHmat(hmat); |
188 |
> |
|
189 |
|
} |
190 |
|
} |
191 |
|
|
213 |
|
double barostat_potential; |
214 |
|
double trEta; |
215 |
|
|
216 |
< |
totalEnergy = tStats->getTotalE(); |
216 |
> |
totalEnergy = thermo.getTotalE(); |
217 |
|
|
218 |
< |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * eConvert); |
218 |
> |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
219 |
|
|
220 |
< |
thermostat_potential = fkBT* integralOfChidt / eConvert; |
220 |
> |
thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
221 |
|
|
222 |
< |
trEta = (eta.transpose() * eta).trace(); |
222 |
> |
SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
223 |
> |
trEta = tmp.trace(); |
224 |
|
|
225 |
< |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * eConvert); |
225 |
> |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
226 |
|
|
227 |
< |
barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /eConvert; |
227 |
> |
barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
228 |
|
|
229 |
|
conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
230 |
|
barostat_kinetic + barostat_potential; |