# | Line 1 | Line 1 | |
---|---|---|
1 | < | #include <math.h> |
2 | < | #include "math/MatVec3.h" |
3 | < | #include "primitives/Atom.hpp" |
4 | < | #include "primitives/SRI.hpp" |
5 | < | #include "primitives/AbstractClasses.hpp" |
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" | |
7 | – | #include "UseTheForce/ForceFields.hpp" |
43 | #include "brains/Thermo.hpp" | |
44 | < | #include "io/ReadWrite.hpp" |
45 | < | #include "integrators/Integrator.hpp" |
44 | > | #include "integrators/IntegratorCreator.hpp" |
45 | > | #include "integrators/NPTxyz.hpp" |
46 | > | #include "primitives/Molecule.hpp" |
47 | > | #include "utils/OOPSEConstant.hpp" |
48 | #include "utils/simError.h" | |
49 | ||
13 | – | #ifdef IS_MPI |
14 | – | #include "brains/mpiSimulation.hpp" |
15 | – | #endif |
16 | – | |
50 | // Basic non-isotropic thermostating and barostating via the Melchionna | |
51 | // modification of the Hoover algorithm: | |
52 | // | |
# | Line 24 | Line 57 | |
57 | // | |
58 | // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. | |
59 | ||
60 | < | template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff): |
28 | < | T( theInfo, the_ff ) |
29 | < | { |
30 | < | GenericData* data; |
31 | < | DoubleArrayData * etaValue; |
32 | < | vector<double> etaArray; |
33 | < | int i,j; |
60 | > | namespace oopse { |
61 | ||
62 | < | for(i = 0; i < 3; i++){ |
63 | < | for (j = 0; j < 3; j++){ |
62 | > | |
63 | > | double NPTxyz::calcConservedQuantity(){ |
64 | ||
65 | < | eta[i][j] = 0.0; |
66 | < | oldEta[i][j] = 0.0; |
67 | < | } |
68 | < | } |
65 | > | // We need NkBT a lot, so just set it here: This is the RAW number |
66 | > | // of integrableObjects, so no subtraction or addition of constraints or |
67 | > | // orientational degrees of freedom: |
68 | > | NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp; |
69 | ||
70 | + | // fkBT is used because the thermostat operates on more degrees of freedom |
71 | + | // than the barostat (when there are particles with orientational degrees |
72 | + | // of freedom). |
73 | + | fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp; |
74 | ||
75 | < | if( theInfo->useInitXSstate ){ |
75 | > | double conservedQuantity; |
76 | > | double totalEnergy; |
77 | > | double thermostat_kinetic; |
78 | > | double thermostat_potential; |
79 | > | double barostat_kinetic; |
80 | > | double barostat_potential; |
81 | > | double trEta; |
82 | ||
83 | < | // retrieve eta array from simInfo if it exists |
47 | < | data = info->getProperty(ETAVALUE_ID); |
48 | < | if(data){ |
49 | < | etaValue = dynamic_cast<DoubleArrayData*>(data); |
50 | < | |
51 | < | if(etaValue){ |
52 | < | etaArray = etaValue->getData(); |
53 | < | |
54 | < | for(i = 0; i < 3; i++){ |
55 | < | for (j = 0; j < 3; j++){ |
56 | < | eta[i][j] = etaArray[3*i+j]; |
57 | < | oldEta[i][j] = eta[i][j]; |
58 | < | } |
59 | < | } |
60 | < | } |
61 | < | } |
62 | < | } |
63 | < | } |
83 | > | totalEnergy = thermo.getTotalE(); |
84 | ||
85 | < | template<typename T> NPTxyz<T>::~NPTxyz() { |
85 | > | thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert); |
86 | ||
87 | < | // empty for now |
68 | < | } |
87 | > | thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert; |
88 | ||
89 | < | template<typename T> void NPTxyz<T>::resetIntegrator() { |
89 | > | SquareMatrix<double, 3> tmp = eta.transpose() * eta; |
90 | > | trEta = tmp.trace(); |
91 | ||
92 | < | int i, j; |
92 | > | barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert); |
93 | ||
94 | < | for(i = 0; i < 3; i++) |
75 | < | for (j = 0; j < 3; j++) |
76 | < | eta[i][j] = 0.0; |
94 | > | barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert; |
95 | ||
96 | < | T::resetIntegrator(); |
97 | < | } |
96 | > | conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
97 | > | barostat_kinetic + barostat_potential; |
98 | ||
81 | – | template<typename T> void NPTxyz<T>::evolveEtaA() { |
99 | ||
100 | < | int i, j; |
100 | > | return conservedQuantity; |
101 | ||
85 | – | for(i = 0; i < 3; i ++){ |
86 | – | for(j = 0; j < 3; j++){ |
87 | – | if( i == j) |
88 | – | eta[i][j] += dt2 * instaVol * |
89 | – | (press[i][j] - targetPressure/p_convert) / (NkBT*tb2); |
90 | – | else |
91 | – | eta[i][j] = 0.0; |
92 | – | } |
102 | } | |
103 | ||
104 | < | for(i = 0; i < 3; i++) |
105 | < | for (j = 0; j < 3; j++) |
97 | < | oldEta[i][j] = eta[i][j]; |
98 | < | } |
104 | > | |
105 | > | void NPTxyz::scaleSimBox(){ |
106 | ||
107 | < | template<typename T> void NPTxyz<T>::evolveEtaB() { |
107 | > | int i,j,k; |
108 | > | Mat3x3d scaleMat; |
109 | > | double eta2ij, scaleFactor; |
110 | > | double bigScale, smallScale, offDiagMax; |
111 | > | Mat3x3d hm; |
112 | > | Mat3x3d hmnew; |
113 | ||
102 | – | int i,j; |
114 | ||
104 | – | for(i = 0; i < 3; i++) |
105 | – | for (j = 0; j < 3; j++) |
106 | – | prevEta[i][j] = eta[i][j]; |
115 | ||
116 | < | for(i = 0; i < 3; i ++){ |
109 | < | for(j = 0; j < 3; j++){ |
110 | < | if( i == j) { |
111 | < | eta[i][j] = oldEta[i][j] + dt2 * instaVol * |
112 | < | (press[i][j] - targetPressure/p_convert) / (NkBT*tb2); |
113 | < | } else { |
114 | < | eta[i][j] = 0.0; |
115 | < | } |
116 | < | } |
117 | < | } |
118 | < | } |
116 | > | // Scale the box after all the positions have been moved: |
117 | ||
118 | < | template<typename T> void NPTxyz<T>::calcVelScale(void) { |
119 | < | int i,j; |
118 | > | // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
119 | > | // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
120 | ||
121 | < | for (i = 0; i < 3; i++ ) { |
122 | < | for (j = 0; j < 3; j++ ) { |
123 | < | vScale[i][j] = eta[i][j]; |
121 | > | bigScale = 1.0; |
122 | > | smallScale = 1.0; |
123 | > | offDiagMax = 0.0; |
124 | ||
125 | < | if (i == j) { |
126 | < | vScale[i][j] += chi; |
125 | > | for(i=0; i<3; i++){ |
126 | > | for(j=0; j<3; j++){ |
127 | > | scaleMat(i, j) = 0.0; |
128 | > | if(i==j) { |
129 | > | scaleMat(i, j) = 1.0; |
130 | > | } |
131 | } | |
132 | } | |
131 | – | } |
132 | – | } |
133 | ||
134 | < | template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) { |
135 | < | matVecMul3( vScale, vel, sc ); |
136 | < | } |
134 | > | for(i=0;i<3;i++){ |
135 | ||
136 | < | template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){ |
139 | < | int j; |
140 | < | double myVel[3]; |
136 | > | // calculate the scaleFactors |
137 | ||
138 | < | for (j = 0; j < 3; j++) |
143 | < | myVel[j] = oldVel[3*index + j]; |
138 | > | scaleFactor = exp(dt*eta(i, i)); |
139 | ||
140 | < | matVecMul3( vScale, myVel, sc ); |
146 | < | } |
140 | > | scaleMat(i, i) = scaleFactor; |
141 | ||
142 | < | template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3], |
143 | < | int index, double sc[3]){ |
144 | < | int j; |
145 | < | double rj[3]; |
142 | > | if (scaleMat(i, i) > bigScale) { |
143 | > | bigScale = scaleMat(i, i); |
144 | > | } |
145 | > | |
146 | > | if (scaleMat(i, i) < smallScale) { |
147 | > | smallScale = scaleMat(i, i); |
148 | > | } |
149 | > | } |
150 | ||
151 | < | for(j=0; j<3; j++) |
152 | < | rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j]; |
151 | > | if ((bigScale > 1.1) || (smallScale < 0.9)) { |
152 | > | sprintf( painCave.errMsg, |
153 | > | "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n" |
154 | > | " Check your tauBarostat, as it is probably too small!\n\n" |
155 | > | " scaleMat = [%lf\t%lf\t%lf]\n" |
156 | > | " [%lf\t%lf\t%lf]\n" |
157 | > | " [%lf\t%lf\t%lf]\n", |
158 | > | scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2), |
159 | > | scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2), |
160 | > | scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2)); |
161 | > | painCave.isFatal = 1; |
162 | > | simError(); |
163 | > | } else { |
164 | ||
165 | < | matVecMul3( eta, rj, sc ); |
166 | < | } |
167 | < | |
159 | < | template<typename T> void NPTxyz<T>::scaleSimBox( void ){ |
160 | < | |
161 | < | int i,j,k; |
162 | < | double scaleMat[3][3]; |
163 | < | double eta2ij, scaleFactor; |
164 | < | double bigScale, smallScale, offDiagMax; |
165 | < | double hm[3][3], hmnew[3][3]; |
166 | < | |
167 | < | |
168 | < | |
169 | < | // Scale the box after all the positions have been moved: |
170 | < | |
171 | < | // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
172 | < | // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
173 | < | |
174 | < | bigScale = 1.0; |
175 | < | smallScale = 1.0; |
176 | < | offDiagMax = 0.0; |
177 | < | |
178 | < | for(i=0; i<3; i++){ |
179 | < | for(j=0; j<3; j++){ |
180 | < | scaleMat[i][j] = 0.0; |
181 | < | if(i==j) scaleMat[i][j] = 1.0; |
165 | > | Mat3x3d hmat = currentSnapshot_->getHmat(); |
166 | > | hmat = hmat *scaleMat; |
167 | > | currentSnapshot_->setHmat(hmat); |
168 | } | |
169 | } | |
170 | ||
171 | < | for(i=0;i<3;i++){ |
172 | < | |
187 | < | // calculate the scaleFactors |
188 | < | |
189 | < | scaleFactor = exp(dt*eta[i][i]); |
190 | < | |
191 | < | scaleMat[i][i] = scaleFactor; |
192 | < | |
193 | < | if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i]; |
194 | < | if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i]; |
171 | > | void NPTxyz::loadEta() { |
172 | > | eta= currentSnapshot_->getEta(); |
173 | } | |
174 | ||
197 | – | // for(i=0; i<3; i++){ |
198 | – | // for(j=0; j<3; j++){ |
199 | – | |
200 | – | // // Calculate the matrix Product of the eta array (we only need |
201 | – | // // the ij element right now): |
202 | – | |
203 | – | // eta2ij = 0.0; |
204 | – | // for(k=0; k<3; k++){ |
205 | – | // eta2ij += eta[i][k] * eta[k][j]; |
206 | – | // } |
207 | – | |
208 | – | // scaleMat[i][j] = 0.0; |
209 | – | // // identity matrix (see above): |
210 | – | // if (i == j) scaleMat[i][j] = 1.0; |
211 | – | // // Taylor expansion for the exponential truncated at second order: |
212 | – | // scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij; |
213 | – | |
214 | – | // if (i != j) |
215 | – | // if (fabs(scaleMat[i][j]) > offDiagMax) |
216 | – | // offDiagMax = fabs(scaleMat[i][j]); |
217 | – | // } |
218 | – | |
219 | – | // if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i]; |
220 | – | // if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i]; |
221 | – | // } |
222 | – | |
223 | – | if ((bigScale > 1.1) || (smallScale < 0.9)) { |
224 | – | sprintf( painCave.errMsg, |
225 | – | "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n" |
226 | – | " Check your tauBarostat, as it is probably too small!\n\n" |
227 | – | " scaleMat = [%lf\t%lf\t%lf]\n" |
228 | – | " [%lf\t%lf\t%lf]\n" |
229 | – | " [%lf\t%lf\t%lf]\n", |
230 | – | scaleMat[0][0],scaleMat[0][1],scaleMat[0][2], |
231 | – | scaleMat[1][0],scaleMat[1][1],scaleMat[1][2], |
232 | – | scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]); |
233 | – | painCave.isFatal = 1; |
234 | – | simError(); |
235 | – | } else { |
236 | – | info->getBoxM(hm); |
237 | – | matMul3(hm, scaleMat, hmnew); |
238 | – | info->setBoxM(hmnew); |
239 | – | } |
175 | } | |
241 | – | |
242 | – | template<typename T> bool NPTxyz<T>::etaConverged() { |
243 | – | int i; |
244 | – | double diffEta, sumEta; |
245 | – | |
246 | – | sumEta = 0; |
247 | – | for(i = 0; i < 3; i++) |
248 | – | sumEta += pow(prevEta[i][i] - eta[i][i], 2); |
249 | – | |
250 | – | diffEta = sqrt( sumEta / 3.0 ); |
251 | – | |
252 | – | return ( diffEta <= etaTolerance ); |
253 | – | } |
254 | – | |
255 | – | template<typename T> double NPTxyz<T>::getConservedQuantity(void){ |
256 | – | |
257 | – | double conservedQuantity; |
258 | – | double totalEnergy; |
259 | – | double thermostat_kinetic; |
260 | – | double thermostat_potential; |
261 | – | double barostat_kinetic; |
262 | – | double barostat_potential; |
263 | – | double trEta; |
264 | – | double a[3][3], b[3][3]; |
265 | – | |
266 | – | totalEnergy = tStats->getTotalE(); |
267 | – | |
268 | – | thermostat_kinetic = fkBT * tt2 * chi * chi / |
269 | – | (2.0 * eConvert); |
270 | – | |
271 | – | thermostat_potential = fkBT* integralOfChidt / eConvert; |
272 | – | |
273 | – | transposeMat3(eta, a); |
274 | – | matMul3(a, eta, b); |
275 | – | trEta = matTrace3(b); |
276 | – | |
277 | – | barostat_kinetic = NkBT * tb2 * trEta / |
278 | – | (2.0 * eConvert); |
279 | – | |
280 | – | barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / |
281 | – | eConvert; |
282 | – | |
283 | – | conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential + |
284 | – | barostat_kinetic + barostat_potential; |
285 | – | |
286 | – | // cout.width(8); |
287 | – | // cout.precision(8); |
288 | – | |
289 | – | // cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << |
290 | – | // "\t" << thermostat_potential << "\t" << barostat_kinetic << |
291 | – | // "\t" << barostat_potential << "\t" << conservedQuantity << endl; |
292 | – | |
293 | – | return conservedQuantity; |
294 | – | |
295 | – | } |
296 | – | |
297 | – | template<typename T> string NPTxyz<T>::getAdditionalParameters(void){ |
298 | – | string parameters; |
299 | – | const int BUFFERSIZE = 2000; // size of the read buffer |
300 | – | char buffer[BUFFERSIZE]; |
301 | – | |
302 | – | sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt); |
303 | – | parameters += buffer; |
304 | – | |
305 | – | for(int i = 0; i < 3; i++){ |
306 | – | sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]); |
307 | – | parameters += buffer; |
308 | – | } |
309 | – | |
310 | – | return parameters; |
311 | – | |
312 | – | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |