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 |
46 |
|
|
47 |
|
namespace oopse { |
48 |
|
|
49 |
< |
NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) { |
49 |
> |
NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) { |
50 |
|
|
51 |
|
Globals* simParams = info_->getSimParams(); |
52 |
|
|
53 |
< |
if (!simParams->getUseInitXSstate()) { |
54 |
< |
Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
55 |
< |
currSnapshot->setChi(0.0); |
56 |
< |
currSnapshot->setIntegralOfChiDt(0.0); |
53 |
> |
if (!simParams->getUseIntialExtendedSystemState()) { |
54 |
> |
Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
55 |
> |
currSnapshot->setChi(0.0); |
56 |
> |
currSnapshot->setIntegralOfChiDt(0.0); |
57 |
|
} |
58 |
|
|
59 |
|
if (!simParams->haveTargetTemp()) { |
60 |
< |
sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n"); |
61 |
< |
painCave.isFatal = 1; |
62 |
< |
painCave.severity = OOPSE_ERROR; |
63 |
< |
simError(); |
60 |
> |
sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n"); |
61 |
> |
painCave.isFatal = 1; |
62 |
> |
painCave.severity = OOPSE_ERROR; |
63 |
> |
simError(); |
64 |
|
} else { |
65 |
< |
targetTemp_ = simParams->getTargetTemp(); |
65 |
> |
targetTemp_ = simParams->getTargetTemp(); |
66 |
|
} |
67 |
|
|
68 |
|
// We must set tauThermostat_. |
69 |
|
|
70 |
|
if (!simParams->haveTauThermostat()) { |
71 |
< |
sprintf(painCave.errMsg, "If you use the constant temperature\n" |
72 |
< |
"\tintegrator, you must set tauThermostat_.\n"); |
71 |
> |
sprintf(painCave.errMsg, "If you use the constant temperature\n" |
72 |
> |
"\tintegrator, you must set tauThermostat_.\n"); |
73 |
|
|
74 |
< |
painCave.severity = OOPSE_ERROR; |
75 |
< |
painCave.isFatal = 1; |
76 |
< |
simError(); |
74 |
> |
painCave.severity = OOPSE_ERROR; |
75 |
> |
painCave.isFatal = 1; |
76 |
> |
simError(); |
77 |
|
} else { |
78 |
< |
tauThermostat_ = simParams->getTauThermostat(); |
78 |
> |
tauThermostat_ = simParams->getTauThermostat(); |
79 |
|
} |
80 |
|
|
81 |
|
update(); |
82 |
< |
} |
82 |
> |
} |
83 |
|
|
84 |
< |
void NVT::doUpdate() { |
84 |
> |
void NVT::doUpdate() { |
85 |
|
oldVel_.resize(info_->getNIntegrableObjects()); |
86 |
|
oldJi_.resize(info_->getNIntegrableObjects()); |
87 |
< |
} |
88 |
< |
void NVT::moveA() { |
87 |
> |
} |
88 |
> |
void NVT::moveA() { |
89 |
|
SimInfo::MoleculeIterator i; |
90 |
|
Molecule::IntegrableObjectIterator j; |
91 |
|
Molecule* mol; |
105 |
|
double instTemp = thermo.getTemperature(); |
106 |
|
|
107 |
|
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
108 |
< |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
109 |
< |
integrableObject = mol->nextIntegrableObject(j)) { |
108 |
> |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
109 |
> |
integrableObject = mol->nextIntegrableObject(j)) { |
110 |
|
|
111 |
|
vel = integrableObject->getVel(); |
112 |
|
pos = integrableObject->getPos(); |
127 |
|
|
128 |
|
if (integrableObject->isDirectional()) { |
129 |
|
|
130 |
< |
//convert the torque to body frame |
131 |
< |
Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
130 |
> |
//convert the torque to body frame |
131 |
> |
Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
132 |
|
|
133 |
< |
// get the angular momentum, and propagate a half step |
133 |
> |
// get the angular momentum, and propagate a half step |
134 |
|
|
135 |
< |
ji = integrableObject->getJ(); |
135 |
> |
ji = integrableObject->getJ(); |
136 |
|
|
137 |
< |
//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); |
138 |
< |
ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji; |
139 |
< |
rotAlgo->rotate(integrableObject, ji, dt); |
137 |
> |
//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); |
138 |
> |
ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji; |
139 |
> |
rotAlgo->rotate(integrableObject, ji, dt); |
140 |
|
|
141 |
< |
integrableObject->setJ(ji); |
141 |
> |
integrableObject->setJ(ji); |
142 |
|
} |
143 |
< |
} |
143 |
> |
} |
144 |
|
|
145 |
|
} |
146 |
|
|
155 |
|
|
156 |
|
currentSnapshot_->setChi(chi); |
157 |
|
currentSnapshot_->setIntegralOfChiDt(integralOfChidt); |
158 |
< |
} |
158 |
> |
} |
159 |
|
|
160 |
< |
void NVT::moveB() { |
160 |
> |
void NVT::moveB() { |
161 |
|
SimInfo::MoleculeIterator i; |
162 |
|
Molecule::IntegrableObjectIterator j; |
163 |
|
Molecule* mol; |
179 |
|
|
180 |
|
index = 0; |
181 |
|
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
182 |
< |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
183 |
< |
integrableObject = mol->nextIntegrableObject(j)) { |
184 |
< |
oldVel_[index] = integrableObject->getVel(); |
185 |
< |
oldJi_[index] = integrableObject->getJ(); |
182 |
> |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
183 |
> |
integrableObject = mol->nextIntegrableObject(j)) { |
184 |
> |
oldVel_[index] = integrableObject->getVel(); |
185 |
> |
oldJi_[index] = integrableObject->getJ(); |
186 |
|
|
187 |
< |
++index; |
188 |
< |
} |
187 |
> |
++index; |
188 |
> |
} |
189 |
|
|
190 |
|
} |
191 |
|
|
192 |
|
// do the iteration: |
193 |
|
|
194 |
|
for(int k = 0; k < maxIterNum_; k++) { |
195 |
< |
index = 0; |
196 |
< |
instTemp = thermo.getTemperature(); |
195 |
> |
index = 0; |
196 |
> |
instTemp = thermo.getTemperature(); |
197 |
|
|
198 |
< |
// evolve chi another half step using the temperature at t + dt/2 |
199 |
< |
|
200 |
< |
prevChi = chi; |
201 |
< |
chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); |
198 |
> |
// evolve chi another half step using the temperature at t + dt/2 |
199 |
|
|
200 |
< |
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
201 |
< |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
205 |
< |
integrableObject = mol->nextIntegrableObject(j)) { |
200 |
> |
prevChi = chi; |
201 |
> |
chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); |
202 |
|
|
203 |
< |
frc = integrableObject->getFrc(); |
204 |
< |
vel = integrableObject->getVel(); |
203 |
> |
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
204 |
> |
for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
205 |
> |
integrableObject = mol->nextIntegrableObject(j)) { |
206 |
|
|
207 |
< |
mass = integrableObject->getMass(); |
207 |
> |
frc = integrableObject->getFrc(); |
208 |
> |
vel = integrableObject->getVel(); |
209 |
|
|
210 |
< |
// velocity half step |
211 |
< |
//for(j = 0; j < 3; j++) |
212 |
< |
// vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi); |
213 |
< |
vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index]; |
210 |
> |
mass = integrableObject->getMass(); |
211 |
> |
|
212 |
> |
// velocity half step |
213 |
> |
//for(j = 0; j < 3; j++) |
214 |
> |
// vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi); |
215 |
> |
vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index]; |
216 |
|
|
217 |
< |
integrableObject->setVel(vel); |
217 |
> |
integrableObject->setVel(vel); |
218 |
|
|
219 |
< |
if (integrableObject->isDirectional()) { |
219 |
> |
if (integrableObject->isDirectional()) { |
220 |
|
|
221 |
< |
// get and convert the torque to body frame |
221 |
> |
// get and convert the torque to body frame |
222 |
|
|
223 |
< |
Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
223 |
> |
Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
224 |
|
|
225 |
< |
//for(j = 0; j < 3; j++) |
226 |
< |
// ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi); |
227 |
< |
ji = oldJi_[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index]; |
225 |
> |
//for(j = 0; j < 3; j++) |
226 |
> |
// ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi); |
227 |
> |
ji = oldJi_[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index]; |
228 |
|
|
229 |
< |
integrableObject->setJ(ji); |
230 |
< |
} |
229 |
> |
integrableObject->setJ(ji); |
230 |
> |
} |
231 |
|
|
232 |
|
|
233 |
< |
++index; |
234 |
< |
} |
235 |
< |
} |
233 |
> |
++index; |
234 |
> |
} |
235 |
> |
} |
236 |
|
|
237 |
|
|
238 |
< |
rattle->constraintB(); |
238 |
> |
rattle->constraintB(); |
239 |
|
|
240 |
< |
if (fabs(prevChi - chi) <= chiTolerance_) |
241 |
< |
break; |
240 |
> |
if (fabs(prevChi - chi) <= chiTolerance_) |
241 |
> |
break; |
242 |
|
|
243 |
|
} |
244 |
|
|
246 |
|
|
247 |
|
currentSnapshot_->setChi(chi); |
248 |
|
currentSnapshot_->setIntegralOfChiDt(integralOfChidt); |
249 |
< |
} |
249 |
> |
} |
250 |
|
|
251 |
+ |
void NVT::resetIntegrator() { |
252 |
+ |
currentSnapshot_->setChi(0.0); |
253 |
+ |
currentSnapshot_->setIntegralOfChiDt(0.0); |
254 |
+ |
} |
255 |
+ |
|
256 |
+ |
double NVT::calcConservedQuantity() { |
257 |
|
|
252 |
– |
double NVT::calcConservedQuantity() { |
253 |
– |
|
258 |
|
double chi = currentSnapshot_->getChi(); |
259 |
|
double integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); |
260 |
|
double conservedQuantity; |
274 |
|
conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; |
275 |
|
|
276 |
|
return conservedQuantity; |
277 |
< |
} |
277 |
> |
} |
278 |
|
|
279 |
|
|
280 |
|
}//end namespace oopse |