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 1969 by gezelter, Wed Feb 26 14:14:50 2014 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, 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
42	–	#include <cmath>
43	–	#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"
50	–
43		#ifdef IS_MPI
44		#include <mpi.h>
45	<	#define TAKE_THIS_TAG_REAL 2
54	<	#endif //is_mpi
45	>	#endif
46
47	<	namespace oopse {
47	>	#include "restraints/ThermoIntegrationForceManager.hpp"
48	>
49	>	namespace OpenMD {
50
51		ThermoIntegrationForceManager::ThermoIntegrationForceManager(SimInfo* info):
52	<	ForceManager(info){
53	<	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
54	<	simParam = info_->getSimParams();
52	>	RestraintForceManager(info){
53	>	currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
54	>	simParam = info_->getSimParams();
55
56	<	if (simParam->haveThermodynamicIntegrationLambda()){
57	<	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
58	<	}
59	<	else{
60	<	tIntLambda_ = 1.0;
61	<	sprintf(painCave.errMsg,
62	<	"ThermoIntegration error: the transformation parameter\n"
63	<	"\t(lambda) was not specified. OOPSE will use a default\n"
64	<	"\tvalue of %f. To set lambda, use the \n"
65	<	"\tthermodynamicIntegrationLambda variable.\n",
66	<	tIntLambda_);
67	<	painCave.isFatal = 0;
68	<	simError();
69	<	}
56	>	if (simParam->haveThermodynamicIntegrationLambda()){
57	>	tIntLambda_ = simParam->getThermodynamicIntegrationLambda();
58	>	}
59	>	else{
60	>	tIntLambda_ = 1.0;
61	>	sprintf(painCave.errMsg,
62	>	"ThermoIntegration error: the transformation parameter\n"
63	>	"\t(lambda) was not specified. OpenMD will use a default\n"
64	>	"\tvalue of %f. To set lambda, use the \n"
65	>	"\tthermodynamicIntegrationLambda variable.\n",
66	>	tIntLambda_);
67	>	painCave.isFatal = 0;
68	>	simError();
69	>	}
70
71	<	if (simParam->haveThermodynamicIntegrationK()){
72	<	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	<
71	>	if (simParam->haveThermodynamicIntegrationK()){
72	>	tIntK_ = simParam->getThermodynamicIntegrationK();
73		}
74	+	else{
75	+	tIntK_ = 1.0;
76	+	sprintf(painCave.errMsg,
77	+	"ThermoIntegration Warning: the tranformation parameter\n"
78	+	"\texponent (k) was not specified. OpenMD will use a default\n"
79	+	"\tvalue of %f. To set k, use the thermodynamicIntegrationK\n"
80	+	"\tvariable.\n",
81	+	tIntK_);
82	+	painCave.isFatal = 0;
83	+	simError();
84	+	}
85	+
86	+	// build the scaling factor used to modulate the forces and torques
87	+	factor_ = pow(tIntLambda_, tIntK_);
88	+	}
89
90		ThermoIntegrationForceManager::~ThermoIntegrationForceManager(){
91		}
92
93	<	void ThermoIntegrationForceManager::calcForces(bool needPotential,
108	<	bool needStress){
93	>	void ThermoIntegrationForceManager::calcForces(){
94		Snapshot* curSnapshot;
95		SimInfo::MoleculeIterator mi;
96		Molecule* mol;
97		Molecule::IntegrableObjectIterator ii;
98	<	StuntDouble* integrableObject;
98	>	StuntDouble* sd;
99		Vector3d frc;
100		Vector3d trq;
101		Mat3x3d tempTau;
102
103		// perform the standard calcForces first
104	<	ForceManager::calcForces(needPotential, needStress);
104	>	ForceManager::calcForces();
105
106		curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
107
108	<	// now scale forces and torques of all the integrableObjects
108	>	// now scale forces and torques of all the sds
109
110		for (mol = info_->beginMolecule(mi); mol != NULL;
111		mol = info_->nextMolecule(mi)) {
112	<	for (integrableObject = mol->beginIntegrableObject(ii);
113	<	integrableObject != NULL;
114	<	integrableObject = mol->nextIntegrableObject(ii)) {
115	<	frc = integrableObject->getFrc();
112	>
113	>	for (sd = mol->beginIntegrableObject(ii); sd != NULL;
114	>	sd = mol->nextIntegrableObject(ii)) {
115	>
116	>	frc = sd->getFrc();
117		frc *= factor_;
118	<	integrableObject->setFrc(frc);
118	>	sd->setFrc(frc);
119
120	<	if (integrableObject->isDirectional()){
121	<	trq = integrableObject->getTrq();
120	>	if (sd->isDirectional()){
121	>	trq = sd->getTrq();
122		trq *= factor_;
123	<	integrableObject->setTrq(trq);
123	>	sd->setTrq(trq);
124		}
125		}
126		}
141	–
142	–	// set vraw to be the unmodulated potential
143	–	lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
144	–	curSnapshot->statData[Stats::VRAW] = lrPot_;
127
128	+	// set rawPotential to be the unmodulated potential
129	+	lrPot_ = curSnapshot->getLongRangePotential();
130	+	curSnapshot->setRawPotential(lrPot_);
131	+
132		// modulate the potential and update the snapshot
133		lrPot_ *= factor_;
134	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
134	>	curSnapshot->setLongRangePotential(lrPot_);
135
136		// scale the pressure tensor
137	<	tempTau = curSnapshot->statData.getTau();
137	>	tempTau = curSnapshot->getStressTensor();
138		tempTau *= factor_;
139	<	curSnapshot->statData.setTau(tempTau);
140	<	#ifndef IS_MPI
141	<	// do the single processor crystal restraint forces for
142	<	// thermodynamic integration
143	<	if (simParam->getUseSolidThermInt()) {
144	<
145	<	lrPot_ += restraint_->Calc_Restraint_Forces();
146	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
147	<
148	<	vHarm_ = restraint_->getVharm();
163	<	curSnapshot->statData[Stats::VHARM] = vHarm_;
139	>	curSnapshot->setStressTensor(tempTau);
140	>
141	>	// now, on to the applied restraining potentials (if needed):
142	>	RealType restPot_local = 0.0;
143	>	RealType vHarm_local = 0.0;
144	>
145	>	if (simParam->getUseRestraints()) {
146	>	// do restraints from RestraintForceManager:
147	>	restPot_local = doRestraints(1.0 - factor_);
148	>	vHarm_local = getUnscaledPotential();
149		}
150	+
151	+	#ifdef IS_MPI
152	+	RealType restPot;
153	+	MPI_Allreduce(&restPot_local, &restPot, 1,
154	+	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
155	+	MPI_Allreduce(&vHarm_local, &vHarm_, 1,
156	+	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
157	+	lrPot_ += restPot;
158		#else
159	<	double tempLRPot = 0.0;
160	<	double tempVHarm = 0.0;
161	<	MPI_Status ierr;
169	<	int nproc;
170	<	MPI_Comm_size(MPI_COMM_WORLD, &nproc);
171	<	vHarm_ = 0.0;
159	>	lrPot_ += restPot_local;
160	>	vHarm_ = vHarm_local;
161	>	#endif
162
163	<	// do the MPI crystal restraint forces for each processor
164	<	if (simParam->getUseSolidThermInt()) {
165	<	tempLRPot = restraint_->Calc_Restraint_Forces();
166	<	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	<
163	>	// give the final values to stats
164	>	curSnapshot->setLongRangePotential(lrPot_);
165	>	curSnapshot->setRestraintPotential(vHarm_);
166	>	}
167		}

Comparing trunk/src/restraints/ThermoIntegrationForceManager.cpp (property svn:keywords):
Revision 1029 by chrisfen, Thu Aug 31 22:34:49 2006 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

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

Diff Legend

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