ViewVC Help

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

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

Comparing trunk/src/integrators/NPTi.cpp (file contents):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

#	Line 1 | Line 1
1	<	#include <math.h>
2	<	#include "Atom.hpp"
3	<	#include "SRI.hpp"
4	<	#include "AbstractClasses.hpp"
5	<	#include "SimInfo.hpp"
6	<	#include "ForceFields.hpp"
7	<	#include "Thermo.hpp"
8	<	#include "ReadWrite.hpp"
9	<	#include "Integrator.hpp"
10	<	#include "simError.h"
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 "NPTi.hpp"
44	>	#include "brains/SimInfo.hpp"
45	>	#include "brains/Thermo.hpp"
46	>	#include "integrators/NPT.hpp"
47	>	#include "primitives/Molecule.hpp"
48	>	#include "utils/PhysicalConstants.hpp"
49	>	#include "utils/simError.h"
50
51	<	#ifdef IS_MPI
13	<	#include "mpiSimulation.hpp"
14	<	#endif
51	>	namespace OpenMD {
52
53	<	// Basic isotropic thermostating and barostating via the Melchionna
54	<	// modification of the Hoover algorithm:
55	<	//
56	<	//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
57	<	//       Molec. Phys., 78, 533.
58	<	//
59	<	//           and
60	<	//
61	<	//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
53	>	// Basic isotropic thermostating and barostating via the Melchionna
54	>	// modification of the Hoover algorithm:
55	>	//
56	>	//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
57	>	//       Molec. Phys., 78, 533.
58	>	//
59	>	//           and
60	>	//
61	>	//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
62
63	<	template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
27	<	T( theInfo, the_ff )
28	<	{
29	<	GenericData* data;
30	<	DoubleArrayData * etaValue;
31	<	vector<double> etaArray;
63	>	NPTi::NPTi ( SimInfo *info) : NPT(info){
64
65	<	eta = 0.0;
34	<	oldEta = 0.0;
65	>	}
66
67	<	if( theInfo->useInitXSstate ){
68	<	// retrieve eta from simInfo if
69	<	data = info->getProperty(ETAVALUE_ID);
70	<	if(data){
40	<	etaValue = dynamic_cast<DoubleArrayData*>(data);
41	<
42	<	if(etaValue){
43	<	etaArray = etaValue->getData();
44	<	eta = etaArray[0];
45	<	oldEta = eta;
46	<	}
47	<	}
67	>	void NPTi::evolveEtaA() {
68	>	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
69	>	(PhysicalConstants::pressureConvert*NkBT*tb2));
70	>	oldEta = eta;
71		}
49	–	}
72
73	<	template<typename T> NPTi<T>::~NPTi() {
52	<	//nothing for now
53	<	}
73	>	void NPTi::evolveEtaB() {
74
75	<	template<typename T> void NPTi<T>::resetIntegrator() {
76	<	eta = 0.0;
77	<	T::resetIntegrator();
78	<	}
75	>	prevEta = eta;
76	>	eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
77	>	(PhysicalConstants::pressureConvert*NkBT*tb2));
78	>	}
79
80	<	template<typename T> void NPTi<T>::evolveEtaA() {
81	<	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
82	<	(p_convert*NkBT*tb2));
63	<	oldEta = eta;
64	<	}
80	>	void NPTi::calcVelScale() {
81	>	vScale = thermostat.first + eta;
82	>	}
83
84	<	template<typename T> void NPTi<T>::evolveEtaB() {
84	>	void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
85	>	sc = vel * vScale;
86	>	}
87
88	<	prevEta = eta;
89	<	eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
90	<	(p_convert*NkBT*tb2));
71	<	}
88	>	void NPTi::getVelScaleB(Vector3d& sc, int index ){
89	>	sc = oldVel[index] * vScale;
90	>	}
91
73	–	template<typename T> void NPTi<T>::calcVelScale(void) {
74	–	vScale = chi + eta;
75	–	}
92
93	<	template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
94	<	int i;
93	>	void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
94	>	int index, Vector3d& sc){
95	>	/**@todo*/
96	>	sc  = (oldPos[index] + pos)/(RealType)2.0 -COM;
97	>	sc *= eta;
98	>	}
99
100	<	for(i=0; i<3; i++) sc[i] = vel[i] * vScale;
81	<	}
100	>	void NPTi::scaleSimBox(){
101
102	<	template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
84	<	int i;
102	>	RealType scaleFactor;
103
104	<	for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * vScale;
87	<	}
104	>	scaleFactor = exp(dt*eta);
105
106	+	if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
107	+	sprintf( painCave.errMsg,
108	+	"NPTi error: Attempting a Box scaling of more than 10 percent"
109	+	" check your tauBarostat, as it is probably too small!\n"
110	+	" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
111	+	);
112	+	painCave.isFatal = 1;
113	+	simError();
114	+	} else {
115	+	Mat3x3d hmat = snap->getHmat();
116	+	hmat *= scaleFactor;
117	+	snap->setHmat(hmat);
118	+	}
119
120	<	template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
91	<	int index, double sc[3]){
92	<	int j;
120	>	}
121
122	<	for(j=0; j<3; j++)
95	<	sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
122	>	bool NPTi::etaConverged() {
123
124	<	for(j=0; j<3; j++)
98	<	sc[j] *= eta;
99	<	}
100	<
101	<	template<typename T> void NPTi<T>::scaleSimBox( void ){
102	<
103	<	double scaleFactor;
104	<
105	<	scaleFactor = exp(dt*eta);
106	<
107	<	if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
108	<	sprintf( painCave.errMsg,
109	<	"NPTi error: Attempting a Box scaling of more than 10 percent"
110	<	" check your tauBarostat, as it is probably too small!\n"
111	<	" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
112	<	);
113	<	painCave.isFatal = 1;
114	<	simError();
115	<	} else {
116	<	info->scaleBox(scaleFactor);
124	>	return ( fabs(prevEta - eta) <= etaTolerance );
125		}
126
127	<	}
127	>	RealType NPTi::calcConservedQuantity(){
128
129	<	template<typename T> bool NPTi<T>::etaConverged() {
129	>	thermostat = snap->getThermostat();
130	>	loadEta();
131	>	// We need NkBT a lot, so just set it here: This is the RAW number
132	>	// of integrableObjects, so no subtraction or addition of constraints or
133	>	// orientational degrees of freedom:
134	>	NkBT = info_->getNGlobalIntegrableObjects()*PhysicalConstants::kB *targetTemp;
135
136	<	return ( fabs(prevEta - eta) <= etaTolerance );
137	<	}
136	>	// fkBT is used because the thermostat operates on more degrees of freedom
137	>	// than the barostat (when there are particles with orientational degrees
138	>	// of freedom).
139	>	fkBT = info_->getNdf()*PhysicalConstants::kB *targetTemp;
140	>
141	>	RealType conservedQuantity;
142	>	RealType Energy;
143	>	RealType thermostat_kinetic;
144	>	RealType thermostat_potential;
145	>	RealType barostat_kinetic;
146	>	RealType barostat_potential;
147
148	<	template<typename T> double NPTi<T>::getConservedQuantity(void){
148	>	Energy =thermo.getTotalEnergy();
149
150	<	double conservedQuantity;
151	<	double Energy;
130	<	double thermostat_kinetic;
131	<	double thermostat_potential;
132	<	double barostat_kinetic;
133	<	double barostat_potential;
150	>	thermostat_kinetic = fkBT* tt2 * thermostat.first *
151	>	thermostat.first / (2.0 * PhysicalConstants::energyConvert);
152
153	<	Energy = tStats->getTotalE();
153	>	thermostat_potential = fkBT* thermostat.second / PhysicalConstants::energyConvert;
154
137	–	thermostat_kinetic = fkBT* tt2 * chi * chi /
138	–	(2.0 * eConvert);
155
156	<	thermostat_potential = fkBT* integralOfChidt / eConvert;
156	>	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * PhysicalConstants::energyConvert);
157
158	+	barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /
159	+	PhysicalConstants::energyConvert;
160
161	<	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
162	<	(2.0 * eConvert);
161	>	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
162	>	barostat_kinetic + barostat_potential;
163	>
164	>	return conservedQuantity;
165	>	}
166
167	<	barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
168	<	eConvert;
167	>	void NPTi::loadEta() {
168	>	Mat3x3d etaMat = snap->getBarostat();
169	>	eta = etaMat(0,0);
170	>	//if (fabs(etaMat(1,1) - eta) >= OpenMD::epsilon || fabs(etaMat(1,1) - eta) >= OpenMD::epsilon || !etaMat.isDiagonal()) {
171	>	//    sprintf( painCave.errMsg,
172	>	//             "NPTi error: the diagonal elements of  eta matrix are not the same or etaMat is not a diagonal matrix");
173	>	//    painCave.isFatal = 1;
174	>	//    simError();
175	>	//}
176	>	}
177
178	<	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
179	<	barostat_kinetic + barostat_potential;
180	<
181	<	//   cout.width(8);
182	<	//   cout.precision(8);
183	<
184	<	//   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
156	<	//       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
157	<	//       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
158	<	return conservedQuantity;
178	>	void NPTi::saveEta() {
179	>	Mat3x3d etaMat(0.0);
180	>	etaMat(0, 0) = eta;
181	>	etaMat(1, 1) = eta;
182	>	etaMat(2, 2) = eta;
183	>	snap->setBarostat(etaMat);
184	>	}
185		}
160	–
161	–	template<typename T> string NPTi<T>::getAdditionalParameters(void){
162	–	string parameters;
163	–	const int BUFFERSIZE = 2000; // size of the read buffer
164	–	char buffer[BUFFERSIZE];
165	–
166	–	sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
167	–	parameters += buffer;
168	–
169	–	sprintf(buffer,"\t%G\t0\t0;", eta);
170	–	parameters += buffer;
171	–
172	–	sprintf(buffer,"\t0\t%G\t0;", eta);
173	–	parameters += buffer;
174	–
175	–	sprintf(buffer,"\t0\t0\t%G;", eta);
176	–	parameters += buffer;
177	–
178	–	return parameters;
179	–
180	–	}

Comparing trunk/src/integrators/NPTi.cpp (property svn:keywords):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

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

Diff Legend

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