ViewVC Help

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

root/OpenMD/branches/development/src/integrators/NPT.cpp

Comparing:
trunk/src/integrators/NPT.cpp (file contents), Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/integrators/NPT.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

#	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. Redistributions of source code must retain the above copyright
10	+	*    notice, this list of conditions and the following disclaimer.
11	+	*
12	+	* 2. Redistributions in binary form must reproduce the above copyright
13	+	*    notice, this list of conditions and the following disclaimer in the
14	+	*    documentation and/or other materials provided with the
15	+	*    distribution.
16	+	*
17	+	* This software is provided "AS IS," without a warranty of any
18	+	* kind. All express or implied conditions, representations and
19	+	* warranties, including any implied warranty of merchantability,
20	+	* fitness for a particular purpose or non-infringement, are hereby
21	+	* excluded.  The University of Notre Dame and its licensors shall not
22	+	* be liable for any damages suffered by licensee as a result of
23	+	* using, modifying or distributing the software or its
24	+	* derivatives. In no event will the University of Notre Dame or its
25	+	* licensors be liable for any lost revenue, profit or data, or for
26	+	* direct, indirect, special, consequential, incidental or punitive
27	+	* damages, however caused and regardless of the theory of liability,
28	+	* arising out of the use of or inability to use software, even if the
29	+	* University of Notre Dame has been advised of the possibility of
30	+	* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	+	*/
42	+
43		#include <math.h>
44
3	–	#include "primitives/Atom.hpp"
4	–	#include "primitives/SRI.hpp"
5	–	#include "primitives/AbstractClasses.hpp"
45		#include "brains/SimInfo.hpp"
7	–	#include "UseTheForce/ForceFields.hpp"
46		#include "brains/Thermo.hpp"
47	<	#include "io/ReadWrite.hpp"
48	<	#include "integrators/Integrator.hpp"
47	>	#include "integrators/NPT.hpp"
48	>	#include "math/SquareMatrix3.hpp"
49	>	#include "primitives/Molecule.hpp"
50	>	#include "utils/PhysicalConstants.hpp"
51		#include "utils/simError.h"
52
13	–	#ifdef IS_MPI
14	–	#include "brains/mpiSimulation.hpp"
15	–	#endif
16	–
17	–
53		// Basic isotropic thermostating and barostating via the Melchionna
54		// modification of the Hoover algorithm:
55		//
#	Line 25 | Line 60
60		//
61		//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
62
63	<	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;
63	>	namespace OpenMD {
64
65	<	chiValue = NULL;
66	<	integralOfChidtValue = NULL;
65	>	NPT::NPT(SimInfo* info) :
66	>	VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
67
68	<	chi = 0.0;
69	<	integralOfChidt = 0.0;
70	<	have_tau_thermostat = 0;
71	<	have_tau_barostat = 0;
72	<	have_target_temp = 0;
73	<	have_target_pressure = 0;
74	<	have_chi_tolerance = 0;
75	<	have_eta_tolerance = 0;
76	<	have_pos_iter_tolerance = 0;
68	>	Globals* simParams = info_->getSimParams();
69	>
70	>	if (!simParams->getUseIntialExtendedSystemState()) {
71	>	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
72	>	currSnapshot->setChi(0.0);
73	>	currSnapshot->setIntegralOfChiDt(0.0);
74	>	currSnapshot->setEta(Mat3x3d(0.0));
75	>	}
76	>
77	>	if (!simParams->haveTargetTemp()) {
78	>	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
79	>	painCave.isFatal = 1;
80	>	painCave.severity = OPENMD_ERROR;
81	>	simError();
82	>	} else {
83	>	targetTemp = simParams->getTargetTemp();
84	>	}
85
86	<	// retrieve chi and integralOfChidt from simInfo
87	<	data = info->getProperty(CHIVALUE_ID);
88	<	if(data){
89	<	chiValue = dynamic_cast<DoubleData*>(data);
52	<	}
86	>	// We must set tauThermostat
87	>	if (!simParams->haveTauThermostat()) {
88	>	sprintf(painCave.errMsg, "If you use the constant temperature\n"
89	>	"\tintegrator, you must set tauThermostat.\n");
90
91	<	data = info->getProperty(INTEGRALOFCHIDT_ID);
92	<	if(data){
93	<	integralOfChidtValue = dynamic_cast<DoubleData*>(data);
94	<	}
91	>	painCave.severity = OPENMD_ERROR;
92	>	painCave.isFatal = 1;
93	>	simError();
94	>	} else {
95	>	tauThermostat = simParams->getTauThermostat();
96	>	}
97
98	<	// chi and integralOfChidt should appear by pair
99	<	if(chiValue && integralOfChidtValue){
100	<	chi = chiValue->getData();
62	<	integralOfChidt = integralOfChidtValue->getData();
63	<	}
98	>	if (!simParams->haveTargetPressure()) {
99	>	sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
100	>	"   without a targetPressure!\n");
101
102	<	oldPos = new double[3*integrableObjects.size()];
103	<	oldVel = new double[3*integrableObjects.size()];
104	<	oldJi = new double[3*integrableObjects.size()];
102	>	painCave.isFatal = 1;
103	>	simError();
104	>	} else {
105	>	targetPressure = simParams->getTargetPressure();
106	>	}
107	>
108	>	if (!simParams->haveTauBarostat()) {
109	>	sprintf(painCave.errMsg,
110	>	"If you use the NPT integrator, you must set tauBarostat.\n");
111	>	painCave.severity = OPENMD_ERROR;
112	>	painCave.isFatal = 1;
113	>	simError();
114	>	} else {
115	>	tauBarostat = simParams->getTauBarostat();
116	>	}
117	>
118	>	tt2 = tauThermostat * tauThermostat;
119	>	tb2 = tauBarostat * tauBarostat;
120
121	<	}
121	>	update();
122	>	}
123
124	<	template<typename T> NPT<T>::~NPT() {
125	<	delete[] oldPos;
73	<	delete[] oldVel;
74	<	delete[] oldJi;
75	<	}
124	>	NPT::~NPT() {
125	>	}
126
127	<	template<typename T> void NPT<T>::moveA() {
127	>	void NPT::doUpdate() {
128
129	<	//new version of NPT
130	<	int i, j, k;
131	<	double Tb[3], ji[3];
82	<	double mass;
83	<	double vel[3], pos[3], frc[3];
84	<	double sc[3];
85	<	double COM[3];
129	>	oldPos.resize(info_->getNIntegrableObjects());
130	>	oldVel.resize(info_->getNIntegrableObjects());
131	>	oldJi.resize(info_->getNIntegrableObjects());
132
133	<	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();
133	>	}
134
135	<	tStats->getCOM(COM);
135	>	void NPT::moveA() {
136	>	SimInfo::MoleculeIterator i;
137	>	Molecule::IntegrableObjectIterator  j;
138	>	Molecule* mol;
139	>	StuntDouble* integrableObject;
140	>	Vector3d Tb, ji;
141	>	RealType mass;
142	>	Vector3d vel;
143	>	Vector3d pos;
144	>	Vector3d frc;
145	>	Vector3d sc;
146	>	int index;
147
148	<	//evolve velocity half step
148	>	chi= currentSnapshot_->getChi();
149	>	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
150	>	loadEta();
151	>
152	>	instaTemp =thermo.getTemperature();
153	>	press = thermo.getPressureTensor();
154	>	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
155	>	instaVol =thermo.getVolume();
156
157	<	calcVelScale();
97	<
98	<	for( i=0; i<integrableObjects.size(); i++ ){
157	>	Vector3d  COM = info_->getCom();
158
159	<	integrableObjects[i]->getVel( vel );
101	<	integrableObjects[i]->getFrc( frc );
159	>	//evolve velocity half step
160
161	<	mass = integrableObjects[i]->getMass();
161	>	calcVelScale();
162
163	<	getVelScaleA( sc, vel );
163	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165	>	integrableObject = mol->nextIntegrableObject(j)) {
166	>
167	>	vel = integrableObject->getVel();
168	>	frc = integrableObject->getFrc();
169
170	<	for (j=0; j < 3; j++) {
170	>	mass = integrableObject->getMass();
171
172	<	// velocity half step  (use chi from previous step here):
110	<	vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
172	>	getVelScaleA(sc, vel);
173
174	<	}
174	>	// velocity half step  (use chi from previous step here):
175	>	//vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
176	>	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
177	>	integrableObject->setVel(vel);
178
179	<	integrableObjects[i]->setVel( vel );
179	>	if (integrableObject->isDirectional()) {
180
181	<	if( integrableObjects[i]->isDirectional() ){
181	>	// get and convert the torque to body frame
182
183	<	// get and convert the torque to body frame
183	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184
185	<	integrableObjects[i]->getTrq( Tb );
121	<	integrableObjects[i]->lab2Body( Tb );
185	>	// get the angular momentum, and propagate a half step
186
187	<	// get the angular momentum, and propagate a half step
187	>	ji = integrableObject->getJ();
188
189	<	integrableObjects[i]->getJ( ji );
189	>	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
190	>	ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
191	>
192	>	rotAlgo->rotate(integrableObject, ji, dt);
193
194	<	for (j=0; j < 3; j++)
195	<	ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
194	>	integrableObject->setJ(ji);
195	>	}
196	>
197	>	}
198	>	}
199	>	// evolve chi and eta  half step
200
201	<	this->rotationPropagation( integrableObjects[i], ji );
201	>	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202	>
203	>	evolveEtaA();
204
205	<	integrableObjects[i]->setJ( ji );
205	>	//calculate the integral of chidt
206	>	integralOfChidt += dt2 * chi;
207	>
208	>	index = 0;
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)) {
212	>	oldPos[index++] = integrableObject->getPos();
213	>	}
214		}
215	<	}
215	>
216	>	//the first estimation of r(t+dt) is equal to  r(t)
217
218	<	// evolve chi and eta  half step
218	>	for(int k = 0; k < maxIterNum_; k++) {
219	>	index = 0;
220	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
221	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
222	>	integrableObject = mol->nextIntegrableObject(j)) {
223
224	<	evolveChiA();
225	<	evolveEtaA();
224	>	vel = integrableObject->getVel();
225	>	pos = integrableObject->getPos();
226
227	<	//calculate the integral of chidt
142	<	integralOfChidt += dt2*chi;
227	>	this->getPosScale(pos, COM, index, sc);
228
229	<	//save the old positions
230	<	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	<	}
229	>	pos = oldPos[index] + dt * (vel + sc);
230	>	integrableObject->setPos(pos);
231
232	<	//the first estimation of r(t+dt) is equal to  r(t)
232	>	++index;
233	>	}
234	>	}
235
236	<	for(k = 0; k < 5; k ++){
236	>	rattle->constraintA();
237	>	}
238
239	<	for(i =0 ; i < integrableObjects.size(); i++){
239	>	// Scale the box after all the positions have been moved:
240
241	<	integrableObjects[i]->getVel(vel);
158	<	integrableObjects[i]->getPos(pos);
241	>	this->scaleSimBox();
242
243	<	this->getPosScale( pos, COM, i, sc );
243	>	currentSnapshot_->setChi(chi);
244	>	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
245
246	<	for(j = 0; j < 3; j++)
163	<	pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
164	<
165	<	integrableObjects[i]->setPos( pos );
166	<	}
167	<
168	<	if(nConstrained)
169	<	constrainA();
246	>	saveEta();
247		}
248
249	+	void NPT::moveB(void) {
250	+	SimInfo::MoleculeIterator i;
251	+	Molecule::IntegrableObjectIterator  j;
252	+	Molecule* mol;
253	+	StuntDouble* integrableObject;
254	+	int index;
255	+	Vector3d Tb;
256	+	Vector3d ji;
257	+	Vector3d sc;
258	+	Vector3d vel;
259	+	Vector3d frc;
260	+	RealType mass;
261
173	–	// Scale the box after all the positions have been moved:
262
263	<	this->scaleSimBox();
264	<	}
263	>	chi= currentSnapshot_->getChi();
264	>	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
265	>	RealType oldChi  = chi;
266	>	RealType prevChi;
267
268	<	template<typename T> void NPT<T>::moveB( void ){
269	<
270	<	//new version of NPT
271	<	int i, j, k;
272	<	double Tb[3], ji[3], sc[3];
273	<	double vel[3], frc[3];
274	<	double mass;
275	<
276	<	// Set things up for the iteration:
277	<
278	<	for( i=0; i<integrableObjects.size(); i++ ){
279	<
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	<
268	>	loadEta();
269	>
270	>	//save velocity and angular momentum
271	>	index = 0;
272	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
273	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
274	>	integrableObject = mol->nextIntegrableObject(j)) {
275	>
276	>	oldVel[index] = integrableObject->getVel();
277	>	oldJi[index] = integrableObject->getJ();
278	>	++index;
279	>	}
280		}
203	–	}
281
282	<	// do the iteration:
282	>	// do the iteration:
283	>	instaVol =thermo.getVolume();
284
285	<	instaVol = tStats->getVolume();
285	>	for(int k = 0; k < maxIterNum_; k++) {
286	>	instaTemp =thermo.getTemperature();
287	>	instaPress =thermo.getPressure();
288
289	<	for (k=0; k < 4; k++) {
289	>	// evolve chi another half step using the temperature at t + dt/2
290	>	prevChi = chi;
291	>	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
292
293	<	instaTemp = tStats->getTemperature();
294	<	instaPress = tStats->getPressure();
293	>	//evolve eta
294	>	this->evolveEtaB();
295	>	this->calcVelScale();
296
297	<	// evolve chi another half step using the temperature at t + dt/2
297	>	index = 0;
298	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
299	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
300	>	integrableObject = mol->nextIntegrableObject(j)) {
301
302	<	this->evolveChiB();
303	<	this->evolveEtaB();
218	<	this->calcVelScale();
302	>	frc = integrableObject->getFrc();
303	>	vel = integrableObject->getVel();
304
305	<	for( i=0; i<integrableObjects.size(); i++ ){
305	>	mass = integrableObject->getMass();
306
307	<	integrableObjects[i]->getFrc( frc );
223	<	integrableObjects[i]->getVel(vel);
307	>	getVelScaleB(sc, index);
308
309	<	mass = integrableObjects[i]->getMass();
309	>	// velocity half step
310	>	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
311	>	vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
312	>	integrableObject->setVel(vel);
313
314	<	getVelScaleB( sc, i );
314	>	if (integrableObject->isDirectional()) {
315	>	// get and convert the torque to body frame
316	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
317
318	<	// velocity half step
319	<	for (j=0; j < 3; j++)
320	<	vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
318	>	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
319	>	ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
320	>	integrableObject->setJ(ji);
321	>	}
322
323	<	integrableObjects[i]->setVel( vel );
324	<
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 );
323	>	++index;
324	>	}
325		}
326	+
327	+	rattle->constraintB();
328	+
329	+	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
330	+	break;
331		}
332
333	<	if(nConstrained)
334	<	constrainB();
333	>	//calculate integral of chidt
334	>	integralOfChidt += dt2 * chi;
335
336	<	if ( this->chiConverged() && this->etaConverged() ) break;
337	<	}
336	>	currentSnapshot_->setChi(chi);
337	>	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
338
339	<	//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;
339	>	saveEta();
340		}
341
342	<	if (!have_target_pressure) {
343	<	sprintf( painCave.errMsg,
344	<	"NPT error: You can't use the NPT integrator\n"
345	<	"   without a targetPressure!\n"
305	<	);
306	<	painCave.isFatal = 1;
307	<	simError();
308	<	return -1;
342	>	void NPT::resetIntegrator(){
343	>	currentSnapshot_->setChi(0.0);
344	>	currentSnapshot_->setIntegralOfChiDt(0.0);
345	>	resetEta();
346		}
347
311	–	// We must set tauThermostat.
348
349	<	if (!have_tau_thermostat) {
350	<	sprintf( painCave.errMsg,
351	<	"NPT error: If you use the NPT\n"
352	<	"   integrator, you must set tauThermostat.\n");
353	<	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;
349	>	void NPT::resetEta() {
350	>	Mat3x3d etaMat(0.0);
351	>	currentSnapshot_->setEta(etaMat);
352	>	}
353	>
354		}

Comparing:
trunk/src/integrators/NPT.cpp (property svn:keywords), Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/integrators/NPT.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

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