ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/NPT.cpp
Revision: 1867
Committed: Tue Dec 7 23:08:14 2004 UTC (19 years, 8 months ago) by tim
File size: 9151 byte(s)
Log Message:
NPT in progress

File Contents

# Content
1 #include <math.h>
2
3 #include "brains/SimInfo.hpp"
4 #include "brains/Thermo.hpp"
5 #include "integrators/NPT.hpp"
6 #include "math/SquareMatrix3.hpp"
7 #include "primitives/Molecule.hpp"
8 #include "utils/OOPSEConstant.hpp"
9 #include "utils/simError.h"
10
11 // Basic isotropic thermostating and barostating via the Melchionna
12 // modification of the Hoover algorithm:
13 //
14 // Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
15 // Molec. Phys., 78, 533.
16 //
17 // and
18 //
19 // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
20
21 namespace oopse {
22
23 NPT::NPT(SimInfo* info) :
24 VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6) {
25
26 Globals* simParams = info_->getSimParams();
27
28 if (!simParams->getUseInitXSstate()) {
29 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
30 currSnapshot->setChi(0.0);
31 currSnapshot->setIntegralOfChiDt(0.0);
32 currSnapshot->setEta(Mat3x3d(0.0));
33 }
34
35 if (!simParams->haveTargetTemp()) {
36 sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
37 painCave.isFatal = 1;
38 painCave.severity = OOPSE_ERROR;
39 simError();
40 } else {
41 targetTemp = simParams->getTargetTemp();
42 }
43
44 // We must set tauThermostat
45 if (!simParams->haveTauThermostat()) {
46 sprintf(painCave.errMsg, "If you use the constant temperature\n"
47 "\tintegrator, you must set tauThermostat_.\n");
48
49 painCave.severity = OOPSE_ERROR;
50 painCave.isFatal = 1;
51 simError();
52 } else {
53 tauThermostat = simParams->getTauThermostat();
54 }
55
56 if (!simParams->haveTargetPressure()) {
57 sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
58 " without a targetPressure!\n");
59
60 painCave.isFatal = 1;
61 simError();
62 } else {
63 targetPressure = simParams->getTargetPressure();
64 }
65
66 if (!simParams->haveTauBarostat()) {
67 sprintf(painCave.errMsg,
68 "If you use the NPT integrator, you must set tauBarostat.\n");
69 painCave.severity = OOPSE_ERROR;
70 painCave.isFatal = 1;
71 simError();
72 } else {
73 tauBarostat = simParams->getTauBarostat();
74 }
75
76 tt2 = tauThermostat * tauThermostat;
77 tb2 = tauBarostat * tauBarostat;
78
79 update();
80 }
81
82 NPT::~NPT() {
83 }
84
85 void NPT::doUpdate() {
86
87 oldPos.resize(info_->getNIntegrableObjects());
88 oldVel.resize(info_->getNIntegrableObjects());
89 oldJi.resize(info_->getNIntegrableObjects());
90
91 }
92
93 void NPT::moveA() {
94 SimInfo::MoleculeIterator i;
95 Molecule::IntegrableObjectIterator j;
96 Molecule* mol;
97 StuntDouble* integrableObject;
98 Vector3d Tb, ji;
99 double mass;
100 Vector3d vel;
101 Vector3d pos;
102 Vector3d frc;
103 Vector3d sc;
104 int index;
105
106 chi= currentSnapshot_->getChi();
107 integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
108 loadEta();
109
110 instaTemp =thermo.getTemperature();
111 press = thermo.getPressureTensor();
112 instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
113 instaVol =thermo.getVolume();
114
115 Vector3d COM = info_->getCom();
116
117 //evolve velocity half step
118
119 calcVelScale();
120
121 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
122 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
123 integrableObject = mol->nextIntegrableObject(j)) {
124
125 vel = integrableObject->getVel();
126 frc = integrableObject->getFrc();
127
128 mass = integrableObject->getMass();
129
130 getVelScaleA(sc, vel);
131
132 // velocity half step (use chi from previous step here):
133 //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
134 vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
135 integrableObject->setVel(vel);
136
137 if (integrableObject->isDirectional()) {
138
139 // get and convert the torque to body frame
140
141 Tb = integrableObject->getTrq();
142 integrableObject->lab2Body(Tb);
143
144 // get the angular momentum, and propagate a half step
145
146 ji = integrableObject->getJ();
147
148 //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
149 ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
150
151 rotAlgo->rotate(integrableObject, ji, dt);
152
153 integrableObject->setJ(ji);
154 }
155
156 }
157 }
158 // evolve chi and eta half step
159
160 chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
161
162 evolveEtaA();
163
164 //calculate the integral of chidt
165 integralOfChidt += dt2 * chi;
166
167 index = 0;
168 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
169 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
170 integrableObject = mol->nextIntegrableObject(j)) {
171 oldPos[index++] = integrableObject->getPos();
172 }
173 }
174
175 //the first estimation of r(t+dt) is equal to r(t)
176
177 for(int k = 0; k < maxIterNum_; k++) {
178 index = 0;
179 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
180 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
181 integrableObject = mol->nextIntegrableObject(j)) {
182
183 vel = integrableObject->getVel();
184 pos = integrableObject->getPos();
185
186 this->getPosScale(pos, COM, index, sc);
187
188 pos = oldPos[index] + dt * (vel + sc);
189 integrableObject->setPos(pos);
190
191 ++index;
192 }
193 }
194
195 //constraintAlgorithm->doConstrainA();
196 }
197
198 // Scale the box after all the positions have been moved:
199
200 this->scaleSimBox();
201
202 currentSnapshot_->setChi(chi);
203 currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
204
205 saveEta();
206 }
207
208 void NPT::moveB(void) {
209 SimInfo::MoleculeIterator i;
210 Molecule::IntegrableObjectIterator j;
211 Molecule* mol;
212 StuntDouble* integrableObject;
213 int index;
214 Vector3d Tb;
215 Vector3d ji;
216 Vector3d sc;
217 Vector3d vel;
218 Vector3d frc;
219 double mass;
220
221
222 chi= currentSnapshot_->getChi();
223 integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
224 double oldChi = chi;
225 double prevChi;
226
227 loadEta();
228
229 //save velocity and angular momentum
230 index = 0;
231 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
232 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
233 integrableObject = mol->nextIntegrableObject(j)) {
234
235 oldVel[index] = integrableObject->getVel();
236 oldJi[index] = integrableObject->getJ();
237 ++index;
238 }
239 }
240
241 // do the iteration:
242 instaVol =thermo.getVolume();
243
244 for(int k = 0; k < maxIterNum_; k++) {
245 instaTemp =thermo.getTemperature();
246 instaPress =thermo.getPressure();
247
248 // evolve chi another half step using the temperature at t + dt/2
249 prevChi = chi;
250 chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
251
252 //evolve eta
253 this->evolveEtaB();
254 this->calcVelScale();
255
256 index = 0;
257 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
258 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
259 integrableObject = mol->nextIntegrableObject(j)) {
260
261 frc = integrableObject->getFrc();
262 vel = integrableObject->getVel();
263
264 mass = integrableObject->getMass();
265
266 getVelScaleB(sc, index);
267
268 // velocity half step
269 //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
270 vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
271 integrableObject->setVel(vel);
272
273 if (integrableObject->isDirectional()) {
274 // get and convert the torque to body frame
275 Tb = integrableObject->getTrq();
276 integrableObject->lab2Body(Tb);
277
278 //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
279 ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
280 integrableObject->setJ(ji);
281 }
282
283 ++index;
284 }
285 }
286
287 //constraintAlgorithm->doConstrainB();
288
289 if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
290 break;
291 }
292
293 //calculate integral of chidt
294 integralOfChidt += dt2 * chi;
295
296 currentSnapshot_->setChi(chi);
297 currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
298
299 saveEta();
300 }
301
302 }