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/NPAT.hpp" |
46 |
#include "primitives/Molecule.hpp" |
47 |
#include "utils/OOPSEConstant.hpp" |
48 |
#include "utils/simError.h" |
49 |
|
50 |
namespace oopse { |
51 |
|
52 |
void NPAT::evolveEtaA() { |
53 |
|
54 |
eta(2,2) += dt2 * instaVol * (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
55 |
oldEta = eta; |
56 |
} |
57 |
|
58 |
void NPAT::evolveEtaB() { |
59 |
|
60 |
prevEta = eta; |
61 |
eta(2,2) = oldEta(2, 2) + dt2 * instaVol * |
62 |
(press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2); |
63 |
} |
64 |
|
65 |
void NPAT::calcVelScale(){ |
66 |
|
67 |
for (int i = 0; i < 3; i++ ) { |
68 |
for (int j = 0; j < 3; j++ ) { |
69 |
vScale(i, j) = eta(i, j); |
70 |
|
71 |
if (i == j) { |
72 |
vScale(i, j) += chi; |
73 |
} |
74 |
} |
75 |
} |
76 |
} |
77 |
|
78 |
void NPAT::getVelScaleA(Vector3d& sc, const Vector3d& vel){ |
79 |
sc = vScale * vel; |
80 |
} |
81 |
|
82 |
void NPAT::getVelScaleB(Vector3d& sc, int index ) { |
83 |
sc = vScale * oldVel[index]; |
84 |
} |
85 |
|
86 |
void NPAT::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) { |
87 |
|
88 |
/**@todo */ |
89 |
Vector3d rj = (oldPos[index] + pos)/2.0 -COM; |
90 |
sc = eta * rj; |
91 |
} |
92 |
|
93 |
void NPAT::scaleSimBox(){ |
94 |
|
95 |
int i; |
96 |
int j; |
97 |
int k; |
98 |
Mat3x3d scaleMat; |
99 |
double eta2ij; |
100 |
double bigScale, smallScale, offDiagMax; |
101 |
Mat3x3d hm; |
102 |
Mat3x3d hmnew; |
103 |
|
104 |
|
105 |
|
106 |
// Scale the box after all the positions have been moved: |
107 |
|
108 |
// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
109 |
// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
110 |
|
111 |
bigScale = 1.0; |
112 |
smallScale = 1.0; |
113 |
offDiagMax = 0.0; |
114 |
|
115 |
for(i=0; i<3; i++){ |
116 |
for(j=0; j<3; j++){ |
117 |
|
118 |
// Calculate the matrix Product of the eta array (we only need |
119 |
// the ij element right now): |
120 |
|
121 |
eta2ij = 0.0; |
122 |
for(k=0; k<3; k++){ |
123 |
eta2ij += eta(i, k) * eta(k, j); |
124 |
} |
125 |
|
126 |
scaleMat(i, j) = 0.0; |
127 |
// identity matrix (see above): |
128 |
if (i == j) scaleMat(i, j) = 1.0; |
129 |
// Taylor expansion for the exponential truncated at second order: |
130 |
scaleMat(i, j) += dt*eta(i, j) + 0.5*dt*dt*eta2ij; |
131 |
|
132 |
|
133 |
if (i != j) |
134 |
if (fabs(scaleMat(i, j)) > offDiagMax) |
135 |
offDiagMax = fabs(scaleMat(i, j)); |
136 |
} |
137 |
|
138 |
if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i); |
139 |
if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i); |
140 |
} |
141 |
|
142 |
if ((bigScale > 1.01) || (smallScale < 0.99)) { |
143 |
sprintf( painCave.errMsg, |
144 |
"NPAT error: Attempting a Box scaling of more than 1 percent.\n" |
145 |
" Check your tauBarostat, as it is probably too small!\n\n" |
146 |
" scaleMat = [%lf\t%lf\t%lf]\n" |
147 |
" [%lf\t%lf\t%lf]\n" |
148 |
" [%lf\t%lf\t%lf]\n" |
149 |
" eta = [%lf\t%lf\t%lf]\n" |
150 |
" [%lf\t%lf\t%lf]\n" |
151 |
" [%lf\t%lf\t%lf]\n", |
152 |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
153 |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
154 |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
155 |
eta(0, 0),eta(0, 1),eta(0, 2), |
156 |
eta(1, 0),eta(1, 1),eta(1, 2), |
157 |
eta(2, 0),eta(2, 1),eta(2, 2)); |
158 |
painCave.isFatal = 1; |
159 |
simError(); |
160 |
} else if (offDiagMax > 0.01) { |
161 |
sprintf( painCave.errMsg, |
162 |
"NPAT error: Attempting an off-diagonal Box scaling of more than 1 percent.\n" |
163 |
" Check your tauBarostat, as it is probably too small!\n\n" |
164 |
" scaleMat = [%lf\t%lf\t%lf]\n" |
165 |
" [%lf\t%lf\t%lf]\n" |
166 |
" [%lf\t%lf\t%lf]\n" |
167 |
" eta = [%lf\t%lf\t%lf]\n" |
168 |
" [%lf\t%lf\t%lf]\n" |
169 |
" [%lf\t%lf\t%lf]\n", |
170 |
scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
171 |
scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
172 |
scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2), |
173 |
eta(0, 0),eta(0, 1),eta(0, 2), |
174 |
eta(1, 0),eta(1, 1),eta(1, 2), |
175 |
eta(2, 0),eta(2, 1),eta(2, 2)); |
176 |
painCave.isFatal = 1; |
177 |
simError(); |
178 |
} else { |
179 |
|
180 |
Mat3x3d hmat = currentSnapshot_->getHmat(); |
181 |
hmat = hmat *scaleMat; |
182 |
currentSnapshot_->setHmat(hmat); |
183 |
|
184 |
} |
185 |
} |
186 |
|
187 |
bool NPAT::etaConverged() { |
188 |
int i; |
189 |
double diffEta, sumEta; |
190 |
|
191 |
sumEta = 0; |
192 |
for(i = 0; i < 3; i++) { |
193 |
sumEta += pow(prevEta(i, i) - eta(i, i), 2); |
194 |
} |
195 |
|
196 |
diffEta = sqrt( sumEta / 3.0 ); |
197 |
|
198 |
return ( diffEta <= etaTolerance ); |
199 |
} |
200 |
|
201 |
double NPAT::calcConservedQuantity(){ |
202 |
|
203 |
chi= currentSnapshot_->getChi(); |
204 |
integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); |
205 |
loadEta(); |
206 |
|
207 |
// We need NkBT a lot, so just set it here: This is the RAW number |
208 |
// of integrableObjects, so no subtraction or addition of constraints or |
209 |
// orientational degrees of freedom: |
210 |
NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp; |
211 |
|
212 |
// fkBT is used because the thermostat operates on more degrees of freedom |
213 |
// than the barostat (when there are particles with orientational degrees |
214 |
// of freedom). |
215 |
fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp; |
216 |
|
217 |
double conservedQuantity; |
218 |
double totalEnergy; |
219 |
double thermostat_kinetic; |
220 |
double thermostat_potential; |
221 |
double barostat_kinetic; |
222 |
double barostat_potential; |
223 |
double trEta; |
224 |
|
225 |
totalEnergy = thermo.getTotalE(); |
226 |
|
227 |
thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
228 |
|
229 |
thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
230 |
|
231 |
SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
232 |
trEta = tmp.trace(); |
233 |
|
234 |
barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
235 |
|
236 |
barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
237 |
|
238 |
conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
239 |
barostat_kinetic + barostat_potential; |
240 |
|
241 |
return conservedQuantity; |
242 |
|
243 |
} |
244 |
|
245 |
void NPAT::loadEta() { |
246 |
eta= currentSnapshot_->getEta(); |
247 |
|
248 |
//if (!eta.isDiagonal()) { |
249 |
// sprintf( painCave.errMsg, |
250 |
// "NPAT error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix"); |
251 |
// painCave.isFatal = 1; |
252 |
// simError(); |
253 |
//} |
254 |
} |
255 |
|
256 |
void NPAT::saveEta() { |
257 |
currentSnapshot_->setEta(eta); |
258 |
} |
259 |
|
260 |
} |
261 |
|