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

Comparing trunk/src/brains/ForceManager.cpp (file contents):
Revision 1788 by gezelter, Wed Aug 29 20:17:07 2012 UTC vs.
Revision 2057 by gezelter, Tue Mar 3 15:22:26 2015 UTC

#	Line 35 \| Line 35
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).
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		*/
#	Line 44 \| Line 44
44		* @file ForceManager.cpp
45		* @author tlin
46		* @date 11/09/2004
47	–	* @time 10:39am
47		* @version 1.0
48		*/
49
#	Line 58 \| Line 57
57		#include "primitives/Torsion.hpp"
58		#include "primitives/Inversion.hpp"
59		#include "nonbonded/NonBondedInteraction.hpp"
60	<	#include "perturbations/ElectricField.hpp"
60	>	#include "perturbations/UniformField.hpp"
61	>	#include "perturbations/UniformGradient.hpp"
62		#include "parallel/ForceMatrixDecomposition.hpp"
63
64		#include <cstdio>
#	Line 68 \| Line 68 \| namespace OpenMD {
68		using namespace std;
69		namespace OpenMD {
70
71	<	ForceManager::ForceManager(SimInfo * info) : info_(info) {
71	>	ForceManager::ForceManager(SimInfo * info) : info_(info), switcher_(NULL),
72	>	initialized_(false) {
73		forceField_ = info_->getForceField();
74		interactionMan_ = new InteractionManager();
75		fDecomp_ = new ForceMatrixDecomposition(info_, interactionMan_);
76	+	thermo = new Thermo(info_);
77		}
78
79	+	ForceManager::~ForceManager() {
80	+	perturbations_.clear();
81	+
82	+	delete switcher_;
83	+	delete interactionMan_;
84	+	delete fDecomp_;
85	+	delete thermo;
86	+	}
87	+
88		/**
89		* setupCutoffs
90		*
91	<	* Sets the values of cutoffRadius, switchingRadius, cutoffMethod,
81	<	* and cutoffPolicy
91	>	* Sets the values of cutoffRadius, switchingRadius, and cutoffMethod
92		*
93		* cutoffRadius : realType
94		* If the cutoffRadius was explicitly set, use that value.
#	Line 88 \| Line 98 \| namespace OpenMD {
98		* simulation for suggested cutoff values (e.g. 2.5 * sigma).
99		* Use the maximum suggested value that was found.
100		*
101	<	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE,
102	<	* or SHIFTED_POTENTIAL)
101	>	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, TAYLOR_SHIFTED,
102	>	* SHIFTED_POTENTIAL, or EWALD_FULL)
103		* If cutoffMethod was explicitly set, use that choice.
104		* If cutoffMethod was not explicitly set, use SHIFTED_FORCE
105		*
96	–	* cutoffPolicy : (one of MIX, MAX, TRADITIONAL)
97	–	* If cutoffPolicy was explicitly set, use that choice.
98	–	* If cutoffPolicy was not explicitly set, use TRADITIONAL
99	–	*
106		* switchingRadius : realType
107		* If the cutoffMethod was set to SWITCHED:
108		* If the switchingRadius was explicitly set, use that value
#	Line 118 \| Line 124 \| namespace OpenMD {
124		else
125		mdFileVersion = 0;
126
127	+	// We need the list of simulated atom types to figure out cutoffs
128	+	// as well as long range corrections.
129	+
130	+	set<AtomType*>::iterator i;
131	+	set<AtomType*> atomTypes_;
132	+	atomTypes_ = info_->getSimulatedAtomTypes();
133	+
134		if (simParams_->haveCutoffRadius()) {
135		rCut_ = simParams_->getCutoffRadius();
136		} else {
#	Line 132 \| Line 145 \| namespace OpenMD {
145		rCut_ = 12.0;
146		} else {
147		RealType thisCut;
148	<	set<AtomType*>::iterator i;
136	<	set<AtomType*> atomTypes;
137	<	atomTypes = info_->getSimulatedAtomTypes();
138	<	for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
148	>	for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
149		thisCut = interactionMan_->getSuggestedCutoffRadius((*i));
150		rCut_ = max(thisCut, rCut_);
151		}
#	Line 149 \| Line 159 \| namespace OpenMD {
159		}
160		}
161
162	<	fDecomp_->setUserCutoff(rCut_);
162	>	fDecomp_->setCutoffRadius(rCut_);
163		interactionMan_->setCutoffRadius(rCut_);
164	+	rCutSq_ = rCut_ * rCut_;
165
166		map<string, CutoffMethod> stringToCutoffMethod;
167		stringToCutoffMethod["HARD"] = HARD;
168		stringToCutoffMethod["SWITCHED"] = SWITCHED;
169		stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
170		stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
171	+	stringToCutoffMethod["TAYLOR_SHIFTED"] = TAYLOR_SHIFTED;
172	+	stringToCutoffMethod["EWALD_FULL"] = EWALD_FULL;
173
174		if (simParams_->haveCutoffMethod()) {
175		string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
#	Line 166 \| Line 179 \| namespace OpenMD {
179		sprintf(painCave.errMsg,
180		"ForceManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
181		"\tShould be one of: "
182	<	"HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
182	>	"HARD, SWITCHED, SHIFTED_POTENTIAL, TAYLOR_SHIFTED,\n"
183	>	"\tSHIFTED_FORCE, or EWALD_FULL\n",
184		cutMeth.c_str());
185		painCave.isFatal = 1;
186		painCave.severity = OPENMD_ERROR;
#	Line 210 \| Line 224 \| namespace OpenMD {
224		cutoffMethod_ = SHIFTED_POTENTIAL;
225		} else if (myMethod == "SHIFTED_FORCE") {
226		cutoffMethod_ = SHIFTED_FORCE;
227	+	} else if (myMethod == "TAYLOR_SHIFTED") {
228	+	cutoffMethod_ = TAYLOR_SHIFTED;
229	+	} else if (myMethod == "EWALD_FULL") {
230	+	cutoffMethod_ = EWALD_FULL;
231		}
232
233		if (simParams_->haveSwitchingRadius())
234		rSwitch_ = simParams_->getSwitchingRadius();
235
236	<	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE") {
236	>	if (myMethod == "SHIFTED_POTENTIAL" \|\| myMethod == "SHIFTED_FORCE" \|\|
237	>	myMethod == "TAYLOR_SHIFTED" \|\| myMethod == "EWALD_FULL") {
238		if (simParams_->haveSwitchingRadius()){
239		sprintf(painCave.errMsg,
240		"ForceManager::setupCutoffs : DEPRECATED ERROR MESSAGE\n"
#	Line 250 \| Line 269 \| namespace OpenMD {
269		}
270		}
271		}
253	–	}
254	–	}
255	–
256	–	map<string, CutoffPolicy> stringToCutoffPolicy;
257	–	stringToCutoffPolicy["MIX"] = MIX;
258	–	stringToCutoffPolicy["MAX"] = MAX;
259	–	stringToCutoffPolicy["TRADITIONAL"] = TRADITIONAL;
260	–
261	–	string cutPolicy;
262	–	if (forceFieldOptions_.haveCutoffPolicy()){
263	–	cutPolicy = forceFieldOptions_.getCutoffPolicy();
264	–	}else if (simParams_->haveCutoffPolicy()) {
265	–	cutPolicy = simParams_->getCutoffPolicy();
266	–	}
267	–
268	–	if (!cutPolicy.empty()){
269	–	toUpper(cutPolicy);
270	–	map<string, CutoffPolicy>::iterator i;
271	–	i = stringToCutoffPolicy.find(cutPolicy);
272	–
273	–	if (i == stringToCutoffPolicy.end()) {
274	–	sprintf(painCave.errMsg,
275	–	"ForceManager::setupCutoffs: Could not find chosen cutoffPolicy %s\n"
276	–	"\tShould be one of: "
277	–	"MIX, MAX, or TRADITIONAL\n",
278	–	cutPolicy.c_str());
279	–	painCave.isFatal = 1;
280	–	painCave.severity = OPENMD_ERROR;
281	–	simError();
282	–	} else {
283	–	cutoffPolicy_ = i->second;
272		}
285	–	} else {
286	–	sprintf(painCave.errMsg,
287	–	"ForceManager::setupCutoffs: No value was set for the cutoffPolicy.\n"
288	–	"\tOpenMD will use TRADITIONAL.\n");
289	–	painCave.isFatal = 0;
290	–	painCave.severity = OPENMD_INFO;
291	–	simError();
292	–	cutoffPolicy_ = TRADITIONAL;
273		}
294	–
295	–	fDecomp_->setCutoffPolicy(cutoffPolicy_);
274
275		// create the switching function object:
276
#	Line 370 \| Line 348 \| namespace OpenMD {
348		}
349		switcher_->setSwitchType(sft_);
350		switcher_->setSwitch(rSwitch_, rCut_);
373	–	interactionMan_->setSwitchingRadius(rSwitch_);
351		}
352
376	–
377	–
378	–
353		void ForceManager::initialize() {
354
355		if (!info_->isTopologyDone()) {
#	Line 384 \| Line 358 \| namespace OpenMD {
358		interactionMan_->setSimInfo(info_);
359		interactionMan_->initialize();
360
361	<	// We want to delay the cutoffs until after the interaction
362	<	// manager has set up the atom-atom interactions so that we can
363	<	// query them for suggested cutoff values
361	>	//! We want to delay the cutoffs until after the interaction
362	>	//! manager has set up the atom-atom interactions so that we can
363	>	//! query them for suggested cutoff values
364		setupCutoffs();
365
366		info_->prepareTopology();
#	Line 394 \| Line 368 \| namespace OpenMD {
368		doParticlePot_ = info_->getSimParams()->getOutputParticlePotential();
369		doHeatFlux_ = info_->getSimParams()->getPrintHeatFlux();
370		if (doHeatFlux_) doParticlePot_ = true;
371	+
372	+	doElectricField_ = info_->getSimParams()->getOutputElectricField();
373	+	doSitePotential_ = info_->getSimParams()->getOutputSitePotential();
374
375		}
376
377		ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
378
379	<	// Force fields can set options on how to scale van der Waals and
380	<	// electrostatic interactions for atoms connected via bonds, bends
381	<	// and torsions in this case the topological distance between
382	<	// atoms is:
383	<	// 0 = topologically unconnected
384	<	// 1 = bonded together
385	<	// 2 = connected via a bend
386	<	// 3 = connected via a torsion
379	>	//! Force fields can set options on how to scale van der Waals and
380	>	//! electrostatic interactions for atoms connected via bonds, bends
381	>	//! and torsions in this case the topological distance between
382	>	//! atoms is:
383	>	//! 0 = topologically unconnected
384	>	//! 1 = bonded together
385	>	//! 2 = connected via a bend
386	>	//! 3 = connected via a torsion
387
388		vdwScale_.reserve(4);
389		fill(vdwScale_.begin(), vdwScale_.end(), 0.0);
#	Line 424 \| Line 401 \| namespace OpenMD {
401		electrostaticScale_[2] = fopts.getelectrostatic13scale();
402		electrostaticScale_[3] = fopts.getelectrostatic14scale();
403
404	<	if (info_->getSimParams()->haveElectricField()) {
405	<	ElectricField* eField = new ElectricField(info_);
404	>	if (info_->getSimParams()->haveUniformField()) {
405	>	UniformField* eField = new UniformField(info_);
406		perturbations_.push_back(eField);
407		}
408	<
408	>	if (info_->getSimParams()->haveUniformGradientStrength() \|\|
409	>	info_->getSimParams()->haveUniformGradientDirection1() \|\|
410	>	info_->getSimParams()->haveUniformGradientDirection2() ) {
411	>	UniformGradient* eGrad = new UniformGradient(info_);
412	>	perturbations_.push_back(eGrad);
413	>	}
414	>
415	>	usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
416	>
417		fDecomp_->distributeInitialData();
418	<
418	>
419		initialized_ = true;
420	<
420	>
421		}
422	<
422	>
423		void ForceManager::calcForces() {
424
425		if (!initialized_) initialize();
426	<
426	>
427		preCalculation();
428		shortRangeInteractions();
429		longRangeInteractions();
#	Line 612 \| Line 597 \| namespace OpenMD {
597		// Collect from all nodes. This should eventually be moved into a
598		// SystemDecomposition, but this is a better place than in
599		// Thermo to do the collection.
600	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bondPotential, 1, MPI::REALTYPE,
601	<	MPI::SUM);
602	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bendPotential, 1, MPI::REALTYPE,
603	<	MPI::SUM);
604	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &torsionPotential, 1,
605	<	MPI::REALTYPE, MPI::SUM);
606	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &inversionPotential, 1,
607	<	MPI::REALTYPE, MPI::SUM);
600	>
601	>	MPI_Allreduce(MPI_IN_PLACE, &bondPotential, 1, MPI_REALTYPE,
602	>	MPI_SUM, MPI_COMM_WORLD);
603	>	MPI_Allreduce(MPI_IN_PLACE, &bendPotential, 1, MPI_REALTYPE,
604	>	MPI_SUM, MPI_COMM_WORLD);
605	>	MPI_Allreduce(MPI_IN_PLACE, &torsionPotential, 1,
606	>	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
607	>	MPI_Allreduce(MPI_IN_PLACE, &inversionPotential, 1,
608	>	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
609		#endif
610
611		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
#	Line 637 \| Line 623 \| namespace OpenMD {
623
624		void ForceManager::longRangeInteractions() {
625
640	–
626		Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
627		DataStorage* config = &(curSnapshot->atomData);
628		DataStorage* cgConfig = &(curSnapshot->cgData);
629	+	int jstart, jend;
630
631		//calculate the center of mass of cutoff group
632
#	Line 648 \| Line 634 \| namespace OpenMD {
634		Molecule* mol;
635		Molecule::CutoffGroupIterator ci;
636		CutoffGroup* cg;
637	<
638	<	if(info_->getNCutoffGroups() > 0){
637	>
638	>	if(info_->getNCutoffGroups() != info_->getNAtoms()){
639		for (mol = info_->beginMolecule(mi); mol != NULL;
640		mol = info_->nextMolecule(mi)) {
641		for(cg = mol->beginCutoffGroup(ci); cg != NULL;
#	Line 674 \| Line 660 \| namespace OpenMD {
660		RealType vij;
661		Vector3d fij, fg, f1;
662		tuple3<RealType, RealType, RealType> cuts;
663	<	RealType rCutSq;
663	>	RealType rCut, rCutSq, rListSq;
664		bool in_switching_region;
665		RealType sw, dswdr, swderiv;
666	<	vector<int> atomListColumn, atomListRow, atomListLocal;
666	>	vector<int> atomListColumn, atomListRow;
667		InteractionData idat;
668		SelfData sdat;
669		RealType mf;
670		RealType vpair;
671		RealType dVdFQ1(0.0);
672		RealType dVdFQ2(0.0);
673	<	Vector3d eField1(0.0);
674	<	Vector3d eField2(0.0);
689	<	potVec longRangePotential(0.0);
673	>	potVec longRangePotential(0.0);
674	>	RealType reciprocalPotential(0.0);
675		potVec workPot(0.0);
676		potVec exPot(0.0);
677	+	Vector3d eField1(0.0);
678	+	Vector3d eField2(0.0);
679	+	RealType sPot1(0.0);
680	+	RealType sPot2(0.0);
681	+	bool newAtom1;
682	+
683		vector<int>::iterator ia, jb;
684
685		int loopStart, loopEnd;
686	<
686	>
687	>	idat.rcut = &rCut_;
688		idat.vdwMult = &vdwMult;
689		idat.electroMult = &electroMult;
690		idat.pot = &workPot;
#	Line 703 \| Line 695 \| namespace OpenMD {
695		idat.dVdFQ1 = &dVdFQ1;
696		idat.dVdFQ2 = &dVdFQ2;
697		idat.eField1 = &eField1;
698	<	idat.eField2 = &eField2;
698	>	idat.eField2 = &eField2;
699	>	idat.sPot1 = &sPot1;
700	>	idat.sPot2 = &sPot2;
701		idat.f1 = &f1;
702		idat.sw = &sw;
703		idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false;
704	<	idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE) ? true : false;
704	>	idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE \|\|
705	>	cutoffMethod_ == TAYLOR_SHIFTED) ? true : false;
706		idat.doParticlePot = doParticlePot_;
707	+	idat.doElectricField = doElectricField_;
708	+	idat.doSitePotential = doSitePotential_;
709		sdat.doParticlePot = doParticlePot_;
710
711		loopEnd = PAIR_LOOP;
#	Line 721 \| Line 718 \| namespace OpenMD {
718
719		if (iLoop == loopStart) {
720		bool update_nlist = fDecomp_->checkNeighborList();
721	<	if (update_nlist)
722	<	neighborList = fDecomp_->buildNeighborList();
723	<	}
724	<
725	<	for (vector<pair<int, int> >::iterator it = neighborList.begin();
726	<	it != neighborList.end(); ++it) {
730	<
731	<	cg1 = (*it).first;
732	<	cg2 = (*it).second;
733	<
734	<	cuts = fDecomp_->getGroupCutoffs(cg1, cg2);
735	<
736	<	d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
721	>	if (update_nlist) {
722	>	if (!usePeriodicBoundaryConditions_)
723	>	Mat3x3d bbox = thermo->getBoundingBox();
724	>	fDecomp_->buildNeighborList(neighborList_, point_);
725	>	}
726	>	}
727
728	<	curSnapshot->wrapVector(d_grp);
729	<	rgrpsq = d_grp.lengthSquare();
730	<	rCutSq = cuts.second;
728	>	for (unsigned int cg1 = 0; cg1 < point_.size() - 1; cg1++) {
729	>
730	>	atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
731	>	newAtom1 = true;
732	>
733	>	for (int m2 = point_[cg1]; m2 < point_[cg1+1]; m2++) {
734
735	<	if (rgrpsq < rCutSq) {
743	<	idat.rcut = &cuts.first;
744	<	if (iLoop == PAIR_LOOP) {
745	<	vij = 0.0;
746	<	fij = V3Zero;
747	<	}
735	>	cg2 = neighborList_[m2];
736
737	<	in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
738	<	rgrp);
739	<
740	<	atomListRow = fDecomp_->getAtomsInGroupRow(cg1);
741	<	atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
742	<
743	<	if (doHeatFlux_)
744	<	gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
745	<
746	<	for (ia = atomListRow.begin();
747	<	ia != atomListRow.end(); ++ia) {
748	<	atom1 = (*ia);
749	<
750	<	for (jb = atomListColumn.begin();
763	<	jb != atomListColumn.end(); ++jb) {
764	<	atom2 = (*jb);
765	<
766	<	if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
767	<
768	<	vpair = 0.0;
769	<	workPot = 0.0;
770	<	exPot = 0.0;
771	<	f1 = V3Zero;
772	<	dVdFQ1 = 0.0;
773	<	dVdFQ2 = 0.0;
774	<
775	<	fDecomp_->fillInteractionData(idat, atom1, atom2);
776	<
777	<	topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
778	<	vdwMult = vdwScale_[topoDist];
779	<	electroMult = electrostaticScale_[topoDist];
780	<
781	<	if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
782	<	idat.d = &d_grp;
783	<	idat.r2 = &rgrpsq;
784	<	if (doHeatFlux_)
785	<	vel2 = gvel2;
786	<	} else {
787	<	d = fDecomp_->getInteratomicVector(atom1, atom2);
788	<	curSnapshot->wrapVector( d );
789	<	r2 = d.lengthSquare();
790	<	idat.d = &d;
791	<	idat.r2 = &r2;
792	<	if (doHeatFlux_)
793	<	vel2 = fDecomp_->getAtomVelocityColumn(atom2);
794	<	}
795	<
796	<	r = sqrt( *(idat.r2) );
797	<	idat.rij = &r;
798	<
799	<	if (iLoop == PREPAIR_LOOP) {
800	<	interactionMan_->doPrePair(idat);
801	<	} else {
802	<	interactionMan_->doPair(idat);
803	<	fDecomp_->unpackInteractionData(idat, atom1, atom2);
804	<	vij += vpair;
805	<	fij += f1;
806	<	stressTensor -= outProduct( *(idat.d), f1);
807	<	if (doHeatFlux_)
808	<	fDecomp_->addToHeatFlux((idat.d) dot(f1, vel2));
809	<	}
810	<	}
737	>	d_grp = fDecomp_->getIntergroupVector(cg1, cg2);
738	>
739	>	// already wrapped in the getIntergroupVector call:
740	>	// curSnapshot->wrapVector(d_grp);
741	>	rgrpsq = d_grp.lengthSquare();
742	>
743	>	if (rgrpsq < rCutSq_) {
744	>	if (iLoop == PAIR_LOOP) {
745	>	vij = 0.0;
746	>	fij.zero();
747	>	eField1.zero();
748	>	eField2.zero();
749	>	sPot1 = 0.0;
750	>	sPot2 = 0.0;
751		}
752	<	}
753	<
754	<	if (iLoop == PAIR_LOOP) {
755	<	if (in_switching_region) {
756	<	swderiv = vij * dswdr / rgrp;
757	<	fg = swderiv * d_grp;
758	<	fij += fg;
759	<
760	<	if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
761	<	stressTensor -= outProduct( *(idat.d), fg);
762	<	if (doHeatFlux_)
763	<	fDecomp_->addToHeatFlux((idat.d) dot(fg, vel2));
752	>
753	>	in_switching_region = switcher_->getSwitch(rgrpsq, sw, dswdr,
754	>	rgrp);
755	>
756	>	atomListColumn = fDecomp_->getAtomsInGroupColumn(cg2);
757	>
758	>	if (doHeatFlux_)
759	>	gvel2 = fDecomp_->getGroupVelocityColumn(cg2);
760	>
761	>	for (ia = atomListRow.begin();
762	>	ia != atomListRow.end(); ++ia) {
763	>	atom1 = (*ia);
764	>
765	>	for (jb = atomListColumn.begin();
766	>	jb != atomListColumn.end(); ++jb) {
767	>	atom2 = (*jb);
768
769	<	}
770	<
771	<	for (ia = atomListRow.begin();
772	<	ia != atomListRow.end(); ++ia) {
773	<	atom1 = (*ia);
774	<	mf = fDecomp_->getMassFactorRow(atom1);
775	<	// fg is the force on atom ia due to cutoff group's
776	<	// presence in switching region
777	<	fg = swderiv * d_grp * mf;
778	<	fDecomp_->addForceToAtomRow(atom1, fg);
779	<	if (atomListRow.size() > 1) {
780	<	if (info_->usesAtomicVirial()) {
781	<	// find the distance between the atom
782	<	// and the center of the cutoff group:
783	<	dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
784	<	stressTensor -= outProduct(dag, fg);
769	>	if (!fDecomp_->skipAtomPair(atom1, atom2, cg1, cg2)) {
770	>
771	>	vpair = 0.0;
772	>	workPot = 0.0;
773	>	exPot = 0.0;
774	>	f1.zero();
775	>	dVdFQ1 = 0.0;
776	>	dVdFQ2 = 0.0;
777	>
778	>	fDecomp_->fillInteractionData(idat, atom1, atom2, newAtom1);
779	>
780	>	topoDist = fDecomp_->getTopologicalDistance(atom1, atom2);
781	>	vdwMult = vdwScale_[topoDist];
782	>	electroMult = electrostaticScale_[topoDist];
783	>
784	>	if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
785	>	idat.d = &d_grp;
786	>	idat.r2 = &rgrpsq;
787		if (doHeatFlux_)
788	<	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
788	>	vel2 = gvel2;
789	>	} else {
790	>	d = fDecomp_->getInteratomicVector(atom1, atom2);
791	>	curSnapshot->wrapVector( d );
792	>	r2 = d.lengthSquare();
793	>	idat.d = &d;
794	>	idat.r2 = &r2;
795	>	if (doHeatFlux_)
796	>	vel2 = fDecomp_->getAtomVelocityColumn(atom2);
797		}
798	+
799	+	r = sqrt( *(idat.r2) );
800	+	idat.rij = &r;
801	+
802	+	if (iLoop == PREPAIR_LOOP) {
803	+	interactionMan_->doPrePair(idat);
804	+	} else {
805	+	interactionMan_->doPair(idat);
806	+	fDecomp_->unpackInteractionData(idat, atom1, atom2);
807	+	vij += vpair;
808	+	fij += f1;
809	+	stressTensor -= outProduct( *(idat.d), f1);
810	+	if (doHeatFlux_)
811	+	fDecomp_->addToHeatFlux((idat.d) dot(f1, vel2));
812	+	}
813		}
814		}
815	<	for (jb = atomListColumn.begin();
816	<	jb != atomListColumn.end(); ++jb) {
817	<	atom2 = (*jb);
818	<	mf = fDecomp_->getMassFactorColumn(atom2);
819	<	// fg is the force on atom jb due to cutoff group's
820	<	// presence in switching region
821	<	fg = -swderiv * d_grp * mf;
822	<	fDecomp_->addForceToAtomColumn(atom2, fg);
823	<
824	<	if (atomListColumn.size() > 1) {
825	<	if (info_->usesAtomicVirial()) {
826	<	// find the distance between the atom
827	<	// and the center of the cutoff group:
828	<	dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
829	<	stressTensor -= outProduct(dag, fg);
830	<	if (doHeatFlux_)
831	<	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
815	>	}
816	>
817	>	if (iLoop == PAIR_LOOP) {
818	>	if (in_switching_region) {
819	>	swderiv = vij * dswdr / rgrp;
820	>	fg = swderiv * d_grp;
821	>	fij += fg;
822	>
823	>	if (atomListRow.size() == 1 && atomListColumn.size() == 1) {
824	>	if (!fDecomp_->skipAtomPair(atomListRow[0],
825	>	atomListColumn[0],
826	>	cg1, cg2)) {
827	>	stressTensor -= outProduct( *(idat.d), fg);
828	>	if (doHeatFlux_)
829	>	fDecomp_->addToHeatFlux((idat.d) dot(fg, vel2));
830	>	}
831	>	}
832	>
833	>	for (ia = atomListRow.begin();
834	>	ia != atomListRow.end(); ++ia) {
835	>	atom1 = (*ia);
836	>	mf = fDecomp_->getMassFactorRow(atom1);
837	>	// fg is the force on atom ia due to cutoff group's
838	>	// presence in switching region
839	>	fg = swderiv * d_grp * mf;
840	>	fDecomp_->addForceToAtomRow(atom1, fg);
841	>	if (atomListRow.size() > 1) {
842	>	if (info_->usesAtomicVirial()) {
843	>	// find the distance between the atom
844	>	// and the center of the cutoff group:
845	>	dag = fDecomp_->getAtomToGroupVectorRow(atom1, cg1);
846	>	stressTensor -= outProduct(dag, fg);
847	>	if (doHeatFlux_)
848	>	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
849	>	}
850		}
851		}
852	+	for (jb = atomListColumn.begin();
853	+	jb != atomListColumn.end(); ++jb) {
854	+	atom2 = (*jb);
855	+	mf = fDecomp_->getMassFactorColumn(atom2);
856	+	// fg is the force on atom jb due to cutoff group's
857	+	// presence in switching region
858	+	fg = -swderiv * d_grp * mf;
859	+	fDecomp_->addForceToAtomColumn(atom2, fg);
860	+
861	+	if (atomListColumn.size() > 1) {
862	+	if (info_->usesAtomicVirial()) {
863	+	// find the distance between the atom
864	+	// and the center of the cutoff group:
865	+	dag = fDecomp_->getAtomToGroupVectorColumn(atom2, cg2);
866	+	stressTensor -= outProduct(dag, fg);
867	+	if (doHeatFlux_)
868	+	fDecomp_->addToHeatFlux( dag * dot(fg, vel2));
869	+	}
870	+	}
871	+	}
872		}
873	+	//if (!info_->usesAtomicVirial()) {
874	+	// stressTensor -= outProduct(d_grp, fij);
875	+	// if (doHeatFlux_)
876	+	// fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
877	+	//}
878		}
867	–	//if (!info_->usesAtomicVirial()) {
868	–	// stressTensor -= outProduct(d_grp, fij);
869	–	// if (doHeatFlux_)
870	–	// fDecomp_->addToHeatFlux( d_grp * dot(fij, vel2));
871	–	//}
879		}
880		}
881	+	newAtom1 = false;
882		}
883	<
883	>
884		if (iLoop == PREPAIR_LOOP) {
885		if (info_->requiresPrepair()) {
886	<
886	>
887		fDecomp_->collectIntermediateData();
888	<
888	>
889		for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
890		fDecomp_->fillSelfData(sdat, atom1);
891		interactionMan_->doPreForce(sdat);
892		}
893	<
893	>
894		fDecomp_->distributeIntermediateData();
895	<
895	>
896		}
897		}
898		}
899	<
899	>
900		// collects pairwise information
901		fDecomp_->collectData();
902	+	if (cutoffMethod_ == EWALD_FULL) {
903	+	interactionMan_->doReciprocalSpaceSum(reciprocalPotential);
904	+
905	+	curSnapshot->setReciprocalPotential(reciprocalPotential);
906	+	}
907
908		if (info_->requiresSelfCorrection()) {
909		for (unsigned int atom1 = 0; atom1 < info_->getNAtoms(); atom1++) {
#	Line 908 \| Line 921 \| namespace OpenMD {
921		curSnapshot->setLongRangePotential(longRangePotential);
922
923		curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) +
924	<	*(fDecomp_->getExcludedPotential()));
924	>	*(fDecomp_->getExcludedPotential()));
925
926		}
927
915	–
928		void ForceManager::postCalculation() {
929
930		vector<Perturbation*>::iterator pi;
#	Line 938 \| Line 950 \| namespace OpenMD {
950		}
951
952		#ifdef IS_MPI
953	<	MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, stressTensor.getArrayPointer(), 9,
954	<	MPI::REALTYPE, MPI::SUM);
953	>	MPI_Allreduce(MPI_IN_PLACE, stressTensor.getArrayPointer(), 9,
954	>	MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
955		#endif
956		curSnapshot->setStressTensor(stressTensor);
957
958	+	if (info_->getSimParams()->getUseLongRangeCorrections()) {
959	+	/*
960	+	RealType vol = curSnapshot->getVolume();
961	+	RealType Elrc(0.0);
962	+	RealType Wlrc(0.0);
963	+
964	+	set<AtomType*>::iterator i;
965	+	set<AtomType*>::iterator j;
966	+
967	+	RealType n_i, n_j;
968	+	RealType rho_i, rho_j;
969	+	pair<RealType, RealType> LRI;
970	+
971	+	for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) {
972	+	n_i = RealType(info_->getGlobalCountOfType(*i));
973	+	rho_i = n_i / vol;
974	+	for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) {
975	+	n_j = RealType(info_->getGlobalCountOfType(*j));
976	+	rho_j = n_j / vol;
977	+
978	+	LRI = interactionMan_->getLongRangeIntegrals( (i), (j) );
979	+
980	+	Elrc += n_i * rho_j * LRI.first;
981	+	Wlrc -= rho_i * rho_j * LRI.second;
982	+	}
983	+	}
984	+	Elrc = 2.0 NumericConstant::PI;
985	+	Wlrc = 2.0 NumericConstant::PI;
986	+
987	+	RealType lrp = curSnapshot->getLongRangePotential();
988	+	curSnapshot->setLongRangePotential(lrp + Elrc);
989	+	stressTensor += Wlrc * SquareMatrix3<RealType>::identity();
990	+	curSnapshot->setStressTensor(stressTensor);
991	+	*/
992	+
993	+	}
994		}
995	<	} //end namespace OpenMD
995	>	}

Diff Legend

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

Comparing trunk/src/brains/ForceManager.cpp (file contents): Revision 1788 by gezelter, Wed Aug 29 20:17:07 2012 UTC vs. Revision 2057 by gezelter, Tue Mar 3 15:22:26 2015 UTC

Diff Legend

Comparing trunk/src/brains/ForceManager.cpp (file contents):
Revision 1788 by gezelter, Wed Aug 29 20:17:07 2012 UTC vs.
Revision 2057 by gezelter, Tue Mar 3 15:22:26 2015 UTC