ViewVC Help

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

root/OpenMD/trunk/src/restraints/ThermoIntegrationForceManager.cpp

Comparing trunk/src/restraints/ThermoIntegrationForceManager.cpp (file contents):
Revision 1029 by chrisfen, Thu Aug 31 22:34:49 2006 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
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.
#	Line 37 | Line 28
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 <cmath>
43	>
44		#include "restraints/ThermoIntegrationForceManager.hpp"
44	–	#include "integrators/Integrator.hpp"
45	–	#include "math/SquareMatrix3.hpp"
46	–	#include "primitives/Molecule.hpp"
47	–	#include "utils/simError.h"
48	–	#include "utils/OOPSEConstant.hpp"
49	–	#include "utils/StringUtils.hpp"
45
46		#ifdef IS_MPI
47		#include <mpi.h>
48	<	#define TAKE_THIS_TAG_REAL 2
54	<	#endif //is_mpi
48	>	#endif
49
50	<	namespace oopse {
50	>	namespace OpenMD {
51
52		ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info):
53	<	ForceManager(info){
54	<	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
55	<	simParam = info_->getSimParams();
53	>	RestraintForceManager(info){
54	>	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
55	>	simParam = info_->getSimParams();
56
57	<	if (simParam->haveThermodynamicIntegrationLambda()){
58	<	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
59	<	}
60	<	else{
61	<	tIntLambda_ = 1.0;
62	<	sprintf(painCave.errMsg,
63	<	"ThermoIntegration error: the transformation parameter\n"
64	<	"\t(lambda) was not specified. OOPSE will use a default\n"
65	<	"\tvalue of %f. To set lambda, use the \n"
66	<	"\tthermodynamicIntegrationLambda variable.\n",
67	<	tIntLambda_);
68	<	painCave.isFatal = 0;
69	<	simError();
70	<	}
57	>	if (simParam->haveThermodynamicIntegrationLambda()){
58	>	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
59	>	}
60	>	else{
61	>	tIntLambda_ = 1.0;
62	>	sprintf(painCave.errMsg,
63	>	"ThermoIntegration error: the transformation parameter\n"
64	>	"\t(lambda) was not specified. OpenMD will use a default\n"
65	>	"\tvalue of %f. To set lambda, use the \n"
66	>	"\tthermodynamicIntegrationLambda variable.\n",
67	>	tIntLambda_);
68	>	painCave.isFatal = 0;
69	>	simError();
70	>	}
71
72	<	if (simParam->haveThermodynamicIntegrationK()){
73	<	tIntK_ = simParam->getThermodynamicIntegrationK();
80	<	}
81	<	else{
82	<	tIntK_ = 1.0;
83	<	sprintf(painCave.errMsg,
84	<	"ThermoIntegration Warning: the tranformation parameter\n"
85	<	"\texponent (k) was not specified. OOPSE will use a default\n"
86	<	"\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
87	<	"\tvariable.\n",
88	<	tIntK_);
89	<	painCave.isFatal = 0;
90	<	simError();
91	<	}
92	<
93	<	if (simParam->getUseSolidThermInt()) {
94	<	// build a restraint object
95	<	restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
96	<
97	<	}
98	<
99	<	// build the scaling factor used to modulate the forces and torques
100	<	factor_ = pow(tIntLambda_, tIntK_);
101	<
72	>	if (simParam->haveThermodynamicIntegrationK()){
73	>	tIntK_ = simParam->getThermodynamicIntegrationK();
74		}
75	+	else{
76	+	tIntK_ = 1.0;
77	+	sprintf(painCave.errMsg,
78	+	"ThermoIntegration Warning: the tranformation parameter\n"
79	+	"\texponent (k) was not specified. OpenMD will use a default\n"
80	+	"\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
81	+	"\tvariable.\n",
82	+	tIntK_);
83	+	painCave.isFatal = 0;
84	+	simError();
85	+	}
86	+
87	+	// build the scaling factor used to modulate the forces and torques
88	+	factor_ = pow(tIntLambda_, tIntK_);
89	+	}
90
91		ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
92		}
93
94	<	void ThermoIntegrationForceManager::calcForces(bool needPotential,
108	<	bool needStress){
94	>	void ThermoIntegrationForceManager::calcForces(){
95		Snapshot* curSnapshot;
96		SimInfo::MoleculeIterator mi;
97		Molecule* mol;
98		Molecule::IntegrableObjectIterator ii;
99	<	StuntDouble* integrableObject;
99	>	StuntDouble* sd;
100		Vector3d frc;
101		Vector3d trq;
102		Mat3x3d tempTau;
103
104		// perform the standard calcForces first
105	<	ForceManager::calcForces(needPotential, needStress);
105	>	ForceManager::calcForces();
106
107		curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
108
109	<	// now scale forces and torques of all the integrableObjects
109	>	// now scale forces and torques of all the sds
110
111		for (mol = info_->beginMolecule(mi); mol != NULL;
112		mol = info_->nextMolecule(mi)) {
113	<	for (integrableObject = mol->beginIntegrableObject(ii);
114	<	integrableObject != NULL;
115	<	integrableObject = mol->nextIntegrableObject(ii)) {
116	<	frc = integrableObject->getFrc();
113	>
114	>	for (sd = mol->beginIntegrableObject(ii); sd != NULL;
115	>	sd = mol->nextIntegrableObject(ii)) {
116	>
117	>	frc = sd->getFrc();
118		frc *= factor_;
119	<	integrableObject->setFrc(frc);
119	>	sd->setFrc(frc);
120
121	<	if (integrableObject->isDirectional()){
122	<	trq = integrableObject->getTrq();
121	>	if (sd->isDirectional()){
122	>	trq = sd->getTrq();
123		trq *= factor_;
124	<	integrableObject->setTrq(trq);
124	>	sd->setTrq(trq);
125		}
126		}
127		}
141	–
142	–	// set vraw to be the unmodulated potential
143	–	lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
144	–	curSnapshot->statData[Stats::VRAW] = lrPot_;
128
129	+	// set rawPotential to be the unmodulated potential
130	+	lrPot_ = curSnapshot->getLongRangePotential();
131	+	curSnapshot->setRawPotential(lrPot_);
132	+
133		// modulate the potential and update the snapshot
134		lrPot_ *= factor_;
135	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
135	>	curSnapshot->setLongRangePotential(lrPot_);
136
137		// scale the pressure tensor
138	<	tempTau = curSnapshot->statData.getTau();
138	>	tempTau = curSnapshot->getStressTensor();
139		tempTau *= factor_;
140	<	curSnapshot->statData.setTau(tempTau);
141	<	#ifndef IS_MPI
142	<	// do the single processor crystal restraint forces for
143	<	// thermodynamic integration
144	<	if (simParam->getUseSolidThermInt()) {
145	<
146	<	lrPot_ += restraint_->Calc_Restraint_Forces();
147	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
148	<
149	<	vHarm_ = restraint_->getVharm();
163	<	curSnapshot->statData[Stats::VHARM] = vHarm_;
140	>	curSnapshot->setStressTensor(tempTau);
141	>
142	>	// now, on to the applied restraining potentials (if needed):
143	>	RealType restPot_local = 0.0;
144	>	RealType vHarm_local = 0.0;
145	>
146	>	if (simParam->getUseRestraints()) {
147	>	// do restraints from RestraintForceManager:
148	>	restPot_local = doRestraints(1.0 - factor_);
149	>	vHarm_local = getUnscaledPotential();
150		}
151	+
152	+	#ifdef IS_MPI
153	+	RealType restPot;
154	+	MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1,
155	+	MPI::REALTYPE, MPI::SUM);
156	+	MPI::COMM_WORLD.Allreduce(&vHarm_local, &vHarm_, 1,
157	+	MPI::REALTYPE, MPI::SUM);
158	+	lrPot_ += restPot;
159		#else
160	<	double tempLRPot = 0.0;
161	<	double tempVHarm = 0.0;
162	<	MPI_Status ierr;
169	<	int nproc;
170	<	MPI_Comm_size(MPI_COMM_WORLD, &nproc);
171	<	vHarm_ = 0.0;
160	>	lrPot_ += restPot_local;
161	>	vHarm_ = vHarm_local;
162	>	#endif
163
164	<	// do the MPI crystal restraint forces for each processor
165	<	if (simParam->getUseSolidThermInt()) {
166	<	tempLRPot = restraint_->Calc_Restraint_Forces();
167	<	tempVHarm = restraint_->getVharm();
177	<	}
178	<
179	<	// master receives and accumulates the restraint info
180	<	if (worldRank == 0) {
181	<	for(int i = 0 ; i < nproc; ++i) {
182	<	if (i == worldRank) {
183	<	lrPot_ += tempLRPot;
184	<	vHarm_ += tempVHarm;
185	<	} else {
186	<	MPI_Recv(&tempLRPot, 1, MPI_REALTYPE, i,
187	<	TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
188	<	MPI_Recv(&tempVHarm, 1, MPI_REALTYPE, i,
189	<	TAKE_THIS_TAG_REAL, MPI_COMM_WORLD, &ierr);
190	<	lrPot_ += tempLRPot;
191	<	vHarm_ += tempVHarm;
192	<	}
193	<	}
194	<
195	<	// give the final values to stats
196	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
197	<	curSnapshot->statData[Stats::VHARM] = vHarm_;
198	<
199	<	} else {
200	<	// pack up and send the appropriate info to the master node
201	<	for(int j = 1; j < nproc; ++j) {
202	<	if (worldRank == j) {
203	<
204	<	MPI_Send(&tempLRPot, 1, MPI_REALTYPE, 0,
205	<	TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
206	<	MPI_Send(&tempVHarm, 1, MPI_REALTYPE, 0,
207	<	TAKE_THIS_TAG_REAL, MPI_COMM_WORLD);
208	<	}
209	<	}
210	<	}
211	<	#endif //is_mpi
212	<	}
213	<
164	>	// give the final values to stats
165	>	curSnapshot->setLongRangePotential(lrPot_);
166	>	curSnapshot->setRestraintPotential(vHarm_);
167	>	}
168		}

Comparing trunk/src/restraints/ThermoIntegrationForceManager.cpp (property svn:keywords):
Revision 1029 by chrisfen, Thu Aug 31 22:34:49 2006 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