1 |
/* |
2 |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
6 |
* redistribute this software in source and binary code form, provided |
7 |
* that the following conditions are met: |
8 |
* |
9 |
* 1. Acknowledgement of the program authors must be made in any |
10 |
* publication of scientific results based in part on use of the |
11 |
* program. An acceptable form of acknowledgement is citation of |
12 |
* the article in which the program was described (Matthew |
13 |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
* Parallel Simulation Engine for Molecular Dynamics," |
16 |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
* |
18 |
* 2. Redistributions of source code must retain the above copyright |
19 |
* notice, this list of conditions and the following disclaimer. |
20 |
* |
21 |
* 3. Redistributions in binary form must reproduce the above copyright |
22 |
* notice, this list of conditions and the following disclaimer in the |
23 |
* documentation and/or other materials provided with the |
24 |
* distribution. |
25 |
* |
26 |
* This software is provided "AS IS," without a warranty of any |
27 |
* kind. All express or implied conditions, representations and |
28 |
* warranties, including any implied warranty of merchantability, |
29 |
* fitness for a particular purpose or non-infringement, are hereby |
30 |
* excluded. The University of Notre Dame and its licensors shall not |
31 |
* be liable for any damages suffered by licensee as a result of |
32 |
* using, modifying or distributing the software or its |
33 |
* derivatives. In no event will the University of Notre Dame or its |
34 |
* licensors be liable for any lost revenue, profit or data, or for |
35 |
* direct, indirect, special, consequential, incidental or punitive |
36 |
* damages, however caused and regardless of the theory of liability, |
37 |
* arising out of the use of or inability to use software, even if the |
38 |
* University of Notre Dame has been advised of the possibility of |
39 |
* such damages. |
40 |
*/ |
41 |
|
42 |
#include "brains/SimInfo.hpp" |
43 |
#include "brains/Thermo.hpp" |
44 |
#include "integrators/IntegratorCreator.hpp" |
45 |
#include "integrators/NPrT.hpp" |
46 |
#include "primitives/Molecule.hpp" |
47 |
#include "utils/OOPSEConstant.hpp" |
48 |
#include "utils/simError.h" |
49 |
|
50 |
namespace oopse { |
51 |
NPrT::NPrT(SimInfo* info) : NPT(info) { |
52 |
Globals* simParams = info_->getSimParams(); |
53 |
if (!simParams->haveTargetStress()) |
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(); |
61 |
} |
62 |
|
63 |
} |
64 |
void NPrT::evolveEtaA() { |
65 |
double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert); |
66 |
double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert); |
67 |
eta(0,0) -= Axy * (sx - targetStress) / (NkBT*tb2); |
68 |
eta(1,1) -= Axy * (sy - targetStress) / (NkBT*tb2); |
69 |
eta(2,2) += dt2 * instaVol * (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
70 |
oldEta = eta; |
71 |
} |
72 |
|
73 |
void NPrT::evolveEtaB() { |
74 |
|
75 |
prevEta = eta; |
76 |
double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert); |
77 |
double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert); |
78 |
eta(0,0) -= Axy * (sx -targetStress) / (NkBT*tb2); |
79 |
eta(1,1) -= Axy * (sy -targetStress) / (NkBT*tb2); |
80 |
eta(2,2) = oldEta(2, 2) + dt2 * instaVol * |
81 |
(press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
82 |
} |
83 |
|
84 |
void NPrT::calcVelScale(){ |
85 |
|
86 |
for (int i = 0; i < 3; i++ ) { |
87 |
for (int j = 0; j < 3; j++ ) { |
88 |
vScale(i, j) = eta(i, j); |
89 |
|
90 |
if (i == j) { |
91 |
vScale(i, j) += chi; |
92 |
} |
93 |
} |
94 |
} |
95 |
} |
96 |
|
97 |
void NPrT::getVelScaleA(Vector3d& sc, const Vector3d& vel){ |
98 |
sc = vScale * vel; |
99 |
} |
100 |
|
101 |
void NPrT::getVelScaleB(Vector3d& sc, int index ) { |
102 |
sc = vScale * oldVel[index]; |
103 |
} |
104 |
|
105 |
void NPrT::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) { |
106 |
|
107 |
/**@todo */ |
108 |
Vector3d rj = (oldPos[index] + pos)/2.0 -COM; |
109 |
sc = eta * rj; |
110 |
} |
111 |
|
112 |
void NPrT::scaleSimBox(){ |
113 |
|
114 |
int i; |
115 |
int j; |
116 |
int k; |
117 |
Mat3x3d scaleMat; |
118 |
double eta2ij; |
119 |
double bigScale, smallScale, offDiagMax; |
120 |
Mat3x3d hm; |
121 |
Mat3x3d hmnew; |
122 |
|
123 |
|
124 |
|
125 |
// Scale the box after all the positions have been moved: |
126 |
|
127 |
// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
128 |
// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
129 |
|
130 |
bigScale = 1.0; |
131 |
smallScale = 1.0; |
132 |
offDiagMax = 0.0; |
133 |
|
134 |
for(i=0; i<3; i++){ |
135 |
for(j=0; j<3; j++){ |
136 |
|
137 |
// Calculate the matrix Product of the eta array (we only need |
138 |
// the ij element right now): |
139 |
|
140 |
eta2ij = 0.0; |
141 |
for(k=0; k<3; k++){ |
142 |
eta2ij += eta(i, k) * eta(k, j); |
143 |
} |
144 |
|
145 |
scaleMat(i, j) = 0.0; |
146 |
// identity matrix (see above): |
147 |
if (i == j) scaleMat(i, j) = 1.0; |
148 |
// Taylor expansion for the exponential truncated at second order: |
149 |
scaleMat(i, j) += dt*eta(i, j) + 0.5*dt*dt*eta2ij; |
150 |
|
151 |
|
152 |
if (i != j) |
153 |
if (fabs(scaleMat(i, j)) > offDiagMax) |
154 |
offDiagMax = fabs(scaleMat(i, j)); |
155 |
} |
156 |
|
157 |
if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i); |
158 |
if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i); |
159 |
} |
160 |
|
161 |
if ((bigScale > 1.01) || (smallScale < 0.99)) { |
162 |
sprintf( painCave.errMsg, |
163 |
"NPrT error: Attempting a Box scaling of more than 1 percent.\n" |
164 |
" Check your tauBarostat, as it is probably too small!\n\n" |
165 |
" scaleMat = [%lf\t%lf\t%lf]\n" |
166 |
" [%lf\t%lf\t%lf]\n" |
167 |
" [%lf\t%lf\t%lf]\n" |
168 |
" eta = [%lf\t%lf\t%lf]\n" |
169 |
" [%lf\t%lf\t%lf]\n" |
170 |
" [%lf\t%lf\t%lf]\n", |
171 |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
172 |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
173 |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
174 |
eta(0, 0),eta(0, 1),eta(0, 2), |
175 |
eta(1, 0),eta(1, 1),eta(1, 2), |
176 |
eta(2, 0),eta(2, 1),eta(2, 2)); |
177 |
painCave.isFatal = 1; |
178 |
simError(); |
179 |
} else if (offDiagMax > 0.01) { |
180 |
sprintf( painCave.errMsg, |
181 |
"NPrT error: Attempting an off-diagonal Box scaling of more than 1 percent.\n" |
182 |
" Check your tauBarostat, as it is probably too small!\n\n" |
183 |
" scaleMat = [%lf\t%lf\t%lf]\n" |
184 |
" [%lf\t%lf\t%lf]\n" |
185 |
" [%lf\t%lf\t%lf]\n" |
186 |
" eta = [%lf\t%lf\t%lf]\n" |
187 |
" [%lf\t%lf\t%lf]\n" |
188 |
" [%lf\t%lf\t%lf]\n", |
189 |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
190 |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
191 |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
192 |
eta(0, 0),eta(0, 1),eta(0, 2), |
193 |
eta(1, 0),eta(1, 1),eta(1, 2), |
194 |
eta(2, 0),eta(2, 1),eta(2, 2)); |
195 |
painCave.isFatal = 1; |
196 |
simError(); |
197 |
} else { |
198 |
|
199 |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
200 |
hmat = hmat *scaleMat; |
201 |
currentSnapshot_->setHmat(hmat); |
202 |
|
203 |
} |
204 |
} |
205 |
|
206 |
bool NPrT::etaConverged() { |
207 |
int i; |
208 |
double diffEta, sumEta; |
209 |
|
210 |
sumEta = 0; |
211 |
for(i = 0; i < 3; i++) { |
212 |
sumEta += pow(prevEta(i, i) - eta(i, i), 2); |
213 |
} |
214 |
|
215 |
diffEta = sqrt( sumEta / 3.0 ); |
216 |
|
217 |
return ( diffEta <= etaTolerance ); |
218 |
} |
219 |
|
220 |
double NPrT::calcConservedQuantity(){ |
221 |
|
222 |
chi= currentSnapshot_->getChi(); |
223 |
integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); |
224 |
loadEta(); |
225 |
|
226 |
// We need NkBT a lot, so just set it here: This is the RAW number |
227 |
// of integrableObjects, so no subtraction or addition of constraints or |
228 |
// orientational degrees of freedom: |
229 |
NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp; |
230 |
|
231 |
// fkBT is used because the thermostat operates on more degrees of freedom |
232 |
// than the barostat (when there are particles with orientational degrees |
233 |
// of freedom). |
234 |
fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp; |
235 |
|
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; |
243 |
|
244 |
totalEnergy = thermo.getTotalE(); |
245 |
|
246 |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
247 |
|
248 |
thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
249 |
|
250 |
SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
251 |
trEta = tmp.trace(); |
252 |
|
253 |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
254 |
|
255 |
barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
256 |
|
257 |
conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
258 |
barostat_kinetic + barostat_potential; |
259 |
|
260 |
return conservedQuantity; |
261 |
|
262 |
} |
263 |
|
264 |
void NPrT::loadEta() { |
265 |
eta= currentSnapshot_->getEta(); |
266 |
|
267 |
//if (!eta.isDiagonal()) { |
268 |
// sprintf( painCave.errMsg, |
269 |
// "NPrT error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix"); |
270 |
// painCave.isFatal = 1; |
271 |
// simError(); |
272 |
//} |
273 |
} |
274 |
|
275 |
void NPrT::saveEta() { |
276 |
currentSnapshot_->setEta(eta); |
277 |
} |
278 |
|
279 |
} |
280 |
|
281 |
|