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 1364 by cli2, Wed Oct 7 20:49:50 2009 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 "config.h"
44		#include <cmath>
45	+
46		#include "restraints/RestraintForceManager.hpp"
47		#include "restraints/MolecularRestraint.hpp"
48		#include "restraints/ObjectRestraint.hpp"
49		#include "io/RestReader.hpp"
50		#include "utils/simError.h"
51	<	#include "utils/OOPSEConstant.hpp"
51	>	#include "utils/PhysicalConstants.hpp"
52		#include "utils/StringUtils.hpp"
53		#include "selection/SelectionEvaluator.hpp"
54		#include "selection/SelectionManager.hpp"
#	Line 54 | Line 57
57		#endif
58
59
60	<	namespace oopse {
60	>	namespace OpenMD {
61
62		RestraintForceManager::RestraintForceManager(SimInfo* info): ForceManager(info) {
63
#	Line 89 | Line 92 | namespace oopse {
92		int nRestraintStamps = simParam->getNRestraintStamps();
93		std::vector<RestraintStamp*> stamp = simParam->getRestraintStamps();
94
95	+	std::vector<int> stuntDoubleIndex;
96	+
97		for (int i = 0; i < nRestraintStamps; i++){
98
99		std::string myType = toUpperCopy(stamp[i]->getType());
#	Line 96 | Line 101 | namespace oopse {
101		if (myType.compare("MOLECULAR")==0){
102
103		int molIndex;
99	–	std::vector<Vector3d> ref;
104		Vector3d refCom;
105
106		if (!stamp[i]->haveMolIndex()) {
#	Line 109 | Line 113 | namespace oopse {
113		molIndex = stamp[i]->getMolIndex();
114		}
115
116	<	Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
113	<
114	<	if (mol == NULL) {
116	>	if (molIndex < 0) {
117		sprintf(painCave.errMsg,
118	<	"Restraint Error: A molecular restraint was specified, but\n"
119	<	"\tno molecule was found with global index %d.\n",
118	<	molIndex);
118	>	"Restraint Error: A molecular restraint was specified\n"
119	>	"\twith a molIndex that was less than 0\n");
120		painCave.isFatal = 1;
121		simError();
122		}
123	+	if (molIndex >= info_->getNGlobalMolecules()) {
124	+	sprintf(painCave.errMsg,
125	+	"Restraint Error: A molecular restraint was specified with\n"
126	+	"\ta molIndex that was greater than the total number of molecules\n");
127	+	painCave.isFatal = 1;
128	+	simError();
129	+	}
130	+
131	+	Molecule* mol = info_->getMoleculeByGlobalIndex(molIndex);
132
133	+	if (mol == NULL) {
134	+	#ifdef IS_MPI
135	+	// getMoleculeByGlobalIndex returns a NULL in parallel if
136	+	// this proc doesn't have the molecule.  Do a quick check to
137	+	// make sure another processor is supposed to have it.
138	+
139	+	int myrank = MPI::COMM_WORLD.Get_rank();
140	+	if (info_->getMolToProc(molIndex) == myrank) {
141	+
142	+	// If we were supposed to have it but got a null, then freak out.
143	+	#endif
144	+
145	+	sprintf(painCave.errMsg,
146	+	"Restraint Error: A molecular restraint was specified, but\n"
147	+	"\tno molecule was found with global index %d.\n",
148	+	molIndex);
149	+	painCave.isFatal = 1;
150	+	simError();
151	+
152	+	#ifdef IS_MPI
153	+	}
154	+	#endif
155	+	}
156	+
157	+
158	+	#ifdef IS_MPI
159	+	// only handle this molecular restraint if this processor owns the
160	+	// molecule
161	+	int myrank = MPI::COMM_WORLD.Get_rank();
162	+	if (info_->getMolToProc(molIndex) == myrank) {
163	+
164	+	#endif
165	+
166		MolecularRestraint* rest = new MolecularRestraint();
167
168		std::string restPre("mol_");
#	Line 155 | Line 198 | namespace oopse {
198		restraints_.push_back(rest);
199		mol->addProperty(new RestraintData("Restraint", rest));
200		restrainedMols_.push_back(mol);
201	<
201	>	#ifdef IS_MPI
202	>	}
203	>	#endif
204		} else if (myType.compare("OBJECT") == 0) {
205
206		std::string objectSelection;
#	Line 191 | Line 236 | namespace oopse {
236
237		for (sd = seleMan.beginSelected(selei); sd != NULL;
238		sd = seleMan.nextSelected(selei)) {
239	<
239	>	stuntDoubleIndex.push_back(sd->getGlobalIntegrableObjectIndex());
240	>
241		ObjectRestraint* rest = new ObjectRestraint();
242
243		if (stamp[i]->haveDisplacementSpringConstant()) {
#	Line 231 | Line 277 | namespace oopse {
277		// are times when it won't use restraints at all, so only open the
278		// restraint file if we are actually using restraints:
279
280	<	if (simParam->getUseRestraints()) {
280	>	if (simParam->getUseRestraints()) {
281		std::string refFile = simParam->getRestraint_file();
282	<	RestReader* rr = new RestReader(info, refFile);
237	<
282	>	RestReader* rr = new RestReader(info, refFile, stuntDoubleIndex);
283		rr->readReferenceStructure();
284		}
285
286		restOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".rest";
287		restOut = new RestWriter(info_, restOutput_.c_str(), restraints_);
243	–
288		if(!restOut){
289		sprintf(painCave.errMsg, "Restraint error: Failed to create RestWriter\n");
290		painCave.isFatal = 1;
#	Line 263 | Line 307 | namespace oopse {
307		currRestTime_ = currSnapshot_->getTime();
308		}
309
310	<	void RestraintForceManager::calcForces(bool needPotential, bool needStress){
310	>	void RestraintForceManager::calcForces(){
311
312	<	ForceManager::calcForces(needPotential, needStress);
312	>	ForceManager::calcForces();
313		RealType restPot_local, restPot;
314
315		restPot_local = doRestraints(1.0);
#	Line 277 | Line 321 | namespace oopse {
321		restPot = restPot_local;
322		#endif
323		currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
324	<	currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += restPot;
325	<	currSnapshot_->statData[Stats::VHARM] = restPot;
324	>	RealType pot = currSnapshot_->getLongRangePotential();
325	>	pot += restPot;
326	>	currSnapshot_->setLongRangePotential(pot);
327	>	currSnapshot_->setRestraintPotential(restPot);
328
329		//write out forces and current positions of restrained molecules
330		if (currSnapshot_->getTime() >= currRestTime_){
#	Line 293 | Line 339 | namespace oopse {
339		Molecule::IntegrableObjectIterator ioi;
340		MolecularRestraint* mRest;
341		StuntDouble* sd;
296	–	RealType pTot;
342
343		std::vector<StuntDouble*>::const_iterator ro;
344		ObjectRestraint* oRest;
#	Line 318 | Line 363 | namespace oopse {
363		if (mRest == NULL) {
364		sprintf( painCave.errMsg,
365		"Can not cast RestraintData to MolecularRestraint\n");
366	<	painCave.severity = OOPSE_ERROR;
366	>	painCave.severity = OPENMD_ERROR;
367		painCave.isFatal = 1;
368		simError();
369		}
370		} else {
371		sprintf( painCave.errMsg,
372		"Can not cast GenericData to RestraintData\n");
373	<	painCave.severity = OOPSE_ERROR;
373	>	painCave.severity = OPENMD_ERROR;
374		painCave.isFatal = 1;
375		simError();
376		}
377		} else {
378		sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
379	<	painCave.severity = OOPSE_ERROR;
379	>	painCave.severity = OPENMD_ERROR;
380		painCave.isFatal = 1;
381		simError();
382		}
#	Line 345 | Line 390 | namespace oopse {
390		std::vector<Vector3d> forces;
391
392		for(sd = (*rm)->beginIntegrableObject(ioi); sd != NULL;
393	<	sd = (*rm)->nextIntegrableObject(ioi)) {
393	>	sd = (*rm)->nextIntegrableObject(ioi)) {
394		struc.push_back(sd->getPos());
395		}
396
#	Line 383 | Line 428 | namespace oopse {
428		if (oRest == NULL) {
429		sprintf( painCave.errMsg,
430		"Can not cast RestraintData to ObjectRestraint\n");
431	<	painCave.severity = OOPSE_ERROR;
431	>	painCave.severity = OPENMD_ERROR;
432		painCave.isFatal = 1;
433		simError();
434		}
435		} else {
436		sprintf( painCave.errMsg,
437		"Can not cast GenericData to RestraintData\n");
438	<	painCave.severity = OOPSE_ERROR;
438	>	painCave.severity = OPENMD_ERROR;
439		painCave.isFatal = 1;
440		simError();
441		}
442		} else {
443		sprintf( painCave.errMsg, "Can not find Restraint for RestrainedObject\n");
444	<	painCave.severity = OOPSE_ERROR;
444	>	painCave.severity = OPENMD_ERROR;
445		painCave.isFatal = 1;
446		simError();
447		}

Comparing trunk/src/restraints/RestraintForceManager.cpp (property svn:keywords):
Revision 1364 by cli2, Wed Oct 7 20:49:50 2009 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