ViewVC Help

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

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

Comparing trunk/src/restraints/RestraintForceManager.cpp (file contents):
Revision 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC vs.
Revision 1464 by gezelter, Fri Jul 9 19:29:05 2010 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]  Vardeman & Gezelter, in progress (2009).
40		*/
41
42		#include <cmath>
#	Line 45 | Line 45
45		#include "restraints/ObjectRestraint.hpp"
46		#include "io/RestReader.hpp"
47		#include "utils/simError.h"
48	<	#include "utils/OOPSEConstant.hpp"
48	>	#include "utils/PhysicalConstants.hpp"
49		#include "utils/StringUtils.hpp"
50		#include "selection/SelectionEvaluator.hpp"
51		#include "selection/SelectionManager.hpp"
#	Line 54 | Line 54
54		#endif
55
56
57	<	namespace oopse {
57	>	namespace OpenMD {
58
59		RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {
60
#	Line 78 | Line 78 | namespace oopse {
78		restTime_ = simParam->getStatusTime();
79		} else {
80		sprintf(painCave.errMsg,
81	<	"Restraint warning: If you use restraints without setting\n",
82	<	"\tstatusTime, no restraint data will be written to the rest\n",
81	>	"Restraint warning: If you use restraints without setting\n"
82	>	"\tstatusTime, no restraint data will be written to the rest\n"
83		"\tfile.\n");
84		painCave.isFatal = 0;
85		simError();
#	Line 89 | Line 89 | namespace oopse {
89		int nRestraintStamps = simParam->getNRestraintStamps();
90		std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();
91
92	+	std::vector<int> stuntDoubleIndex;
93	+
94		for (int i = 0; i < nRestraintStamps; i++){
95
96		std::string myType = toUpperCopy(stamp[i]->getType());
#	Line 109 | Line 111 | namespace oopse {
111		molIndex = stamp[i]->getMolIndex();
112		}
113
114	<	Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
113	<
114	<	if (mol == NULL) {
114	>	if (molIndex < 0) {
115		sprintf(painCave.errMsg,
116	<	"Restraint Error: A molecular restraint was specified, but\n"
117	<	"\tno molecule was found with global index %d.\n",
118	<	molIndex);
116	>	"Restraint Error: A molecular restraint was specified\n"
117	>	"\twith a molIndex that was less than 0\n");
118		painCave.isFatal = 1;
119		simError();
120		}
121	+	if (molIndex >= info_->getNGlobalMolecules()) {
122	+	sprintf(painCave.errMsg,
123	+	"Restraint Error: A molecular restraint was specified with\n"
124	+	"\ta molIndex that was greater than the total number of molecules\n");
125	+	painCave.isFatal = 1;
126	+	simError();
127	+	}
128	+
129	+	Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
130
131	+	if (mol == NULL) {
132	+	#ifdef IS_MPI
133	+	// getMoleculeByGlobalIndex returns a NULL in parallel if
134	+	// this proc doesn't have the molecule.  Do a quick check to
135	+	// make sure another processor is supposed to have it.
136	+
137	+	int myrank = MPI::COMM_WORLD.Get_rank();
138	+	if (info_->getMolToProc(molIndex) == myrank) {
139	+
140	+	// If we were supposed to have it but got a null, then freak out.
141	+	#endif
142	+
143	+	sprintf(painCave.errMsg,
144	+	"Restraint Error: A molecular restraint was specified, but\n"
145	+	"\tno molecule was found with global index %d.\n",
146	+	molIndex);
147	+	painCave.isFatal = 1;
148	+	simError();
149	+
150	+	#ifdef IS_MPI
151	+	}
152	+	#endif
153	+	}
154	+
155	+
156	+	#ifdef IS_MPI
157	+	// only handle this molecular restraint if this processor owns the
158	+	// molecule
159	+	int myrank = MPI::COMM_WORLD.Get_rank();
160	+	if (info_->getMolToProc(molIndex) == myrank) {
161	+
162	+	#endif
163	+
164		MolecularRestraint* rest = new MolecularRestraint();
165
166		std::string restPre("mol_");
#	Line 148 | Line 189 | namespace oopse {
189		if (stamp[i]->haveRestrainedSwingXAngle()) {
190		rest->setRestrainedSwingXAngle(stamp[i]->getRestrainedSwingXAngle() * M_PI/180.0);
191		}
192	+	if (stamp[i]->havePrint()) {
193	+	rest->setPrintRestraint(stamp[i]->getPrint());
194	+	}
195
196		restraints_.push_back(rest);
197		mol->addProperty(new RestraintData("Restraint", rest));
198		restrainedMols_.push_back(mol);
199	<
199	>	#ifdef IS_MPI
200	>	}
201	>	#endif
202		} else if (myType.compare("OBJECT") == 0) {
203
204		std::string objectSelection;
#	Line 170 | Line 216 | namespace oopse {
216
217		SelectionEvaluator evaluator(info);
218		SelectionManager seleMan(info);
219	<
219	>
220		evaluator.loadScriptString(objectSelection);
221		seleMan.setSelectionSet(evaluator.evaluate());
222		int selectionCount = seleMan.getSelectionCount();
#	Line 188 | Line 234 | namespace oopse {
234
235		for (sd = seleMan.beginSelected(selei); sd != NULL;
236		sd = seleMan.nextSelected(selei)) {
237	<
237	>	stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());
238	>
239		ObjectRestraint* rest = new ObjectRestraint();
240
241		if (stamp[i]->haveDisplacementSpringConstant()) {
#	Line 211 | Line 258 | namespace oopse {
258		}
259		if (stamp[i]->haveRestrainedSwingYAngle()) {
260		rest->setRestrainedSwingYAngle(stamp[i]->getRestrainedSwingYAngle());
261	<	}
261	>	}
262	>	if (stamp[i]->havePrint()) {
263	>	rest->setPrintRestraint(stamp[i]->getPrint());
264	>	}
265	>
266		restraints_.push_back(rest);
267		sd->addProperty(new RestraintData("Restraint", rest));
268		restrainedObjs_.push_back(sd);
#	Line 224 | Line 275 | namespace oopse {
275		// are times when it won't use restraints at all, so only open the
276		// restraint file if we are actually using restraints:
277
278	<	if (simParam->getUseRestraints()) {
278	>	if (simParam->getUseRestraints()) {
279		std::string refFile = simParam->getRestraint_file();
280	<	RestReader* rr = new RestReader(info, refFile);
230	<
280	>	RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex);
281		rr->readReferenceStructure();
282		}
283
284		restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
285		restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
236	–
286		if(!restOut){
287		sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
288		painCave.isFatal = 1;
#	Line 256 | Line 305 | namespace oopse {
305		currRestTime_ = currSnapshot_->getTime();
306		}
307
308	<	void RestraintForceManager::calcForces(bool needPotential, bool needStress){
308	>	void RestraintForceManager::calcForces(){
309
310	<	ForceManager::calcForces(needPotential, needStress);
310	>	ForceManager::calcForces();
311		RealType restPot_local, restPot;
312
313		restPot_local = doRestraints(1.0);
#	Line 311 | Line 360 | namespace oopse {
360		if (mRest == NULL) {
361		sprintf( painCave.errMsg,
362		"Can not cast RestraintData to MolecularRestraint\n");
363	<	painCave.severity = OOPSE_ERROR;
363	>	painCave.severity = OPENMD_ERROR;
364		painCave.isFatal = 1;
365		simError();
366		}
367		} else {
368		sprintf( painCave.errMsg,
369		"Can not cast GenericData to RestraintData\n");
370	<	painCave.severity = OOPSE_ERROR;
370	>	painCave.severity = OPENMD_ERROR;
371		painCave.isFatal = 1;
372		simError();
373		}
374		} else {
375		sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
376	<	painCave.severity = OOPSE_ERROR;
376	>	painCave.severity = OPENMD_ERROR;
377		painCave.isFatal = 1;
378		simError();
379		}
#	Line 338 | Line 387 | namespace oopse {
387		std::vector<Vector3d> forces;
388
389		for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
390	<	sd = (*rm)->nextIntegrableObject(ioi)) {
390	>	sd = (*rm)->nextIntegrableObject(ioi)) {
391		struc.push_back(sd->getPos());
392		}
393
#	Line 357 | Line 406 | namespace oopse {
406		unscaledPotential_ += mRest->getUnscaledPotential();
407
408		restInfo = mRest->getRestraintInfo();
409	<	restInfo_.push_back(restInfo);
409	>
410	>	// only collect data on restraints that we're going to print:
411	>	if (mRest->getPrintRestraint())
412	>	restInfo_.push_back(restInfo);
413		}
414
415		for(ro=restrainedObjs_.begin(); ro != restrainedObjs_.end(); ++ro){
#	Line 373 | Line 425 | namespace oopse {
425		if (oRest == NULL) {
426		sprintf( painCave.errMsg,
427		"Can not cast RestraintData to ObjectRestraint\n");
428	<	painCave.severity = OOPSE_ERROR;
428	>	painCave.severity = OPENMD_ERROR;
429		painCave.isFatal = 1;
430		simError();
431		}
432		} else {
433		sprintf( painCave.errMsg,
434		"Can not cast GenericData to RestraintData\n");
435	<	painCave.severity = OOPSE_ERROR;
435	>	painCave.severity = OPENMD_ERROR;
436		painCave.isFatal = 1;
437		simError();
438		}
439		} else {
440		sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
441	<	painCave.severity = OOPSE_ERROR;
441	>	painCave.severity = OPENMD_ERROR;
442		painCave.isFatal = 1;
443		simError();
444		}
#	Line 418 | Line 470 | namespace oopse {
470		unscaledPotential_ += oRest->getUnscaledPotential();
471
472		restInfo = oRest->getRestraintInfo();
473	<	restInfo_.push_back(restInfo);
473	>
474	>	// only collect data on restraints that we're going to print:
475	>	if (oRest->getPrintRestraint())
476	>	restInfo_.push_back(restInfo);
477		}
478
479		return unscaledPotential_ * scalingFactor;

Comparing trunk/src/restraints/RestraintForceManager.cpp (property svn:keywords):
Revision 1360 by cli2, Mon Sep 7 16:31:51 2009 UTC vs.
Revision 1464 by gezelter, Fri Jul 9 19:29:05 2010 UTC

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

Diff Legend

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