[ViewVC] Diff of: OpenMD/branches/development/src/brains/ForceManager.cpp

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1576 by gezelter, Wed Jun 8 16:05:07 2011 UTC vs.
Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

#	Line 36 \| Line 36
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).
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		/**
44		* @file ForceManager.cpp
45		* @author tlin
46		* @date 11/09/2004
46	–	* @time 10:39am
47		* @version 1.0
48		*/
49
#	Line 57 \| Line 57
57		#include "primitives/Torsion.hpp"
58		#include "primitives/Inversion.hpp"
59		#include "nonbonded/NonBondedInteraction.hpp"
60	+	#include "perturbations/ElectricField.hpp"
61		#include "parallel/ForceMatrixDecomposition.hpp"
62
63	+	#include <cstdio>
64	+	#include <iostream>
65	+	#include <iomanip>
66	+
67		using namespace std;
68		namespace OpenMD {
69
70		ForceManager::ForceManager(SimInfo * info) : info_(info) {
71		forceField_ = info_->getForceField();
72	<	fDecomp_ = new ForceMatrixDecomposition(info_);
72	>	interactionMan_ = new InteractionManager();
73	>	fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
74		}
75
76		/**
77		* setupCutoffs
78		*
79	<	* Sets the values of cutoffRadius, cutoffMethod, and cutoffPolicy
79	>	* Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
80	>	* and cutoffPolicy
81		*
82		* cutoffRadius : realType
83		* If the cutoffRadius was explicitly set, use that value.
#	Line 80 \| Line 87 \| namespace OpenMD {
87		* simulation for suggested cutoff values (e.g. 2.5 * sigma).
88		* Use the maximum suggested value that was found.
89		*
90	<	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
90	>	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
91	>	* or SHIFTED_POTENTIAL)
92		* If cutoffMethod was explicitly set, use that choice.
93		* If cutoffMethod was not explicitly set, use SHIFTED_FORCE
94		*
95		* cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
96		* If cutoffPolicy was explicitly set, use that choice.
97		* If cutoffPolicy was not explicitly set, use TRADITIONAL
98	+	*
99	+	* switchingRadius : realType
100	+	* If the cutoffMethod was set to SWITCHED:
101	+	* If the switchingRadius was explicitly set, use that value
102	+	* (but do a sanity check first).
103	+	* If the switchingRadius was not explicitly set: use 0.85 *
104	+	* cutoffRadius_
105	+	* If the cutoffMethod was not set to SWITCHED:
106	+	* Set switchingRadius equal to cutoffRadius for safety.
107		*/
108		void ForceManager::setupCutoffs() {
109
110		Globals* simParams_ = info_->getSimParams();
111		ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
112	+	int mdFileVersion;
113	+	rCut_ = 0.0; //Needs a value for a later max() call;
114
115	+	if (simParams_->haveMDfileVersion())
116	+	mdFileVersion = simParams_->getMDfileVersion();
117	+	else
118	+	mdFileVersion = 0;
119	+
120		if (simParams_->haveCutoffRadius()) {
121		rCut_ = simParams_->getCutoffRadius();
122		} else {
#	Line 121 \| Line 145 \| namespace OpenMD {
145		painCave.isFatal = 0;
146		painCave.severity = OPENMD_INFO;
147		simError();
148	<	}
148	>	}
149		}
150
151	+	fDecomp_->setUserCutoff(rCut_);
152	+	interactionMan_->setCutoffRadius(rCut_);
153	+
154		map<string, CutoffMethod> stringToCutoffMethod;
155		stringToCutoffMethod["HARD"] = HARD;
156		stringToCutoffMethod["SWITCHED"] = SWITCHED;
#	Line 147 \| Line 174 \| namespace OpenMD {
174		cutoffMethod_ = i->second;
175		}
176		} else {
177	<	sprintf(painCave.errMsg,
178	<	"ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
179	<	"\tOpenMD will use SHIFTED_FORCE.\n");
180	<	painCave.isFatal = 0;
181	<	painCave.severity = OPENMD_INFO;
182	<	simError();
183	<	cutoffMethod_ = SHIFTED_FORCE;
177	>	if (mdFileVersion > 1) {
178	>	sprintf(painCave.errMsg,
179	>	"ForceManager::setupCutoffs: No value was set for the cutoffMethod.\n"
180	>	"\tOpenMD will use SHIFTED_FORCE.\n");
181	>	painCave.isFatal = 0;
182	>	painCave.severity = OPENMD_INFO;
183	>	simError();
184	>	cutoffMethod_ = SHIFTED_FORCE;
185	>	} else {
186	>	// handle the case where the old file version was in play
187	>	// (there should be no cutoffMethod, so we have to deduce it
188	>	// from other data).
189	>
190	>	sprintf(painCave.errMsg,
191	>	"ForceManager::setupCutoffs : DEPRECATED FILE FORMAT!\n"
192	>	"\tOpenMD found a file which does not set a cutoffMethod.\n"
193	>	"\tOpenMD will attempt to deduce a cutoffMethod using the\n"
194	>	"\tbehavior of the older (version 1) code. To remove this\n"
195	>	"\twarning, add an explicit cutoffMethod and change the top\n"
196	>	"\tof the file so that it begins with <OpenMD version=2>\n");
197	>	painCave.isFatal = 0;
198	>	painCave.severity = OPENMD_WARNING;
199	>	simError();
200	>
201	>	// The old file version tethered the shifting behavior to the
202	>	// electrostaticSummationMethod keyword.
203	>
204	>	if (simParams_->haveElectrostaticSummationMethod()) {
205	>	string myMethod = simParams_->getElectrostaticSummationMethod();
206	>	toUpper(myMethod);
207	>
208	>	if (myMethod == "SHIFTED_POTENTIAL") {
209	>	cutoffMethod_ = SHIFTED_POTENTIAL;
210	>	} else if (myMethod == "SHIFTED_FORCE") {
211	>	cutoffMethod_ = SHIFTED_FORCE;
212	>	}
213	>
214	>	if (simParams_->haveSwitchingRadius())
215	>	rSwitch_ = simParams_->getSwitchingRadius();
216	>
217	>	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE") {
218	>	if (simParams_->haveSwitchingRadius()){
219	>	sprintf(painCave.errMsg,
220	>	"ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n"
221	>	"\tA value was set for the switchingRadius\n"
222	>	"\teven though the electrostaticSummationMethod was\n"
223	>	"\tset to %s\n", myMethod.c_str());
224	>	painCave.severity = OPENMD_WARNING;
225	>	painCave.isFatal = 1;
226	>	simError();
227	>	}
228	>	}
229	>	if (abs(rCut_ - rSwitch_) < 0.0001) {
230	>	if (cutoffMethod_ == SHIFTED_FORCE) {
231	>	sprintf(painCave.errMsg,
232	>	"ForceManager::setupCutoffs : DEPRECATED BEHAVIOR\n"
233	>	"\tcutoffRadius and switchingRadius are set to the\n"
234	>	"\tsame value. OpenMD will use shifted force\n"
235	>	"\tpotentials instead of switching functions.\n");
236	>	painCave.isFatal = 0;
237	>	painCave.severity = OPENMD_WARNING;
238	>	simError();
239	>	} else {
240	>	cutoffMethod_ = SHIFTED_POTENTIAL;
241	>	sprintf(painCave.errMsg,
242	>	"ForceManager::setupCutoffs : DEPRECATED BEHAVIOR\n"
243	>	"\tcutoffRadius and switchingRadius are set to the\n"
244	>	"\tsame value. OpenMD will use shifted potentials\n"
245	>	"\tinstead of switching functions.\n");
246	>	painCave.isFatal = 0;
247	>	painCave.severity = OPENMD_WARNING;
248	>	simError();
249	>	}
250	>	}
251	>	}
252	>	}
253		}
254
255		map<string, CutoffPolicy> stringToCutoffPolicy;
#	Line 161 \| Line 257 \| namespace OpenMD {
257		stringToCutoffPolicy["MAX"] = MAX;
258		stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;
259
260	<	std::string cutPolicy;
260	>	string cutPolicy;
261		if (forceFieldOptions_.haveCutoffPolicy()){
262		cutPolicy = forceFieldOptions_.getCutoffPolicy();
263		}else if (simParams_->haveCutoffPolicy()) {
#	Line 194 \| Line 290 \| namespace OpenMD {
290		simError();
291		cutoffPolicy_ = TRADITIONAL;
292		}
197	–	}
293
294	<	/**
295	<	* setupSwitching
296	<	*
297	<	* Sets the values of switchingRadius and
298	<	* If the switchingRadius was explicitly set, use that value (but check it)
299	<	* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
300	<	*/
301	<	void ForceManager::setupSwitching() {
302	<	Globals* simParams_ = info_->getSimParams();
303	<
304	<	if (simParams_->haveSwitchingRadius()) {
305	<	rSwitch_ = simParams_->getSwitchingRadius();
306	<	if (rSwitch_ > rCut_) {
294	>	fDecomp_->setCutoffPolicy(cutoffPolicy_);
295	>
296	>	// create the switching function object:
297	>
298	>	switcher_ = new SwitchingFunction();
299	>
300	>	if (cutoffMethod_ == SWITCHED) {
301	>	if (simParams_->haveSwitchingRadius()) {
302	>	rSwitch_ = simParams_->getSwitchingRadius();
303	>	if (rSwitch_ > rCut_) {
304	>	sprintf(painCave.errMsg,
305	>	"ForceManager::setupCutoffs: switchingRadius (%f) is larger "
306	>	"than the cutoffRadius(%f)\n", rSwitch_, rCut_);
307	>	painCave.isFatal = 1;
308	>	painCave.severity = OPENMD_ERROR;
309	>	simError();
310	>	}
311	>	} else {
312	>	rSwitch_ = 0.85 * rCut_;
313		sprintf(painCave.errMsg,
314	<	"ForceManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
315	<	rSwitch_, rCut_);
316	<	painCave.isFatal = 1;
317	<	painCave.severity = OPENMD_ERROR;
314	>	"ForceManager::setupCutoffs: No value was set for the switchingRadius.\n"
315	>	"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
316	>	"\tswitchingRadius = %f. for this simulation\n", rSwitch_);
317	>	painCave.isFatal = 0;
318	>	painCave.severity = OPENMD_WARNING;
319		simError();
320		}
321	<	} else {
322	<	rSwitch_ = 0.85 * rCut_;
323	<	sprintf(painCave.errMsg,
324	<	"ForceManager::setupSwitching: No value was set for the switchingRadius.\n"
325	<	"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
326	<	"\tswitchingRadius = %f. for this simulation\n", rSwitch_);
327	<	painCave.isFatal = 0;
328	<	painCave.severity = OPENMD_WARNING;
329	<	simError();
330	<	}
321	>	} else {
322	>	if (mdFileVersion > 1) {
323	>	// throw an error if we define a switching radius and don't need one.
324	>	// older file versions should not do this.
325	>	if (simParams_->haveSwitchingRadius()) {
326	>	map<string, CutoffMethod>::const_iterator it;
327	>	string theMeth;
328	>	for (it = stringToCutoffMethod.begin();
329	>	it != stringToCutoffMethod.end(); ++it) {
330	>	if (it->second == cutoffMethod_) {
331	>	theMeth = it->first;
332	>	break;
333	>	}
334	>	}
335	>	sprintf(painCave.errMsg,
336	>	"ForceManager::setupCutoffs: the cutoffMethod (%s)\n"
337	>	"\tis not set to SWITCHED, so switchingRadius value\n"
338	>	"\twill be ignored for this simulation\n", theMeth.c_str());
339	>	painCave.isFatal = 0;
340	>	painCave.severity = OPENMD_WARNING;
341	>	simError();
342	>	}
343	>	}
344	>	rSwitch_ = rCut_;
345	>	}
346
347	+	// Default to cubic switching function.
348	+	sft_ = cubic;
349		if (simParams_->haveSwitchingFunctionType()) {
350		string funcType = simParams_->getSwitchingFunctionType();
351		toUpper(funcType);
#	Line 251 \| Line 370 \| namespace OpenMD {
370		switcher_->setSwitchType(sft_);
371		switcher_->setSwitch(rSwitch_, rCut_);
372		}
373	+
374	+
375	+
376
377		void ForceManager::initialize() {
378
379		if (!info_->isTopologyDone()) {
380	+
381		info_->update();
382		interactionMan_->setSimInfo(info_);
383		interactionMan_->initialize();
#	Line 262 \| Line 385 \| namespace OpenMD {
385		// We want to delay the cutoffs until after the interaction
386		// manager has set up the atom-atom interactions so that we can
387		// query them for suggested cutoff values
265	–
388		setupCutoffs();
267	–	setupSwitching();
389
390		info_->prepareTopology();
391	+
392	+	doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
393	+	doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
394	+	if (doHeatFlux_) doParticlePot_ = true;
395	+
396	+	doElectricField_ = info_->getSimParams()->getOutputElectricField();
397	+
398		}
399
400		ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
401
402	<	// Force fields can set options on how to scale van der Waals and electrostatic
403	<	// interactions for atoms connected via bonds, bends and torsions
404	<	// in this case the topological distance between atoms is:
402	>	// Force fields can set options on how to scale van der Waals and
403	>	// electrostatic interactions for atoms connected via bonds, bends
404	>	// and torsions in this case the topological distance between
405	>	// atoms is:
406		// 0 = topologically unconnected
407		// 1 = bonded together
408		// 2 = connected via a bend
#	Line 295 \| Line 424 \| namespace OpenMD {
424		electrostaticScale_[2] = fopts.getelectrostatic13scale();
425		electrostaticScale_[3] = fopts.getelectrostatic14scale();
426
427	+	if (info_->getSimParams()->haveElectricField()) {
428	+	ElectricField* eField = new ElectricField(info_);
429	+	perturbations_.push_back(eField);
430	+	}
431	+
432		fDecomp_->distributeInitialData();
433
434		initialized_ = true;
#	Line 321 \| Line 455 \| namespace OpenMD {
455		Molecule::CutoffGroupIterator ci;
456		CutoffGroup* cg;
457
458	<	// forces are zeroed here, before any are accumulated.
458	>	// forces and potentials are zeroed here, before any are
459	>	// accumulated.
460
461	+	Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
462	+
463	+	snap->setBondPotential(0.0);
464	+	snap->setBendPotential(0.0);
465	+	snap->setTorsionPotential(0.0);
466	+	snap->setInversionPotential(0.0);
467	+
468	+	potVec zeroPot(0.0);
469	+	snap->setLongRangePotential(zeroPot);
470	+	snap->setExcludedPotentials(zeroPot);
471	+
472	+	snap->setRestraintPotential(0.0);
473	+	snap->setRawPotential(0.0);
474	+
475		for (mol = info_->beginMolecule(mi); mol != NULL;
476		mol = info_->nextMolecule(mi)) {
477	<	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
477	>	for(atom = mol->beginAtom(ai); atom != NULL;
478	>	atom = mol->nextAtom(ai)) {
479		atom->zeroForcesAndTorques();
480		}
481	<
481	>
482		//change the positions of atoms which belong to the rigidbodies
483		for (rb = mol->beginRigidBody(rbIter); rb != NULL;
484		rb = mol->nextRigidBody(rbIter)) {
485		rb->zeroForcesAndTorques();
486		}
487	<
487	>
488		if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
489		for(cg = mol->beginCutoffGroup(ci); cg != NULL;
490		cg = mol->nextCutoffGroup(ci)) {
#	Line 343 \| Line 493 \| namespace OpenMD {
493		}
494		}
495		}
346	–
347	–	// Zero out the stress tensor
348	–	tau *= 0.0;
496
497	+	// Zero out the stress tensor
498	+	stressTensor *= 0.0;
499	+	// Zero out the heatFlux
500	+	fDecomp_->setHeatFlux( Vector3d(0.0) );
501		}
502
503		void ForceManager::shortRangeInteractions() {
#	Line 379 \| Line 530 \| namespace OpenMD {
530
531		for (bond = mol->beginBond(bondIter); bond != NULL;
532		bond = mol->nextBond(bondIter)) {
533	<	bond->calcForce();
533	>	bond->calcForce(doParticlePot_);
534		bondPotential += bond->getPotential();
535		}
536
#	Line 387 \| Line 538 \| namespace OpenMD {
538		bend = mol->nextBend(bendIter)) {
539
540		RealType angle;
541	<	bend->calcForce(angle);
541	>	bend->calcForce(angle, doParticlePot_);
542		RealType currBendPot = bend->getPotential();
543
544		bendPotential += bend->getPotential();
#	Line 397 \| Line 548 \| namespace OpenMD {
548		dataSet.prev.angle = dataSet.curr.angle = angle;
549		dataSet.prev.potential = dataSet.curr.potential = currBendPot;
550		dataSet.deltaV = 0.0;
551	<	bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend, dataSet));
551	>	bendDataSets.insert(map<Bend*, BendDataSet>::value_type(bend,
552	>	dataSet));
553		}else {
554		i->second.prev.angle = i->second.curr.angle;
555		i->second.prev.potential = i->second.curr.potential;
#	Line 411 \| Line 563 \| namespace OpenMD {
563		for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
564		torsion = mol->nextTorsion(torsionIter)) {
565		RealType angle;
566	<	torsion->calcForce(angle);
566	>	torsion->calcForce(angle, doParticlePot_);
567		RealType currTorsionPot = torsion->getPotential();
568		torsionPotential += torsion->getPotential();
569		map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
#	Line 435 \| Line 587 \| namespace OpenMD {
587		inversion != NULL;
588		inversion = mol->nextInversion(inversionIter)) {
589		RealType angle;
590	<	inversion->calcForce(angle);
590	>	inversion->calcForce(angle, doParticlePot_);
591		RealType currInversionPot = inversion->getPotential();
592		inversionPotential += inversion->getPotential();
593		map<Inversion*, InversionDataSet>::iterator i = inversionDataSets.find(inversion);
#	Line 455 \| Line 607 \| namespace OpenMD {
607		}
608		}
609		}
610	<
611	<	RealType shortRangePotential = bondPotential + bendPotential +
612	<	torsionPotential + inversionPotential;
610	>
611	>	#ifdef IS_MPI
612	>	// Collect from all nodes. This should eventually be moved into a
613	>	// SystemDecomposition, but this is a better place than in
614	>	// Thermo to do the collection.
615	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
616	>	MPI::SUM);
617	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
618	>	MPI::SUM);
619	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
620	>	MPI::REALTYPE, MPI::SUM);
621	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
622	>	MPI::REALTYPE, MPI::SUM);
623	>	#endif
624	>
625		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
626	<	curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
627	<	curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
628	<	curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
629	<	curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
630	<	curSnapshot->statData[Stats::INVERSION_POTENTIAL] = inversionPotential;
626	>
627	>	curSnapshot->setBondPotential(bondPotential);
628	>	curSnapshot->setBendPotential(bendPotential);
629	>	curSnapshot->setTorsionPotential(torsionPotential);
630	>	curSnapshot->setInversionPotential(inversionPotential);
631	>
632	>	// RealType shortRangePotential = bondPotential + bendPotential +
633	>	// torsionPotential + inversionPotential;
634	>
635	>	// curSnapshot->setShortRangePotential(shortRangePotential);
636		}
637
638		void ForceManager::longRangeInteractions() {
639
640	<	// some of this initial stuff will go away:
640	>
641		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
642		DataStorage* config = &(curSnapshot->atomData);
643		DataStorage* cgConfig = &(curSnapshot->cgData);
475	–	RealType* frc = config->getArrayPointer(DataStorage::dslForce);
476	–	RealType* pos = config->getArrayPointer(DataStorage::dslPosition);
477	–	RealType* trq = config->getArrayPointer(DataStorage::dslTorque);
478	–	RealType* A = config->getArrayPointer(DataStorage::dslAmat);
479	–	RealType* electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame);
480	–	RealType* particlePot = config->getArrayPointer(DataStorage::dslParticlePot);
481	–	RealType* rc;
644
645	<	if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){
646	<	rc = cgConfig->getArrayPointer(DataStorage::dslPosition);
645	>	//calculate the center of mass of cutoff group
646	>
647	>	SimInfo::MoleculeIterator mi;
648	>	Molecule* mol;
649	>	Molecule::CutoffGroupIterator ci;
650	>	CutoffGroup* cg;
651	>
652	>	if(info_->getNCutoffGroups() > 0){
653	>	for (mol = info_->beginMolecule(mi); mol != NULL;
654	>	mol = info_->nextMolecule(mi)) {
655	>	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
656	>	cg = mol->nextCutoffGroup(ci)) {
657	>	cg->updateCOM();
658	>	}
659	>	}
660		} else {
661		// center of mass of the group is the same as position of the atom
662		// if cutoff group does not exist
663	<	rc = pos;
663	>	cgConfig->position = config->position;
664	>	cgConfig->velocity = config->velocity;
665		}
666	<
491	<	// new stuff starts here:
666	>
667		fDecomp_->zeroWorkArrays();
668		fDecomp_->distributeData();
669	<
670	<	int cg1, cg2, atom1, atom2;
671	<	Vector3d d_grp, dag;
672	<	RealType rgrpsq, rgrp;
669	>
670	>	int cg1, cg2, atom1, atom2, topoDist;
671	>	Vector3d d_grp, dag, d, gvel2, vel2;
672	>	RealType rgrpsq, rgrp, r2, r;
673	>	RealType electroMult, vdwMult;
674		RealType vij;
675	<	Vector3d fij, fg;
675	>	Vector3d fij, fg, f1;
676		tuple3<RealType, RealType, RealType> cuts;
677		RealType rCutSq;
678		bool in_switching_region;
#	Line 505 \| Line 681 \| namespace OpenMD {
681		InteractionData idat;
682		SelfData sdat;
683		RealType mf;
684	<	potVec pot(0.0);
684	>	RealType vpair;
685	>	RealType dVdFQ1(0.0);
686	>	RealType dVdFQ2(0.0);
687		potVec longRangePotential(0.0);
688	<	RealType lrPot;
688	>	potVec workPot(0.0);
689	>	potVec exPot(0.0);
690	>	Vector3d eField1(0.0);
691	>	Vector3d eField2(0.0);
692	>	vector<int>::iterator ia, jb;
693
694		int loopStart, loopEnd;
695
696	+	idat.vdwMult = &vdwMult;
697	+	idat.electroMult = &electroMult;
698	+	idat.pot = &workPot;
699	+	idat.excludedPot = &exPot;
700	+	sdat.pot = fDecomp_->getEmbeddingPotential();
701	+	sdat.excludedPot = fDecomp_->getExcludedSelfPotential();
702	+	idat.vpair = &vpair;
703	+	idat.dVdFQ1 = &dVdFQ1;
704	+	idat.dVdFQ2 = &dVdFQ2;
705	+	idat.eField1 = &eField1;
706	+	idat.eField2 = &eField2;
707	+	idat.f1 = &f1;
708	+	idat.sw = &sw;
709	+	idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
710	+	idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
711	+	idat.doParticlePot = doParticlePot_;
712	+	idat.doElectricField = doElectricField_;
713	+	sdat.doParticlePot = doParticlePot_;
714	+
715		loopEnd = PAIR_LOOP;
716		if (info_->requiresPrepair() ) {
717		loopStart = PREPAIR_LOOP;
718		} else {
719		loopStart = PAIR_LOOP;
720		}
721	<
722	<	for (int iLoop = loopStart; iLoop < loopEnd; iLoop++) {
522	<
721	>	for (int iLoop = loopStart; iLoop <= loopEnd; iLoop++) {
722	>
723		if (iLoop == loopStart) {
724		bool update_nlist = fDecomp_->checkNeighborList();
725		if (update_nlist)
726		neighborList = fDecomp_->buildNeighborList();
727	<	}
727	>	}
728
729		for (vector<pair<int, int> >::iterator it = neighborList.begin();
730		it != neighborList.end(); ++it) {
731	<
731	>
732		cg1 = (*it).first;
733		cg2 = (*it).second;
734
735		cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
736
737		d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
738	+
739		curSnapshot->wrapVector(d_grp);
740		rgrpsq = d_grp.lengthSquare();
540	–
741		rCutSq = cuts.second;
742
743		if (rgrpsq < rCutSq) {
744	<	*(idat.rcut) = cuts.first;
744	>	idat.rcut = &cuts.first;
745		if (iLoop == PAIR_LOOP) {
746	<	vij *= 0.0;
746	>	vij = 0.0;
747		fij = V3Zero;
748		}
749
750	<	in_switching_region = switcher_->getSwitch(rgrpsq, *(idat.sw), dswdr,
750	>	in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
751		rgrp);
752	<
752	>
753		atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
754		atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
755
756	<	for (vector<int>::iterator ia = atomListRow.begin();
756	>	if (doHeatFlux_)
757	>	gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
758	>
759	>	for (ia = atomListRow.begin();
760		ia != atomListRow.end(); ++ia) {
761		atom1 = (*ia);
762	<
763	<	for (vector<int>::iterator jb = atomListColumn.begin();
762	>
763	>	for (jb = atomListColumn.begin();
764		jb != atomListColumn.end(); ++jb) {
765		atom2 = (*jb);
563	–
564	–	if (!fDecomp_->skipAtomPair(atom1, atom2)) {
565	–
566	–	pot *= 0.0;
766
767	<	idat = fDecomp_->fillInteractionData(atom1, atom2);
569	<	*(idat.pot) = pot;
767	>	if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
768
769	+	vpair = 0.0;
770	+	workPot = 0.0;
771	+	exPot = 0.0;
772	+	f1 = V3Zero;
773	+	dVdFQ1 = 0.0;
774	+	dVdFQ2 = 0.0;
775	+
776	+	fDecomp_->fillInteractionData(idat, atom1, atom2);
777	+
778	+	topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
779	+	vdwMult = vdwScale_[topoDist];
780	+	electroMult = electrostaticScale_[topoDist];
781	+
782		if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
783	<	*(idat.d) = d_grp;
784	<	*(idat.r2) = rgrpsq;
783	>	idat.d = &d_grp;
784	>	idat.r2 = &rgrpsq;
785	>	if (doHeatFlux_)
786	>	vel2 = gvel2;
787		} else {
788	<	*(idat.d) = fDecomp_->getInteratomicVector(atom1, atom2);
789	<	curSnapshot->wrapVector( *(idat.d) );
790	<	*(idat.r2) = idat.d->lengthSquare();
788	>	d = fDecomp_->getInteratomicVector(atom1, atom2);
789	>	curSnapshot->wrapVector( d );
790	>	r2 = d.lengthSquare();
791	>	idat.d = &d;
792	>	idat.r2 = &r2;
793	>	if (doHeatFlux_)
794	>	vel2 = fDecomp_->getAtomVelocityColumn(atom2);
795		}
579	–
580	–	(idat.rij) = sqrt( (idat.r2) );
796
797	+	r = sqrt( *(idat.r2) );
798	+	idat.rij = &r;
799	+
800		if (iLoop == PREPAIR_LOOP) {
801		interactionMan_->doPrePair(idat);
802		} else {
803		interactionMan_->doPair(idat);
804		fDecomp_->unpackInteractionData(idat, atom1, atom2);
805	<	vij += *(idat.vpair);
806	<	fij += *(idat.f1);
807	<	tau -= outProduct( (idat.d), (idat.f1));
805	>	vij += vpair;
806	>	fij += f1;
807	>	stressTensor -= outProduct( *(idat.d), f1);
808	>	if (doHeatFlux_)
809	>	fDecomp_->addToHeatFlux((idat.d) dot(f1, vel2));
810		}
811		}
812		}
#	Line 596 \| Line 816 \| namespace OpenMD {
816		if (in_switching_region) {
817		swderiv = vij * dswdr / rgrp;
818		fg = swderiv * d_grp;
599	–
819		fij += fg;
820
821		if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
822	<	tau -= outProduct( *(idat.d), fg);
822	>	stressTensor -= outProduct( *(idat.d), fg);
823	>	if (doHeatFlux_)
824	>	fDecomp_->addToHeatFlux((idat.d) dot(fg, vel2));
825	>
826		}
827
828	<	for (vector<int>::iterator ia = atomListRow.begin();
828	>	for (ia = atomListRow.begin();
829		ia != atomListRow.end(); ++ia) {
830		atom1 = (*ia);
831		mf = fDecomp_->getMassFactorRow(atom1);
#	Line 611 \| Line 833 \| namespace OpenMD {
833		// presence in switching region
834		fg = swderiv * d_grp * mf;
835		fDecomp_->addForceToAtomRow(atom1, fg);
614	–
836		if (atomListRow.size() > 1) {
837		if (info_->usesAtomicVirial()) {
838		// find the distance between the atom
839		// and the center of the cutoff group:
840		dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
841	<	tau -= outProduct(dag, fg);
841	>	stressTensor -= outProduct(dag, fg);
842	>	if (doHeatFlux_)
843	>	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
844		}
845		}
846		}
847	<	for (vector<int>::iterator jb = atomListColumn.begin();
847	>	for (jb = atomListColumn.begin();
848		jb != atomListColumn.end(); ++jb) {
849		atom2 = (*jb);
850		mf = fDecomp_->getMassFactorColumn(atom2);
#	Line 635 \| Line 858 \| namespace OpenMD {
858		// find the distance between the atom
859		// and the center of the cutoff group:
860		dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
861	<	tau -= outProduct(dag, fg);
861	>	stressTensor -= outProduct(dag, fg);
862	>	if (doHeatFlux_)
863	>	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
864		}
865		}
866		}
867		}
868	<	//if (!SIM_uses_AtomicVirial) {
869	<	// tau -= outProduct(d_grp, fij);
868	>	//if (!info_->usesAtomicVirial()) {
869	>	// stressTensor -= outProduct(d_grp, fij);
870	>	// if (doHeatFlux_)
871	>	// fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
872		//}
873		}
874		}
875		}
876
877		if (iLoop == PREPAIR_LOOP) {
878	<	if (info_->requiresPrepair()) {
878	>	if (info_->requiresPrepair()) {
879	>
880		fDecomp_->collectIntermediateData();
881
882	<	for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
883	<	sdat = fDecomp_->fillSelfData(atom1);
882	>	for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
883	>	fDecomp_->fillSelfData(sdat, atom1);
884		interactionMan_->doPreForce(sdat);
885		}
886
887	<	fDecomp_->distributeIntermediateData();
887	>	fDecomp_->distributeIntermediateData();
888	>
889		}
890		}
662	–
891		}
892
893	+	// collects pairwise information
894		fDecomp_->collectData();
666	–
667	–	if ( info_->requiresSkipCorrection() ) {
668	–
669	–	for (int atom1 = 0; atom1 < fDecomp_->getNAtomsInRow(); atom1++) {
670	–
671	–	vector<int> skipList = fDecomp_->getSkipsForRowAtom( atom1 );
895
673	–	for (vector<int>::iterator jb = skipList.begin();
674	–	jb != skipList.end(); ++jb) {
675	–
676	–	atom2 = (*jb);
677	–	idat = fDecomp_->fillSkipData(atom1, atom2);
678	–	interactionMan_->doSkipCorrection(idat);
679	–
680	–	}
681	–	}
682	–	}
683	–
896		if (info_->requiresSelfCorrection()) {
897	<
898	<	for (int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
687	<	sdat = fDecomp_->fillSelfData(atom1);
897	>	for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
898	>	fDecomp_->fillSelfData(sdat, atom1);
899		interactionMan_->doSelfCorrection(sdat);
900		}
690	–
901		}
902
903	<	longRangePotential = fDecomp_->getLongRangePotential();
904	<	lrPot = longRangePotential.sum();
903	>	// collects single-atom information
904	>	fDecomp_->collectSelfData();
905
906	<	//store the tau and long range potential
907	<	curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
908	<	curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VANDERWAALS_FAMILY];
909	<	curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_FAMILY];
906	>	longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
907	>	*(fDecomp_->getPairwisePotential());
908	>
909	>	curSnapshot->setLongRangePotential(longRangePotential);
910	>
911	>	// collects single-atom information
912	>	fDecomp_->collectSelfData();
913	>
914	>	longRangePotential = *(fDecomp_->getEmbeddingPotential()) +
915	>	*(fDecomp_->getPairwisePotential());
916	>
917	>	curSnapshot->setLongRangePotential(longRangePotential);
918	>
919	>	curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
920	>	*(fDecomp_->getExcludedPotential()));
921	>
922		}
923
924
925		void ForceManager::postCalculation() {
926	+
927	+	vector<Perturbation*>::iterator pi;
928	+	for (pi = perturbations_.begin(); pi != perturbations_.end(); ++pi) {
929	+	(*pi)->applyPerturbation();
930	+	}
931	+
932		SimInfo::MoleculeIterator mi;
933		Molecule* mol;
934		Molecule::RigidBodyIterator rbIter;
935		RigidBody* rb;
936		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
937	<
937	>
938		// collect the atomic forces onto rigid bodies
939
940		for (mol = info_->beginMolecule(mi); mol != NULL;
#	Line 714 \| Line 942 \| namespace OpenMD {
942		for (rb = mol->beginRigidBody(rbIter); rb != NULL;
943		rb = mol->nextRigidBody(rbIter)) {
944		Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial();
945	<	tau += rbTau;
945	>	stressTensor += rbTau;
946		}
947		}
948
949		#ifdef IS_MPI
950	<	Mat3x3d tmpTau(tau);
951	<	MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(),
724	<	9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
950	>	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
951	>	MPI::REALTYPE, MPI::SUM);
952		#endif
953	<	curSnapshot->statData.setTau(tau);
953	>	curSnapshot->setStressTensor(stressTensor);
954	>
955		}
728	–
956		} //end namespace OpenMD

Diff Legend

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

Comparing branches/development/src/brains/ForceManager.cpp (file contents): Revision 1576 by gezelter, Wed Jun 8 16:05:07 2011 UTC vs. Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC

Diff Legend

Comparing branches/development/src/brains/ForceManager.cpp (file contents):
Revision 1576 by gezelter, Wed Jun 8 16:05:07 2011 UTC vs.
Revision 1808 by gezelter, Mon Oct 22 20:42:10 2012 UTC