ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/integrators/NVT.cpp

Comparing trunk/src/integrators/NVT.cpp (file contents):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

#	Line 35 | Line 35
35		*
36		* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37		* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	<	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	<	* [4]  Vardeman & Gezelter, in progress (2009).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	>	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include "integrators/NVT.hpp"
#	Line 51 | Line 52 | namespace OpenMD {
52		Globals* simParams = info_->getSimParams();
53
54		if (!simParams->getUseIntialExtendedSystemState()) {
55	<	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
56	<	currSnapshot->setChi(0.0);
56	<	currSnapshot->setIntegralOfChiDt(0.0);
55	>	Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
56	>	snap->setThermostat(make_pair(0.0, 0.0));
57		}
58
59		if (!simParams->haveTargetTemp()) {
#	Line 78 | Line 78 | namespace OpenMD {
78		tauThermostat_ = simParams->getTauThermostat();
79		}
80
81	<	update();
81	>	updateSizes();
82		}
83
84	<	void NVT::doUpdate() {
84	>	void NVT::doUpdateSizes() {
85		oldVel_.resize(info_->getNIntegrableObjects());
86	<	oldJi_.resize(info_->getNIntegrableObjects());
86	>	oldJi_.resize(info_->getNIntegrableObjects());
87		}
88	+
89		void NVT::moveA() {
90		SimInfo::MoleculeIterator i;
91		Molecule::IntegrableObjectIterator  j;
92		Molecule* mol;
93	<	StuntDouble* integrableObject;
93	>	StuntDouble* sd;
94		Vector3d Tb;
95		Vector3d ji;
96		RealType mass;
#	Line 97 | Line 98 | namespace OpenMD {
98		Vector3d pos;
99		Vector3d frc;
100
101	<	RealType chi = currentSnapshot_->getChi();
102	<	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
102	<
101	>	pair<RealType, RealType> thermostat = snap->getThermostat();
102	>
103		// We need the temperature at time = t for the chi update below:
104
105		RealType 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)) {
107	>	for (mol = info_->beginMolecule(i); mol != NULL;
108	>	mol = info_->nextMolecule(i)) {
109
110	<	vel = integrableObject->getVel();
111	<	pos = integrableObject->getPos();
113	<	frc = integrableObject->getFrc();
110	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
111	>	sd = mol->nextIntegrableObject(j)) {
112
113	<	mass = integrableObject->getMass();
113	>	vel = sd->getVel();
114	>	pos = sd->getPos();
115	>	frc = sd->getFrc();
116
117	<	// velocity half step  (use chi from previous step here):
118	<	//vel[j] += dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - vel[j]*chi);
119	<	vel += dt2 *PhysicalConstants::energyConvert/mass*frc - dt2*chi*vel;
117	>	mass = sd->getMass();
118	>
119	>	// velocity half step (use chi from previous step here):
120	>	vel += dt2 *PhysicalConstants::energyConvert/mass*frc
121	>	- dt2*thermostat.first*vel;
122
123		// position whole step
122	–	//pos[j] += dt * vel[j];
124		pos += dt * vel;
125
126	<	integrableObject->setVel(vel);
127	<	integrableObject->setPos(pos);
126	>	sd->setVel(vel);
127	>	sd->setPos(pos);
128
129	<	if (integrableObject->isDirectional()) {
129	>	if (sd->isDirectional()) {
130
131		//convert the torque to body frame
132	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
132	>	Tb = sd->lab2Body(sd->getTrq());
133
134		// get the angular momentum, and propagate a half step
135
136	<	ji = integrableObject->getJ();
136	>	ji = sd->getJ();
137
138	<	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
139	<	ji += dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *ji;
139	<	rotAlgo->rotate(integrableObject, ji, dt);
138	>	ji += dt2*PhysicalConstants::energyConvert*Tb
139	>	- dt2*thermostat.first *ji;
140
141	<	integrableObject->setJ(ji);
141	>	rotAlgo_->rotate(sd, ji, dt);
142	>
143	>	sd->setJ(ji);
144		}
145		}
146
147		}
148
149	<	rattle->constraintA();
149	>	flucQ_->moveA();
150	>	rattle_->constraintA();
151
152		// Finally, evolve chi a half step (just like a velocity) using
153		// temperature at time t, not time t+dt/2
154
155	<
156	<	chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
157	<	integralOfChidt += chi * dt2;
155	>	thermostat.first += dt2 * (instTemp / targetTemp_ - 1.0)
156	>	/ (tauThermostat_ * tauThermostat_);
157	>	thermostat.second += thermostat.first * dt2;
158
159	<	currentSnapshot_->setChi(chi);
157	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
159	>	snap->setThermostat(thermostat);
160		}
161
162		void NVT::moveB() {
163		SimInfo::MoleculeIterator i;
164		Molecule::IntegrableObjectIterator  j;
165		Molecule* mol;
166	<	StuntDouble* integrableObject;
166	>	StuntDouble* sd;
167
168		Vector3d Tb;
169		Vector3d ji;
#	Line 172 | Line 174 | namespace OpenMD {
174		int index;
175		// Set things up for the iteration:
176
177	<	RealType chi = currentSnapshot_->getChi();
178	<	RealType oldChi = chi;
177	>	pair<RealType, RealType> thermostat = snap->getThermostat();
178	>	RealType oldChi = thermostat.first;
179		RealType  prevChi;
178	–	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
180
181		index = 0;
182	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
183	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
183	<	integrableObject = mol->nextIntegrableObject(j)) {
184	<	oldVel_[index] = integrableObject->getVel();
185	<	oldJi_[index] = integrableObject->getJ();
182	>	for (mol = info_->beginMolecule(i); mol != NULL;
183	>	mol = info_->nextMolecule(i)) {
184
185	+	for (sd = mol->beginIntegrableObject(j); sd != NULL;
186	+	sd = mol->nextIntegrableObject(j)) {
187	+
188	+	oldVel_[index] = sd->getVel();
189	+
190	+	if (sd->isDirectional())
191	+	oldJi_[index] = sd->getJ();
192	+
193		++index;
194	<	}
189	<
194	>	}
195		}
196
197		// do the iteration:
#	Line 197 | Line 202 | namespace OpenMD {
202
203		// evolve chi another half step using the temperature at t + dt/2
204
205	<	prevChi = chi;
206	<	chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
205	>	prevChi = thermostat.first;
206	>	thermostat.first = oldChi + dt2 * (instTemp / targetTemp_ - 1.0)
207	>	/ (tauThermostat_ * tauThermostat_);
208
209	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
210	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
211	<	integrableObject = mol->nextIntegrableObject(j)) {
209	>	for (mol = info_->beginMolecule(i); mol != NULL;
210	>	mol = info_->nextMolecule(i)) {
211	>
212	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
213	>	sd = mol->nextIntegrableObject(j)) {
214
215	<	frc = integrableObject->getFrc();
216	<	vel = integrableObject->getVel();
215	>	frc = sd->getFrc();
216	>	mass = sd->getMass();
217
210	–	mass = integrableObject->getMass();
211	–
218		// velocity half step
219	<	//for(j = 0; j < 3; j++)
220	<	//    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - oldVel_[3*i + j]*chi);
221	<	vel = oldVel_[index] + dt2/mass*PhysicalConstants::energyConvert * frc - dt2*chi*oldVel_[index];
219	>
220	>	vel = oldVel_[index]
221	>	+ dt2/mass*PhysicalConstants::energyConvert * frc
222	>	- dt2*thermostat.first*oldVel_[index];
223
224	<	integrableObject->setVel(vel);
224	>	sd->setVel(vel);
225
226	<	if (integrableObject->isDirectional()) {
226	>	if (sd->isDirectional()) {
227
228		// get and convert the torque to body frame
229
230	<	Tb =  integrableObject->lab2Body(integrableObject->getTrq());
230	>	Tb =  sd->lab2Body(sd->getTrq());
231
232	<	//for(j = 0; j < 3; j++)
233	<	//    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi_[3*i+j]*chi);
227	<	ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *oldJi_[index];
232	>	ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb
233	>	- dt2*thermostat.first *oldJi_[index];
234
235	<	integrableObject->setJ(ji);
235	>	sd->setJ(ji);
236		}
237
238
#	Line 234 | Line 240 | namespace OpenMD {
240		}
241		}
242
243	+	rattle_->constraintB();
244
245	<	rattle->constraintB();
239	<
240	<	if (fabs(prevChi - chi) <= chiTolerance_)
245	>	if (fabs(prevChi - thermostat.first) <= chiTolerance_)
246		break;
247
248		}
249
250	<	integralOfChidt += dt2 * chi;
250	>	flucQ_->moveB();
251
252	<	currentSnapshot_->setChi(chi);
253	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
252	>	thermostat.second += dt2 * thermostat.first;
253	>	snap->setThermostat(thermostat);
254		}
255
256		void NVT::resetIntegrator() {
257	<	currentSnapshot_->setChi(0.0);
253	<	currentSnapshot_->setIntegralOfChiDt(0.0);
257	>	snap->setThermostat(make_pair(0.0, 0.0));
258		}
259
260		RealType NVT::calcConservedQuantity() {
261
262	<	RealType chi = currentSnapshot_->getChi();
259	<	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
262	>	pair<RealType, RealType> thermostat = snap->getThermostat();
263		RealType conservedQuantity;
264		RealType fkBT;
265		RealType Energy;
#	Line 265 | Line 268 | namespace OpenMD {
268
269		fkBT = info_->getNdf() *PhysicalConstants::kB *targetTemp_;
270
271	<	Energy = thermo.getTotalE();
271	>	Energy = thermo.getTotalEnergy();
272
273	<	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * PhysicalConstants::energyConvert);
273	>	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * thermostat.first * thermostat.first / (2.0 * PhysicalConstants::energyConvert);
274
275	<	thermostat_potential = fkBT * integralOfChidt / PhysicalConstants::energyConvert;
275	>	thermostat_potential = fkBT * thermostat.second / PhysicalConstants::energyConvert;
276
277		conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
278

Comparing trunk/src/integrators/NVT.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines