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 665 by tim, Thu Oct 13 22:26:47 2005 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 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
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	<	namespace oopse {
46	>	#ifdef IS_MPI
47	>	#include <mpi.h>
48	>	#endif
49	>
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 (lambda) was\n"
64	<	"\tnot specified. OOPSE will use a default value of %f. To set\n"
65	<	"\tlambda, use the thermodynamicIntegrationLambda variable.\n",
66	<	tIntLambda_);
67	<	painCave.isFatal = 0;
68	<	simError();
69	<	}
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();
74	<	}
75	<	else{
76	<	tIntK_ = 1.0;
77	<	sprintf(painCave.errMsg,
78	<	"ThermoIntegration Warning: the tranformation parameter exponent\n"
79	<	"\t(k) was not specified. OOPSE will use a default value of %f.\n"
80	<	"\tTo set k, use the thermodynamicIntegrationK variable.\n",
81	<	tIntK_);
82	<	painCave.isFatal = 0;
83	<	simError();
84	<	}
85	<
86	<	if (simParam->getUseSolidThermInt()) {
87	<	// build a restraint object
88	<	restraint_ =  new Restraints(info_, tIntLambda_, tIntK_);
89	<
90	<	}
91	<
92	<	// build the scaling factor used to modulate the forces and torques
93	<	factor_ = pow(tIntLambda_, tIntK_);
94	<
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,
101	<	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		}
134	–
135	–	// set vraw to be the unmodulated potential
136	–	lrPot_ = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL];
137	–	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);
140	>	curSnapshot->setStressTensor(tempTau);
141
142	<	// do crystal restraint forces for thermodynamic integration
143	<	if (simParam->getUseSolidThermInt()) {
144	<
151	<	lrPot_ += restraint_->Calc_Restraint_Forces();
152	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot_;
153	<
154	<	vHarm_ = restraint_->getVharm();
155	<	curSnapshot->statData[Stats::VHARM] = vHarm_;
156	<	}
142	>	// now, on to the applied restraining potentials (if needed):
143	>	RealType restPot_local = 0.0;
144	>	RealType vHarm_local = 0.0;
145
146	<	}
147	<
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	>	lrPot_ += restPot_local;
161	>	vHarm_ = vHarm_local;
162	>	#endif
163	>
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 665 by tim, Thu Oct 13 22:26:47 2005 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