# | Line 1 | Line 1 | |
---|---|---|
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 <math.h> | |
43 | ||
3 | – | #include "primitives/Atom.hpp" |
4 | – | #include "primitives/SRI.hpp" |
5 | – | #include "primitives/AbstractClasses.hpp" |
44 | #include "brains/SimInfo.hpp" | |
7 | – | #include "UseTheForce/ForceFields.hpp" |
45 | #include "brains/Thermo.hpp" | |
46 | < | #include "io/ReadWrite.hpp" |
47 | < | #include "integrators/Integrator.hpp" |
46 | > | #include "integrators/NPT.hpp" |
47 | > | #include "math/SquareMatrix3.hpp" |
48 | > | #include "primitives/Molecule.hpp" |
49 | > | #include "utils/OOPSEConstant.hpp" |
50 | #include "utils/simError.h" | |
51 | ||
13 | – | #ifdef IS_MPI |
14 | – | #include "brains/mpiSimulation.hpp" |
15 | – | #endif |
16 | – | |
17 | – | |
52 | // Basic isotropic thermostating and barostating via the Melchionna | |
53 | // modification of the Hoover algorithm: | |
54 | // | |
# | Line 25 | Line 59 | |
59 | // | |
60 | // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. | |
61 | ||
62 | < | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff): |
29 | < | T( theInfo, the_ff ) |
30 | < | { |
31 | < | GenericData* data; |
32 | < | DoubleData * chiValue; |
33 | < | DoubleData * integralOfChidtValue; |
62 | > | namespace oopse { |
63 | ||
64 | < | chiValue = NULL; |
65 | < | integralOfChidtValue = NULL; |
64 | > | NPT::NPT(SimInfo* info) : |
65 | > | VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) { |
66 | ||
67 | < | chi = 0.0; |
68 | < | integralOfChidt = 0.0; |
69 | < | have_tau_thermostat = 0; |
70 | < | have_tau_barostat = 0; |
71 | < | have_target_temp = 0; |
72 | < | have_target_pressure = 0; |
73 | < | have_chi_tolerance = 0; |
74 | < | have_eta_tolerance = 0; |
75 | < | have_pos_iter_tolerance = 0; |
67 | > | Globals* simParams = info_->getSimParams(); |
68 | > | |
69 | > | if (!simParams->getUseIntialExtendedSystemState()) { |
70 | > | Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
71 | > | currSnapshot->setChi(0.0); |
72 | > | currSnapshot->setIntegralOfChiDt(0.0); |
73 | > | currSnapshot->setEta(Mat3x3d(0.0)); |
74 | > | } |
75 | > | |
76 | > | if (!simParams->haveTargetTemp()) { |
77 | > | sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n"); |
78 | > | painCave.isFatal = 1; |
79 | > | painCave.severity = OOPSE_ERROR; |
80 | > | simError(); |
81 | > | } else { |
82 | > | targetTemp = simParams->getTargetTemp(); |
83 | > | } |
84 | ||
85 | < | // retrieve chi and integralOfChidt from simInfo |
86 | < | data = info->getProperty(CHIVALUE_ID); |
87 | < | if(data){ |
88 | < | chiValue = dynamic_cast<DoubleData*>(data); |
52 | < | } |
85 | > | // We must set tauThermostat |
86 | > | if (!simParams->haveTauThermostat()) { |
87 | > | sprintf(painCave.errMsg, "If you use the constant temperature\n" |
88 | > | "\tintegrator, you must set tauThermostat_.\n"); |
89 | ||
90 | < | data = info->getProperty(INTEGRALOFCHIDT_ID); |
91 | < | if(data){ |
92 | < | integralOfChidtValue = dynamic_cast<DoubleData*>(data); |
93 | < | } |
90 | > | painCave.severity = OOPSE_ERROR; |
91 | > | painCave.isFatal = 1; |
92 | > | simError(); |
93 | > | } else { |
94 | > | tauThermostat = simParams->getTauThermostat(); |
95 | > | } |
96 | ||
97 | < | // chi and integralOfChidt should appear by pair |
98 | < | if(chiValue && integralOfChidtValue){ |
99 | < | chi = chiValue->getData(); |
62 | < | integralOfChidt = integralOfChidtValue->getData(); |
63 | < | } |
97 | > | if (!simParams->haveTargetPressure()) { |
98 | > | sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n" |
99 | > | " without a targetPressure!\n"); |
100 | ||
101 | < | oldPos = new double[3*integrableObjects.size()]; |
102 | < | oldVel = new double[3*integrableObjects.size()]; |
103 | < | oldJi = new double[3*integrableObjects.size()]; |
101 | > | painCave.isFatal = 1; |
102 | > | simError(); |
103 | > | } else { |
104 | > | targetPressure = simParams->getTargetPressure(); |
105 | > | } |
106 | > | |
107 | > | if (!simParams->haveTauBarostat()) { |
108 | > | sprintf(painCave.errMsg, |
109 | > | "If you use the NPT integrator, you must set tauBarostat.\n"); |
110 | > | painCave.severity = OOPSE_ERROR; |
111 | > | painCave.isFatal = 1; |
112 | > | simError(); |
113 | > | } else { |
114 | > | tauBarostat = simParams->getTauBarostat(); |
115 | > | } |
116 | > | |
117 | > | tt2 = tauThermostat * tauThermostat; |
118 | > | tb2 = tauBarostat * tauBarostat; |
119 | ||
120 | < | } |
120 | > | update(); |
121 | > | } |
122 | ||
123 | < | template<typename T> NPT<T>::~NPT() { |
124 | < | delete[] oldPos; |
73 | < | delete[] oldVel; |
74 | < | delete[] oldJi; |
75 | < | } |
123 | > | NPT::~NPT() { |
124 | > | } |
125 | ||
126 | < | template<typename T> void NPT<T>::moveA() { |
126 | > | void NPT::doUpdate() { |
127 | ||
128 | < | //new version of NPT |
129 | < | int i, j, k; |
130 | < | double Tb[3], ji[3]; |
82 | < | double mass; |
83 | < | double vel[3], pos[3], frc[3]; |
84 | < | double sc[3]; |
85 | < | double COM[3]; |
128 | > | oldPos.resize(info_->getNIntegrableObjects()); |
129 | > | oldVel.resize(info_->getNIntegrableObjects()); |
130 | > | oldJi.resize(info_->getNIntegrableObjects()); |
131 | ||
132 | < | instaTemp = tStats->getTemperature(); |
88 | < | tStats->getPressureTensor( press ); |
89 | < | instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; |
90 | < | instaVol = tStats->getVolume(); |
132 | > | } |
133 | ||
134 | < | tStats->getCOM(COM); |
134 | > | void NPT::moveA() { |
135 | > | SimInfo::MoleculeIterator i; |
136 | > | Molecule::IntegrableObjectIterator j; |
137 | > | Molecule* mol; |
138 | > | StuntDouble* integrableObject; |
139 | > | Vector3d Tb, ji; |
140 | > | double mass; |
141 | > | Vector3d vel; |
142 | > | Vector3d pos; |
143 | > | Vector3d frc; |
144 | > | Vector3d sc; |
145 | > | int index; |
146 | ||
147 | < | //evolve velocity half step |
147 | > | chi= currentSnapshot_->getChi(); |
148 | > | integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); |
149 | > | loadEta(); |
150 | > | |
151 | > | instaTemp =thermo.getTemperature(); |
152 | > | press = thermo.getPressureTensor(); |
153 | > | instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0; |
154 | > | instaVol =thermo.getVolume(); |
155 | ||
156 | < | calcVelScale(); |
97 | < | |
98 | < | for( i=0; i<integrableObjects.size(); i++ ){ |
156 | > | Vector3d COM = info_->getCom(); |
157 | ||
158 | < | integrableObjects[i]->getVel( vel ); |
101 | < | integrableObjects[i]->getFrc( frc ); |
158 | > | //evolve velocity half step |
159 | ||
160 | < | mass = integrableObjects[i]->getMass(); |
160 | > | calcVelScale(); |
161 | ||
162 | < | getVelScaleA( sc, vel ); |
162 | > | for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
163 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
164 | > | integrableObject = mol->nextIntegrableObject(j)) { |
165 | > | |
166 | > | vel = integrableObject->getVel(); |
167 | > | frc = integrableObject->getFrc(); |
168 | ||
169 | < | for (j=0; j < 3; j++) { |
169 | > | mass = integrableObject->getMass(); |
170 | ||
171 | < | // velocity half step (use chi from previous step here): |
110 | < | vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]); |
171 | > | getVelScaleA(sc, vel); |
172 | ||
173 | < | } |
173 | > | // velocity half step (use chi from previous step here): |
174 | > | //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]); |
175 | > | vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc; |
176 | > | integrableObject->setVel(vel); |
177 | ||
178 | < | integrableObjects[i]->setVel( vel ); |
178 | > | if (integrableObject->isDirectional()) { |
179 | ||
180 | < | if( integrableObjects[i]->isDirectional() ){ |
180 | > | // get and convert the torque to body frame |
181 | ||
182 | < | // get and convert the torque to body frame |
182 | > | Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
183 | ||
184 | < | integrableObjects[i]->getTrq( Tb ); |
121 | < | integrableObjects[i]->lab2Body( Tb ); |
184 | > | // get the angular momentum, and propagate a half step |
185 | ||
186 | < | // get the angular momentum, and propagate a half step |
186 | > | ji = integrableObject->getJ(); |
187 | ||
188 | < | integrableObjects[i]->getJ( ji ); |
188 | > | //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); |
189 | > | ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji; |
190 | > | |
191 | > | rotAlgo->rotate(integrableObject, ji, dt); |
192 | ||
193 | < | for (j=0; j < 3; j++) |
194 | < | ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); |
193 | > | integrableObject->setJ(ji); |
194 | > | } |
195 | > | |
196 | > | } |
197 | > | } |
198 | > | // evolve chi and eta half step |
199 | ||
200 | < | this->rotationPropagation( integrableObjects[i], ji ); |
200 | > | chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2; |
201 | > | |
202 | > | evolveEtaA(); |
203 | ||
204 | < | integrableObjects[i]->setJ( ji ); |
204 | > | //calculate the integral of chidt |
205 | > | integralOfChidt += dt2 * chi; |
206 | > | |
207 | > | index = 0; |
208 | > | for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
209 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
210 | > | integrableObject = mol->nextIntegrableObject(j)) { |
211 | > | oldPos[index++] = integrableObject->getPos(); |
212 | > | } |
213 | } | |
214 | < | } |
214 | > | |
215 | > | //the first estimation of r(t+dt) is equal to r(t) |
216 | ||
217 | < | // evolve chi and eta half step |
217 | > | for(int k = 0; k < maxIterNum_; k++) { |
218 | > | index = 0; |
219 | > | for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
220 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
221 | > | integrableObject = mol->nextIntegrableObject(j)) { |
222 | ||
223 | < | evolveChiA(); |
224 | < | evolveEtaA(); |
223 | > | vel = integrableObject->getVel(); |
224 | > | pos = integrableObject->getPos(); |
225 | ||
226 | < | //calculate the integral of chidt |
142 | < | integralOfChidt += dt2*chi; |
226 | > | this->getPosScale(pos, COM, index, sc); |
227 | ||
228 | < | //save the old positions |
229 | < | for(i = 0; i < integrableObjects.size(); i++){ |
146 | < | integrableObjects[i]->getPos(pos); |
147 | < | for(j = 0; j < 3; j++) |
148 | < | oldPos[i*3 + j] = pos[j]; |
149 | < | } |
228 | > | pos = oldPos[index] + dt * (vel + sc); |
229 | > | integrableObject->setPos(pos); |
230 | ||
231 | < | //the first estimation of r(t+dt) is equal to r(t) |
231 | > | ++index; |
232 | > | } |
233 | > | } |
234 | ||
235 | < | for(k = 0; k < 5; k ++){ |
235 | > | rattle->constraintA(); |
236 | > | } |
237 | ||
238 | < | for(i =0 ; i < integrableObjects.size(); i++){ |
238 | > | // Scale the box after all the positions have been moved: |
239 | ||
240 | < | integrableObjects[i]->getVel(vel); |
158 | < | integrableObjects[i]->getPos(pos); |
240 | > | this->scaleSimBox(); |
241 | ||
242 | < | this->getPosScale( pos, COM, i, sc ); |
243 | < | |
162 | < | for(j = 0; j < 3; j++) |
163 | < | pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]); |
242 | > | currentSnapshot_->setChi(chi); |
243 | > | currentSnapshot_->setIntegralOfChiDt(integralOfChidt); |
244 | ||
245 | < | integrableObjects[i]->setPos( pos ); |
166 | < | } |
167 | < | |
168 | < | if(nConstrained) |
169 | < | constrainA(); |
245 | > | saveEta(); |
246 | } | |
247 | ||
248 | + | void NPT::moveB(void) { |
249 | + | SimInfo::MoleculeIterator i; |
250 | + | Molecule::IntegrableObjectIterator j; |
251 | + | Molecule* mol; |
252 | + | StuntDouble* integrableObject; |
253 | + | int index; |
254 | + | Vector3d Tb; |
255 | + | Vector3d ji; |
256 | + | Vector3d sc; |
257 | + | Vector3d vel; |
258 | + | Vector3d frc; |
259 | + | double mass; |
260 | ||
173 | – | // Scale the box after all the positions have been moved: |
261 | ||
262 | < | this->scaleSimBox(); |
263 | < | } |
262 | > | chi= currentSnapshot_->getChi(); |
263 | > | integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); |
264 | > | double oldChi = chi; |
265 | > | double prevChi; |
266 | ||
267 | < | template<typename T> void NPT<T>::moveB( void ){ |
268 | < | |
269 | < | //new version of NPT |
270 | < | int i, j, k; |
271 | < | double Tb[3], ji[3], sc[3]; |
272 | < | double vel[3], frc[3]; |
273 | < | double mass; |
274 | < | |
275 | < | // Set things up for the iteration: |
276 | < | |
277 | < | for( i=0; i<integrableObjects.size(); i++ ){ |
278 | < | |
190 | < | integrableObjects[i]->getVel( vel ); |
191 | < | |
192 | < | for (j=0; j < 3; j++) |
193 | < | oldVel[3*i + j] = vel[j]; |
194 | < | |
195 | < | if( integrableObjects[i]->isDirectional() ){ |
196 | < | |
197 | < | integrableObjects[i]->getJ( ji ); |
198 | < | |
199 | < | for (j=0; j < 3; j++) |
200 | < | oldJi[3*i + j] = ji[j]; |
201 | < | |
267 | > | loadEta(); |
268 | > | |
269 | > | //save velocity and angular momentum |
270 | > | index = 0; |
271 | > | for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
272 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
273 | > | integrableObject = mol->nextIntegrableObject(j)) { |
274 | > | |
275 | > | oldVel[index] = integrableObject->getVel(); |
276 | > | oldJi[index] = integrableObject->getJ(); |
277 | > | ++index; |
278 | > | } |
279 | } | |
203 | – | } |
280 | ||
281 | < | // do the iteration: |
281 | > | // do the iteration: |
282 | > | instaVol =thermo.getVolume(); |
283 | ||
284 | < | instaVol = tStats->getVolume(); |
284 | > | for(int k = 0; k < maxIterNum_; k++) { |
285 | > | instaTemp =thermo.getTemperature(); |
286 | > | instaPress =thermo.getPressure(); |
287 | ||
288 | < | for (k=0; k < 4; k++) { |
288 | > | // evolve chi another half step using the temperature at t + dt/2 |
289 | > | prevChi = chi; |
290 | > | chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2; |
291 | ||
292 | < | instaTemp = tStats->getTemperature(); |
293 | < | instaPress = tStats->getPressure(); |
292 | > | //evolve eta |
293 | > | this->evolveEtaB(); |
294 | > | this->calcVelScale(); |
295 | ||
296 | < | // evolve chi another half step using the temperature at t + dt/2 |
296 | > | index = 0; |
297 | > | for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
298 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
299 | > | integrableObject = mol->nextIntegrableObject(j)) { |
300 | ||
301 | < | this->evolveChiB(); |
302 | < | this->evolveEtaB(); |
218 | < | this->calcVelScale(); |
301 | > | frc = integrableObject->getFrc(); |
302 | > | vel = integrableObject->getVel(); |
303 | ||
304 | < | for( i=0; i<integrableObjects.size(); i++ ){ |
304 | > | mass = integrableObject->getMass(); |
305 | ||
306 | < | integrableObjects[i]->getFrc( frc ); |
223 | < | integrableObjects[i]->getVel(vel); |
306 | > | getVelScaleB(sc, index); |
307 | ||
308 | < | mass = integrableObjects[i]->getMass(); |
308 | > | // velocity half step |
309 | > | //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]); |
310 | > | vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc; |
311 | > | integrableObject->setVel(vel); |
312 | ||
313 | < | getVelScaleB( sc, i ); |
313 | > | if (integrableObject->isDirectional()) { |
314 | > | // get and convert the torque to body frame |
315 | > | Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
316 | ||
317 | < | // velocity half step |
318 | < | for (j=0; j < 3; j++) |
319 | < | vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]); |
317 | > | //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi); |
318 | > | ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index]; |
319 | > | integrableObject->setJ(ji); |
320 | > | } |
321 | ||
322 | < | integrableObjects[i]->setVel( vel ); |
323 | < | |
235 | < | if( integrableObjects[i]->isDirectional() ){ |
236 | < | |
237 | < | // get and convert the torque to body frame |
238 | < | |
239 | < | integrableObjects[i]->getTrq( Tb ); |
240 | < | integrableObjects[i]->lab2Body( Tb ); |
241 | < | |
242 | < | for (j=0; j < 3; j++) |
243 | < | ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); |
244 | < | |
245 | < | integrableObjects[i]->setJ( ji ); |
322 | > | ++index; |
323 | > | } |
324 | } | |
325 | + | |
326 | + | rattle->constraintB(); |
327 | + | |
328 | + | if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged()) |
329 | + | break; |
330 | } | |
331 | ||
332 | < | if(nConstrained) |
333 | < | constrainB(); |
332 | > | //calculate integral of chidt |
333 | > | integralOfChidt += dt2 * chi; |
334 | ||
335 | < | if ( this->chiConverged() && this->etaConverged() ) break; |
336 | < | } |
335 | > | currentSnapshot_->setChi(chi); |
336 | > | currentSnapshot_->setIntegralOfChiDt(integralOfChidt); |
337 | ||
338 | < | //calculate integral of chida |
256 | < | integralOfChidt += dt2*chi; |
257 | < | |
258 | < | |
259 | < | } |
260 | < | |
261 | < | template<typename T> void NPT<T>::resetIntegrator() { |
262 | < | chi = 0.0; |
263 | < | T::resetIntegrator(); |
264 | < | } |
265 | < | |
266 | < | template<typename T> void NPT<T>::evolveChiA() { |
267 | < | chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
268 | < | oldChi = chi; |
269 | < | } |
270 | < | |
271 | < | template<typename T> void NPT<T>::evolveChiB() { |
272 | < | |
273 | < | prevChi = chi; |
274 | < | chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
275 | < | } |
276 | < | |
277 | < | template<typename T> bool NPT<T>::chiConverged() { |
278 | < | |
279 | < | return ( fabs( prevChi - chi ) <= chiTolerance ); |
280 | < | } |
281 | < | |
282 | < | template<typename T> int NPT<T>::readyCheck() { |
283 | < | |
284 | < | //check parent's readyCheck() first |
285 | < | if (T::readyCheck() == -1) |
286 | < | return -1; |
287 | < | |
288 | < | // First check to see if we have a target temperature. |
289 | < | // Not having one is fatal. |
290 | < | |
291 | < | if (!have_target_temp) { |
292 | < | sprintf( painCave.errMsg, |
293 | < | "NPT error: You can't use the NPT integrator\n" |
294 | < | " without a targetTemp!\n" |
295 | < | ); |
296 | < | painCave.isFatal = 1; |
297 | < | simError(); |
298 | < | return -1; |
338 | > | saveEta(); |
339 | } | |
340 | ||
341 | < | if (!have_target_pressure) { |
342 | < | sprintf( painCave.errMsg, |
343 | < | "NPT error: You can't use the NPT integrator\n" |
344 | < | " without a targetPressure!\n" |
305 | < | ); |
306 | < | painCave.isFatal = 1; |
307 | < | simError(); |
308 | < | return -1; |
341 | > | void NPT::resetIntegrator(){ |
342 | > | currentSnapshot_->setChi(0.0); |
343 | > | currentSnapshot_->setIntegralOfChiDt(0.0); |
344 | > | resetEta(); |
345 | } | |
346 | ||
311 | – | // We must set tauThermostat. |
347 | ||
348 | < | if (!have_tau_thermostat) { |
349 | < | sprintf( painCave.errMsg, |
350 | < | "NPT error: If you use the NPT\n" |
351 | < | " integrator, you must set tauThermostat.\n"); |
352 | < | painCave.isFatal = 1; |
318 | < | simError(); |
319 | < | return -1; |
320 | < | } |
321 | < | |
322 | < | // We must set tauBarostat. |
323 | < | |
324 | < | if (!have_tau_barostat) { |
325 | < | sprintf( painCave.errMsg, |
326 | < | "If you use the NPT integrator, you must set tauBarostat.\n"); |
327 | < | painCave.severity = OOPSE_ERROR; |
328 | < | painCave.isFatal = 1; |
329 | < | simError(); |
330 | < | return -1; |
331 | < | } |
332 | < | |
333 | < | if (!have_chi_tolerance) { |
334 | < | sprintf( painCave.errMsg, |
335 | < | "Setting chi tolerance to 1e-6 in NPT integrator\n"); |
336 | < | chiTolerance = 1e-6; |
337 | < | have_chi_tolerance = 1; |
338 | < | painCave.severity = OOPSE_INFO; |
339 | < | painCave.isFatal = 0; |
340 | < | simError(); |
341 | < | } |
342 | < | |
343 | < | if (!have_eta_tolerance) { |
344 | < | sprintf( painCave.errMsg, |
345 | < | "Setting eta tolerance to 1e-6 in NPT integrator"); |
346 | < | etaTolerance = 1e-6; |
347 | < | have_eta_tolerance = 1; |
348 | < | painCave.severity = OOPSE_INFO; |
349 | < | painCave.isFatal = 0; |
350 | < | simError(); |
351 | < | } |
352 | < | |
353 | < | // We need NkBT a lot, so just set it here: This is the RAW number |
354 | < | // of integrableObjects, so no subtraction or addition of constraints or |
355 | < | // orientational degrees of freedom: |
356 | < | |
357 | < | NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp; |
358 | < | |
359 | < | // fkBT is used because the thermostat operates on more degrees of freedom |
360 | < | // than the barostat (when there are particles with orientational degrees |
361 | < | // of freedom). |
362 | < | |
363 | < | fkBT = (double)(info->getNDF()) * kB * targetTemp; |
364 | < | |
365 | < | tt2 = tauThermostat * tauThermostat; |
366 | < | tb2 = tauBarostat * tauBarostat; |
367 | < | |
368 | < | return 1; |
348 | > | void NPT::resetEta() { |
349 | > | Mat3x3d etaMat(0.0); |
350 | > | currentSnapshot_->setEta(etaMat); |
351 | > | } |
352 | > | |
353 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |