1 |
#include "brains/SimInfo.hpp" |
2 |
#include "brains/Thermo.hpp" |
3 |
#include "integrators/IntegratorCreator.hpp" |
4 |
#include "integrators/NPTxyz.hpp" |
5 |
#include "primitives/Molecule.hpp" |
6 |
#include "utils/OOPSEConstant.hpp" |
7 |
#include "utils/simError.h" |
8 |
|
9 |
// Basic non-isotropic thermostating and barostating via the Melchionna |
10 |
// modification of the Hoover algorithm: |
11 |
// |
12 |
// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993, |
13 |
// Molec. Phys., 78, 533. |
14 |
// |
15 |
// and |
16 |
// |
17 |
// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. |
18 |
|
19 |
namespace oopse { |
20 |
|
21 |
static IntegratorBuilder<NPTxyz>* NPTxyzCreator = new IntegratorBuilder<NPTxyz>("NPTxyz"); |
22 |
|
23 |
double NPTxyz::calcConservedQuantity(){ |
24 |
|
25 |
double conservedQuantity; |
26 |
double totalEnergy; |
27 |
double thermostat_kinetic; |
28 |
double thermostat_potential; |
29 |
double barostat_kinetic; |
30 |
double barostat_potential; |
31 |
double trEta; |
32 |
|
33 |
totalEnergy = thermo.getTotalE(); |
34 |
|
35 |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
36 |
|
37 |
thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
38 |
|
39 |
SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
40 |
trEta = tmp.trace(); |
41 |
|
42 |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
43 |
|
44 |
barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
45 |
|
46 |
conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
47 |
barostat_kinetic + barostat_potential; |
48 |
|
49 |
|
50 |
return conservedQuantity; |
51 |
|
52 |
} |
53 |
|
54 |
|
55 |
void NPTxyz::scaleSimBox(){ |
56 |
|
57 |
int i,j,k; |
58 |
Mat3x3d scaleMat; |
59 |
double eta2ij, scaleFactor; |
60 |
double bigScale, smallScale, offDiagMax; |
61 |
Mat3x3d hm; |
62 |
Mat3x3d hmnew; |
63 |
|
64 |
|
65 |
|
66 |
// Scale the box after all the positions have been moved: |
67 |
|
68 |
// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
69 |
// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
70 |
|
71 |
bigScale = 1.0; |
72 |
smallScale = 1.0; |
73 |
offDiagMax = 0.0; |
74 |
|
75 |
for(i=0; i<3; i++){ |
76 |
for(j=0; j<3; j++){ |
77 |
scaleMat(i, j) = 0.0; |
78 |
if(i==j) scaleMat(i, j) = 1.0; |
79 |
} |
80 |
} |
81 |
|
82 |
for(i=0;i<3;i++){ |
83 |
|
84 |
// calculate the scaleFactors |
85 |
|
86 |
scaleFactor = exp(dt*eta(i, i)); |
87 |
|
88 |
scaleMat(i, i) = scaleFactor; |
89 |
|
90 |
if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i); |
91 |
if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i); |
92 |
} |
93 |
|
94 |
if ((bigScale > 1.1) || (smallScale < 0.9)) { |
95 |
sprintf( painCave.errMsg, |
96 |
"NPTxyz error: Attempting a Box scaling of more than 10 percent.\n" |
97 |
" Check your tauBarostat, as it is probably too small!\n\n" |
98 |
" scaleMat = [%lf\t%lf\t%lf]\n" |
99 |
" [%lf\t%lf\t%lf]\n" |
100 |
" [%lf\t%lf\t%lf]\n", |
101 |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
102 |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
103 |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2)); |
104 |
painCave.isFatal = 1; |
105 |
simError(); |
106 |
} else { |
107 |
|
108 |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
109 |
hmat = hmat *scaleMat; |
110 |
currentSnapshot_->setHmat(hmat); |
111 |
} |
112 |
} |
113 |
|
114 |
void NPTf::loadEta() { |
115 |
eta= currentSnapshot_->getEta(); |
116 |
} |
117 |
|
118 |
} |