[ViewVC] Diff of: OpenMD/branches/development/src/nonbonded/Electrostatic.cpp

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
Revision 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC

#	Line 34 \| Line 34
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).
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		*/
#	Line 47 \| Line 47
47		#include "utils/simError.h"
48		#include "types/NonBondedInteractionType.hpp"
49		#include "types/DirectionalAtomType.hpp"
50	+	#include "io/Globals.hpp"
51
51	–
52		namespace OpenMD {
53
54		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
55	<	forceField_(NULL) {}
55	>	forceField_(NULL), info_(NULL),
56	>	haveCutoffRadius_(false),
57	>	haveDampingAlpha_(false),
58	>	haveDielectric_(false),
59	>	haveElectroSpline_(false)
60	>	{}
61
62		void Electrostatic::initialize() {
63	+
64	+	Globals* simParams_ = info_->getSimParams();
65	+
66	+	summationMap_["HARD"] = esm_HARD;
67	+	summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
68	+	summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
69	+	summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
70	+	summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD;
71	+	summationMap_["EWALD_FULL"] = esm_EWALD_FULL;
72	+	summationMap_["EWALD_PME"] = esm_EWALD_PME;
73	+	summationMap_["EWALD_SPME"] = esm_EWALD_SPME;
74	+	screeningMap_["DAMPED"] = DAMPED;
75	+	screeningMap_["UNDAMPED"] = UNDAMPED;
76	+
77		// these prefactors convert the multipole interactions into kcal / mol
78		// all were computed assuming distances are measured in angstroms
79		// Charge-Charge, assuming charges are measured in electrons
#	Line 79 \| Line 98 \| namespace OpenMD {
98
99		// variables to handle different summation methods for long-range
100		// electrostatics:
101	<	summationMethod_ = NONE;
101	>	summationMethod_ = esm_HARD;
102		screeningMethod_ = UNDAMPED;
103		dielectric_ = 1.0;
104		one_third_ = 1.0 / 3.0;
86	–	haveDefaultCutoff_ = false;
87	–	haveDampingAlpha_ = false;
88	–	haveDielectric_ = false;
89	–	haveElectroSpline_ = false;
105
106	+	// check the summation method:
107	+	if (simParams_->haveElectrostaticSummationMethod()) {
108	+	string myMethod = simParams_->getElectrostaticSummationMethod();
109	+	toUpper(myMethod);
110	+	map<string, ElectrostaticSummationMethod>::iterator i;
111	+	i = summationMap_.find(myMethod);
112	+	if ( i != summationMap_.end() ) {
113	+	summationMethod_ = (*i).second;
114	+	} else {
115	+	// throw error
116	+	sprintf( painCave.errMsg,
117	+	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
118	+	"\t(Input file specified %s .)\n"
119	+	"\telectrostaticSummationMethod must be one of: \"none\",\n"
120	+	"\t\"shifted_potential\", \"shifted_force\", or \n"
121	+	"\t\"reaction_field\".\n", myMethod.c_str() );
122	+	painCave.isFatal = 1;
123	+	simError();
124	+	}
125	+	} else {
126	+	// set ElectrostaticSummationMethod to the cutoffMethod:
127	+	if (simParams_->haveCutoffMethod()){
128	+	string myMethod = simParams_->getCutoffMethod();
129	+	toUpper(myMethod);
130	+	map<string, ElectrostaticSummationMethod>::iterator i;
131	+	i = summationMap_.find(myMethod);
132	+	if ( i != summationMap_.end() ) {
133	+	summationMethod_ = (*i).second;
134	+	}
135	+	}
136	+	}
137	+
138	+	if (summationMethod_ == esm_REACTION_FIELD) {
139	+	if (!simParams_->haveDielectric()) {
140	+	// throw warning
141	+	sprintf( painCave.errMsg,
142	+	"SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
143	+	"\tA default value of %f will be used for the dielectric.\n", dielectric_);
144	+	painCave.isFatal = 0;
145	+	painCave.severity = OPENMD_INFO;
146	+	simError();
147	+	} else {
148	+	dielectric_ = simParams_->getDielectric();
149	+	}
150	+	haveDielectric_ = true;
151	+	}
152	+
153	+	if (simParams_->haveElectrostaticScreeningMethod()) {
154	+	string myScreen = simParams_->getElectrostaticScreeningMethod();
155	+	toUpper(myScreen);
156	+	map<string, ElectrostaticScreeningMethod>::iterator i;
157	+	i = screeningMap_.find(myScreen);
158	+	if ( i != screeningMap_.end()) {
159	+	screeningMethod_ = (*i).second;
160	+	} else {
161	+	sprintf( painCave.errMsg,
162	+	"SimInfo error: Unknown electrostaticScreeningMethod.\n"
163	+	"\t(Input file specified %s .)\n"
164	+	"\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
165	+	"or \"damped\".\n", myScreen.c_str() );
166	+	painCave.isFatal = 1;
167	+	simError();
168	+	}
169	+	}
170	+
171	+	// check to make sure a cutoff value has been set:
172	+	if (!haveCutoffRadius_) {
173	+	sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
174	+	"Cutoff value!\n");
175	+	painCave.severity = OPENMD_ERROR;
176	+	painCave.isFatal = 1;
177	+	simError();
178	+	}
179	+
180	+	if (screeningMethod_ == DAMPED) {
181	+	if (!simParams_->haveDampingAlpha()) {
182	+	// first set a cutoff dependent alpha value
183	+	// we assume alpha depends linearly with rcut from 0 to 20.5 ang
184	+	dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
185	+	if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
186	+
187	+	// throw warning
188	+	sprintf( painCave.errMsg,
189	+	"Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
190	+	"\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
191	+	dampingAlpha_, cutoffRadius_);
192	+	painCave.severity = OPENMD_INFO;
193	+	painCave.isFatal = 0;
194	+	simError();
195	+	} else {
196	+	dampingAlpha_ = simParams_->getDampingAlpha();
197	+	}
198	+	haveDampingAlpha_ = true;
199	+	}
200	+
201		// find all of the Electrostatic atom Types:
202		ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
203		ForceField::AtomTypeContainer::MapTypeIterator i;
204		AtomType* at;
205	<
205	>
206		for (at = atomTypes->beginType(i); at != NULL;
207		at = atomTypes->nextType(i)) {
208
#	Line 100 \| Line 210 \| namespace OpenMD {
210		addType(at);
211		}
212
103	–	// check to make sure a cutoff value has been set:
104	–	if (!haveDefaultCutoff_) {
105	–	sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
106	–	"Cutoff value!\n");
107	–	painCave.severity = OPENMD_ERROR;
108	–	painCave.isFatal = 1;
109	–	simError();
110	–	}
213
214	<	defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
215	<	rcuti_ = 1.0 / defaultCutoff_;
214	>	cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
215	>	rcuti_ = 1.0 / cutoffRadius_;
216		rcuti2_ = rcuti_ * rcuti_;
217		rcuti3_ = rcuti2_ * rcuti_;
218		rcuti4_ = rcuti2_ * rcuti2_;
219
220		if (screeningMethod_ == DAMPED) {
221	<	if (!haveDampingAlpha_) {
120	<	sprintf( painCave.errMsg, "Electrostatic::initialize has no "
121	<	"DampingAlpha value!\n");
122	<	painCave.severity = OPENMD_ERROR;
123	<	painCave.isFatal = 1;
124	<	simError();
125	<	}
126	<
221	>
222		alpha2_ = dampingAlpha_ * dampingAlpha_;
223		alpha4_ = alpha2_ * alpha2_;
224		alpha6_ = alpha4_ * alpha2_;
225		alpha8_ = alpha4_ * alpha4_;
226
227	<	constEXP_ = exp(-alpha2_ * defaultCutoff2_);
227	>	constEXP_ = exp(-alpha2_ * cutoffRadius2_);
228		invRootPi_ = 0.56418958354775628695;
229		alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
230
231	<	c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
231	>	c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
232		c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
233		c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
234		c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
#	Line 148 \| Line 243 \| namespace OpenMD {
243		c6c_ = 9.0 * c5c_ * rcuti2_;
244		}
245
246	<	if (summationMethod_ == REACTION_FIELD) {
247	<	if (haveDielectric_) {
248	<	preRF_ = (dielectric_ - 1.0) /
249	<	((2.0 * dielectric_ + 1.0) * defaultCutoff2_ * defaultCutoff_);
155	<	preRF2_ = 2.0 * preRF_;
156	<	} else {
157	<	sprintf( painCave.errMsg, "Electrostatic::initialize has no Dielectric"
158	<	" value!\n");
159	<	painCave.severity = OPENMD_ERROR;
160	<	painCave.isFatal = 1;
161	<	simError();
162	<	}
246	>	if (summationMethod_ == esm_REACTION_FIELD) {
247	>	preRF_ = (dielectric_ - 1.0) /
248	>	((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
249	>	preRF2_ = 2.0 * preRF_;
250		}
251	<
252	<	RealType dx = defaultCutoff_ / RealType(np_ - 1);
251	>
252	>	RealType dx = cutoffRadius_ / RealType(np_ - 1);
253		RealType rval;
254		vector<RealType> rvals;
255		vector<RealType> yvals;
#	Line 283 \| Line 370 \| namespace OpenMD {
370		simError();
371		}
372
373	+	// Quadrupoles in OpenMD are set as the diagonal elements
374	+	// of the diagonalized traceless quadrupole moment tensor.
375	+	// The column vectors of the unitary matrix that diagonalizes
376	+	// the quadrupole moment tensor become the eFrame (or the
377	+	// electrostatic version of the body-fixed frame.
378	+
379		Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
380		if (v3dData == NULL) {
381		sprintf( painCave.errMsg,
#	Line 315 \| Line 408 \| namespace OpenMD {
408		return;
409		}
410
411	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
412	<	RealType theRSW ) {
413	<	defaultCutoff_ = theECR;
414	<	rrf_ = defaultCutoff_;
322	<	rt_ = theRSW;
323	<	haveDefaultCutoff_ = true;
411	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
412	>	cutoffRadius_ = rCut;
413	>	rrf_ = cutoffRadius_;
414	>	haveCutoffRadius_ = true;
415		}
416	+
417	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
418	+	rt_ = rSwitch;
419	+	}
420		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
421		summationMethod_ = esm;
422		}
#	Line 337 \| Line 432 \| namespace OpenMD {
432		haveDielectric_ = true;
433		}
434
435	<	void Electrostatic::calcForce(InteractionData idat) {
435	>	void Electrostatic::calcForce(InteractionData &idat) {
436
437		// utility variables. Should clean these up and use the Vector3d and
438		// Mat3x3d to replace as many as we can in future versions:
#	Line 351 \| Line 446 \| namespace OpenMD {
446		RealType ct_i, ct_j, ct_ij, a1;
447		RealType riji, ri, ri2, ri3, ri4;
448		RealType pref, vterm, epot, dudr;
449	+	RealType vpair(0.0);
450		RealType scale, sc2;
451		RealType pot_term, preVal, rfVal;
452		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
453		RealType preSw, preSwSc;
454		RealType c1, c2, c3, c4;
455	<	RealType erfcVal, derfcVal;
455	>	RealType erfcVal(1.0), derfcVal(0.0);
456		RealType BigR;
457
458		Vector3d Q_i, Q_j;
#	Line 367 \| Line 463 \| namespace OpenMD {
463		Vector3d rhatdot2, rhatc4;
464		Vector3d dVdr;
465
466	+	// variables for indirect (reaction field) interactions for excluded pairs:
467	+	RealType indirect_Pot(0.0);
468	+	RealType indirect_vpair(0.0);
469	+	Vector3d indirect_dVdr(V3Zero);
470	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
471	+
472		pair<RealType, RealType> res;
473
474		if (!initialized_) initialize();
475
476	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
477	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
476	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
477	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
478
479		// some variables we'll need independent of electrostatic type:
480
481	<	riji = 1.0 / idat.rij;
482	<	Vector3d rhat = idat.d * riji;
481	>	riji = 1.0 / *(idat.rij) ;
482	>	Vector3d rhat = (idat.d) riji;
483
484		// logicals
485
#	Line 391 \| Line 493 \| namespace OpenMD {
493		bool j_is_SplitDipole = data2.is_SplitDipole;
494		bool j_is_Quadrupole = data2.is_Quadrupole;
495
496	<	if (i_is_Charge)
496	>	if (i_is_Charge) {
497		q_i = data1.charge;
498	+	if (idat.excluded) {
499	+	*(idat.skippedCharge2) += q_i;
500	+	}
501	+	}
502
503		if (i_is_Dipole) {
504		mu_i = data1.dipole_moment;
505	<	uz_i = idat.eFrame1.getColumn(2);
505	>	uz_i = idat.eFrame1->getColumn(2);
506
507		ct_i = dot(uz_i, rhat);
508
#	Line 412 \| Line 518 \| namespace OpenMD {
518		qyy_i = Q_i.y();
519		qzz_i = Q_i.z();
520
521	<	ux_i = idat.eFrame1.getColumn(0);
522	<	uy_i = idat.eFrame1.getColumn(1);
523	<	uz_i = idat.eFrame1.getColumn(2);
521	>	ux_i = idat.eFrame1->getColumn(0);
522	>	uy_i = idat.eFrame1->getColumn(1);
523	>	uz_i = idat.eFrame1->getColumn(2);
524
525		cx_i = dot(ux_i, rhat);
526		cy_i = dot(uy_i, rhat);
#	Line 425 \| Line 531 \| namespace OpenMD {
531		duduz_i = V3Zero;
532		}
533
534	<	if (j_is_Charge)
534	>	if (j_is_Charge) {
535		q_j = data2.charge;
536	+	if (idat.excluded) {
537	+	*(idat.skippedCharge1) += q_j;
538	+	}
539	+	}
540
541	+
542		if (j_is_Dipole) {
543		mu_j = data2.dipole_moment;
544	<	uz_j = idat.eFrame2.getColumn(2);
544	>	uz_j = idat.eFrame2->getColumn(2);
545
546		ct_j = dot(uz_j, rhat);
547
#	Line 446 \| Line 557 \| namespace OpenMD {
557		qyy_j = Q_j.y();
558		qzz_j = Q_j.z();
559
560	<	ux_j = idat.eFrame2.getColumn(0);
561	<	uy_j = idat.eFrame2.getColumn(1);
562	<	uz_j = idat.eFrame2.getColumn(2);
560	>	ux_j = idat.eFrame2->getColumn(0);
561	>	uy_j = idat.eFrame2->getColumn(1);
562	>	uz_j = idat.eFrame2->getColumn(2);
563
564		cx_j = dot(ux_j, rhat);
565		cy_j = dot(uy_j, rhat);
#	Line 467 \| Line 578 \| namespace OpenMD {
578		if (j_is_Charge) {
579		if (screeningMethod_ == DAMPED) {
580		// assemble the damping variables
581	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
581	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
582		erfcVal = res.first;
583		derfcVal = res.second;
584		c1 = erfcVal * riji;
#	Line 477 \| Line 588 \| namespace OpenMD {
588		c2 = c1 * riji;
589		}
590
591	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
591	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
592
593	<	if (summationMethod_ == SHIFTED_POTENTIAL) {
593	>	if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
594		vterm = preVal * (c1 - c1c_);
595	<	dudr = -idat.sw * preVal * c2;
595	>	dudr = - (idat.sw) preVal * c2;
596
597	<	} else if (summationMethod_ == SHIFTED_FORCE) {
598	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
599	<	dudr = idat.sw * preVal * (c2c_ - c2);
597	>	} else if (summationMethod_ == esm_SHIFTED_FORCE) {
598	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
599	>	dudr = (idat.sw) preVal * (c2c_ - c2);
600
601	<	} else if (summationMethod_ == REACTION_FIELD) {
602	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
601	>	} else if (summationMethod_ == esm_REACTION_FIELD) {
602	>	rfVal = preRF_ * (idat.rij) *(idat.rij);
603	>
604		vterm = preVal * ( riji + rfVal );
605	<	dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
605	>	dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
606	>
607	>	// if this is an excluded pair, there are still indirect
608	>	// interactions via the reaction field we must worry about:
609
610	+	if (idat.excluded) {
611	+	indirect_vpair += preVal * rfVal;
612	+	indirect_Pot += (idat.sw) preVal * rfVal;
613	+	indirect_dVdr += (idat.sw) preVal * 2.0 * rfVal * riji * rhat;
614	+	}
615	+
616		} else {
496	–	vterm = preVal * riji * erfcVal;
617
618	<	dudr = - idat.sw * preVal * c2;
618	>	vterm = preVal * riji * erfcVal;
619	>	dudr = - (idat.sw) preVal * c2;
620
621		}
501	–
502	–	idat.vpair += vterm;
503	–	epot += idat.sw * vterm;
622
623	<	dVdr += dudr * rhat;
623	>	vpair += vterm;
624	>	epot += (idat.sw) vterm;
625	>	dVdr += dudr * rhat;
626		}
627
628		if (j_is_Dipole) {
629		// pref is used by all the possible methods
630	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
631	<	preSw = idat.sw * pref;
630	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
631	>	preSw = (idat.sw) pref;
632
633	<	if (summationMethod_ == REACTION_FIELD) {
633	>	if (summationMethod_ == esm_REACTION_FIELD) {
634		ri2 = riji * riji;
635		ri3 = ri2 * riji;
636
637	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
638	<	idat.vpair += vterm;
639	<	epot += idat.sw * vterm;
637	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
638	>	vpair += vterm;
639	>	epot += (idat.sw) vterm;
640
641		dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
642	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
642	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
643
644	+	// Even if we excluded this pair from direct interactions,
645	+	// we still have the reaction-field-mediated charge-dipole
646	+	// interaction:
647	+
648	+	if (idat.excluded) {
649	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
650	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
651	+	indirect_dVdr += preSw * preRF2_ * uz_j;
652	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
653	+	}
654	+
655		} else {
656		// determine the inverse r used if we have split dipoles
657		if (j_is_SplitDipole) {
658	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
658	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
659		ri = 1.0 / BigR;
660	<	scale = idat.rij * ri;
660	>	scale = (idat.rij) ri;
661		} else {
662		ri = riji;
663		scale = 1.0;
#	Line 536 \| Line 667 \| namespace OpenMD {
667
668		if (screeningMethod_ == DAMPED) {
669		// assemble the damping variables
670	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
670	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
671		erfcVal = res.first;
672		derfcVal = res.second;
673		c1 = erfcVal * ri;
#	Line 553 \| Line 684 \| namespace OpenMD {
684		// calculate the potential
685		pot_term = scale * c2;
686		vterm = -pref * ct_j * pot_term;
687	<	idat.vpair += vterm;
688	<	epot += idat.sw * vterm;
687	>	vpair += vterm;
688	>	epot += (idat.sw) vterm;
689
690		// calculate derivatives for forces and torques
691
#	Line 569 \| Line 700 \| namespace OpenMD {
700		cx2 = cx_j * cx_j;
701		cy2 = cy_j * cy_j;
702		cz2 = cz_j * cz_j;
703	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
703	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
704
705		if (screeningMethod_ == DAMPED) {
706		// assemble the damping variables
707	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
707	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
708		erfcVal = res.first;
709		derfcVal = res.second;
710		c1 = erfcVal * riji;
#	Line 588 \| Line 719 \| namespace OpenMD {
719		}
720
721		// precompute variables for convenience
722	<	preSw = idat.sw * pref;
722	>	preSw = (idat.sw) pref;
723		c2ri = c2 * riji;
724		c3ri = c3 * riji;
725	<	c4rij = c4 * idat.rij;
725	>	c4rij = c4 * *(idat.rij) ;
726		rhatdot2 = 2.0 * rhat * c3;
727		rhatc4 = rhat * c4rij;
728
#	Line 600 \| Line 731 \| namespace OpenMD {
731		qyy_j * (cy2*c3 - c2ri) +
732		qzz_j * (cz2*c3 - c2ri) );
733		vterm = pref * pot_term;
734	<	idat.vpair += vterm;
735	<	epot += idat.sw * vterm;
734	>	vpair += vterm;
735	>	epot += (idat.sw) vterm;
736
737		// calculate derivatives for the forces and torques
738
#	Line 619 \| Line 750 \| namespace OpenMD {
750
751		if (j_is_Charge) {
752		// variables used by all the methods
753	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
754	<	preSw = idat.sw * pref;
753	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
754	>	preSw = (idat.sw) pref;
755
756	<	if (summationMethod_ == REACTION_FIELD) {
756	>	if (summationMethod_ == esm_REACTION_FIELD) {
757
758		ri2 = riji * riji;
759		ri3 = ri2 * riji;
760
761	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
762	<	idat.vpair += vterm;
763	<	epot += idat.sw * vterm;
761	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
762	>	vpair += vterm;
763	>	epot += (idat.sw) vterm;
764
765		dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
766
767	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
767	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
768	>
769	>	// Even if we excluded this pair from direct interactions,
770	>	// we still have the reaction-field-mediated charge-dipole
771	>	// interaction:
772	>
773	>	if (idat.excluded) {
774	>	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
775	>	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
776	>	indirect_dVdr += -preSw * preRF2_ * uz_i;
777	>	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
778	>	}
779
780		} else {
781
782		// determine inverse r if we are using split dipoles
783		if (i_is_SplitDipole) {
784	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
784	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
785		ri = 1.0 / BigR;
786	<	scale = idat.rij * ri;
786	>	scale = (idat.rij) ri;
787		} else {
788		ri = riji;
789		scale = 1.0;
#	Line 651 \| Line 793 \| namespace OpenMD {
793
794		if (screeningMethod_ == DAMPED) {
795		// assemble the damping variables
796	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
796	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
797		erfcVal = res.first;
798		derfcVal = res.second;
799		c1 = erfcVal * ri;
#	Line 668 \| Line 810 \| namespace OpenMD {
810		// calculate the potential
811		pot_term = c2 * scale;
812		vterm = pref * ct_i * pot_term;
813	<	idat.vpair += vterm;
814	<	epot += idat.sw * vterm;
813	>	vpair += vterm;
814	>	epot += (idat.sw) vterm;
815
816		// calculate derivatives for the forces and torques
817		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
#	Line 681 \| Line 823 \| namespace OpenMD {
823		// variables used by all methods
824		ct_ij = dot(uz_i, uz_j);
825
826	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
827	<	preSw = idat.sw * pref;
826	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
827	>	preSw = (idat.sw) pref;
828
829	<	if (summationMethod_ == REACTION_FIELD) {
829	>	if (summationMethod_ == esm_REACTION_FIELD) {
830		ri2 = riji * riji;
831		ri3 = ri2 * riji;
832		ri4 = ri2 * ri2;
833
834		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
835		preRF2_ * ct_ij );
836	<	idat.vpair += vterm;
837	<	epot += idat.sw * vterm;
836	>	vpair += vterm;
837	>	epot += (idat.sw) vterm;
838
839		a1 = 5.0 * ct_i * ct_j - ct_ij;
840
#	Line 701 \| Line 843 \| namespace OpenMD {
843		duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
844		duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
845
846	+	if (idat.excluded) {
847	+	indirect_vpair += - pref * preRF2_ * ct_ij;
848	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
849	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
850	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
851	+	}
852	+
853		} else {
854
855		if (i_is_SplitDipole) {
856		if (j_is_SplitDipole) {
857	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
857	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
858		} else {
859	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
859	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
860		}
861		ri = 1.0 / BigR;
862	<	scale = idat.rij * ri;
862	>	scale = (idat.rij) ri;
863		} else {
864		if (j_is_SplitDipole) {
865	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
865	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
866		ri = 1.0 / BigR;
867	<	scale = idat.rij * ri;
867	>	scale = (idat.rij) ri;
868		} else {
869		ri = riji;
870		scale = 1.0;
#	Line 723 \| Line 872 \| namespace OpenMD {
872		}
873		if (screeningMethod_ == DAMPED) {
874		// assemble damping variables
875	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
875	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
876		erfcVal = res.first;
877		derfcVal = res.second;
878		c1 = erfcVal * ri;
#	Line 745 \| Line 894 \| namespace OpenMD {
894		preSwSc = preSw * scale;
895		c2ri = c2 * ri;
896		c3ri = c3 * ri;
897	<	c4rij = c4 * idat.rij;
897	>	c4rij = c4 * *(idat.rij) ;
898
899		// calculate the potential
900		pot_term = (ct_ij * c2ri - ctidotj * c3);
901		vterm = pref * pot_term;
902	<	idat.vpair += vterm;
903	<	epot += idat.sw * vterm;
902	>	vpair += vterm;
903	>	epot += (idat.sw) vterm;
904
905		// calculate derivatives for the forces and torques
906		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 770 \| Line 919 \| namespace OpenMD {
919		cy2 = cy_i * cy_i;
920		cz2 = cz_i * cz_i;
921
922	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
922	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
923
924		if (screeningMethod_ == DAMPED) {
925		// assemble the damping variables
926	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
926	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
927		erfcVal = res.first;
928		derfcVal = res.second;
929		c1 = erfcVal * riji;
#	Line 789 \| Line 938 \| namespace OpenMD {
938		}
939
940		// precompute some variables for convenience
941	<	preSw = idat.sw * pref;
941	>	preSw = (idat.sw) pref;
942		c2ri = c2 * riji;
943		c3ri = c3 * riji;
944	<	c4rij = c4 * idat.rij;
944	>	c4rij = c4 * *(idat.rij) ;
945		rhatdot2 = 2.0 * rhat * c3;
946		rhatc4 = rhat * c4rij;
947
#	Line 802 \| Line 951 \| namespace OpenMD {
951		qzz_i * (cz2 * c3 - c2ri) );
952
953		vterm = pref * pot_term;
954	<	idat.vpair += vterm;
955	<	epot += idat.sw * vterm;
954	>	vpair += vterm;
955	>	epot += (idat.sw) vterm;
956
957		// calculate the derivatives for the forces and torques
958
#	Line 817 \| Line 966 \| namespace OpenMD {
966		}
967		}
968
820	–	idat.pot += epot;
821	–	idat.f1 += dVdr;
969
970	<	if (i_is_Dipole \|\| i_is_Quadrupole)
971	<	idat.t1 -= cross(uz_i, duduz_i);
972	<	if (i_is_Quadrupole) {
973	<	idat.t1 -= cross(ux_i, dudux_i);
974	<	idat.t1 -= cross(uy_i, duduy_i);
975	<	}
970	>	if (!idat.excluded) {
971	>	*(idat.vpair) += vpair;
972	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
973	>	*(idat.f1) += dVdr;
974	>
975	>	if (i_is_Dipole \|\| i_is_Quadrupole)
976	>	*(idat.t1) -= cross(uz_i, duduz_i);
977	>	if (i_is_Quadrupole) {
978	>	*(idat.t1) -= cross(ux_i, dudux_i);
979	>	*(idat.t1) -= cross(uy_i, duduy_i);
980	>	}
981	>
982	>	if (j_is_Dipole \|\| j_is_Quadrupole)
983	>	*(idat.t2) -= cross(uz_j, duduz_j);
984	>	if (j_is_Quadrupole) {
985	>	*(idat.t2) -= cross(uz_j, dudux_j);
986	>	*(idat.t2) -= cross(uz_j, duduy_j);
987	>	}
988
989	<	if (j_is_Dipole \|\| j_is_Quadrupole)
990	<	idat.t2 -= cross(uz_j, duduz_j);
991	<	if (j_is_Quadrupole) {
992	<	idat.t2 -= cross(uz_j, dudux_j);
993	<	idat.t2 -= cross(uz_j, duduy_j);
989	>	} else {
990	>
991	>	// only accumulate the forces and torques resulting from the
992	>	// indirect reaction field terms.
993	>	*(idat.vpair) += indirect_vpair;
994	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
995	>	*(idat.f1) += indirect_dVdr;
996	>
997	>	if (i_is_Dipole)
998	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
999	>	if (j_is_Dipole)
1000	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1001		}
1002
1003	+
1004		return;
1005		}
839	–
840	–	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
841	–
842	–	if (!initialized_) initialize();
1006
1007	<	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
845	<	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
846	<
847	<	// logicals
848	<
849	<	bool i_is_Charge = data1.is_Charge;
850	<	bool i_is_Dipole = data1.is_Dipole;
851	<
852	<	bool j_is_Charge = data2.is_Charge;
853	<	bool j_is_Dipole = data2.is_Dipole;
854	<
855	<	RealType q_i, q_j;
856	<
857	<	// The skippedCharge computation is needed by the real-space cutoff methods
858	<	// (i.e. shifted force and shifted potential)
859	<
860	<	if (i_is_Charge) {
861	<	q_i = data1.charge;
862	<	skdat.skippedCharge2 += q_i;
863	<	}
864	<
865	<	if (j_is_Charge) {
866	<	q_j = data2.charge;
867	<	skdat.skippedCharge1 += q_j;
868	<	}
869	<
870	<	// the rest of this function should only be necessary for reaction field.
871	<
872	<	if (summationMethod_ == REACTION_FIELD) {
873	<	RealType riji, ri2, ri3;
874	<	RealType q_i, mu_i, ct_i;
875	<	RealType q_j, mu_j, ct_j;
876	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
877	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
878	<
879	<	// some variables we'll need independent of electrostatic type:
880	<
881	<	riji = 1.0 / skdat.rij;
882	<	rhat = skdat.d * riji;
883	<
884	<	if (i_is_Dipole) {
885	<	mu_i = data1.dipole_moment;
886	<	uz_i = skdat.eFrame1.getColumn(2);
887	<	ct_i = dot(uz_i, rhat);
888	<	duduz_i = V3Zero;
889	<	}
890	<
891	<	if (j_is_Dipole) {
892	<	mu_j = data2.dipole_moment;
893	<	uz_j = skdat.eFrame2.getColumn(2);
894	<	ct_j = dot(uz_j, rhat);
895	<	duduz_j = V3Zero;
896	<	}
897	<
898	<	if (i_is_Charge) {
899	<	if (j_is_Charge) {
900	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
901	<	rfVal = preRF_ * skdat.rij * skdat.rij;
902	<	vterm = preVal * rfVal;
903	<	myPot += skdat.sw * vterm;
904	<	dudr = skdat.sw * preVal * 2.0 * rfVal * riji;
905	<	dVdr += dudr * rhat;
906	<	}
907	<
908	<	if (j_is_Dipole) {
909	<	ri2 = riji * riji;
910	<	ri3 = ri2 * riji;
911	<	pref = skdat.electroMult * pre12_ * q_i * mu_j;
912	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
913	<	myPot += skdat.sw * vterm;
914	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
915	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
916	<	}
917	<	}
918	<	if (i_is_Dipole) {
919	<	if (j_is_Charge) {
920	<	ri2 = riji * riji;
921	<	ri3 = ri2 * riji;
922	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
923	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
924	<	myPot += skdat.sw * vterm;
925	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
926	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
927	<	}
928	<	}
929	<
930	<	// accumulate the forces and torques resulting from the self term
931	<	skdat.pot += myPot;
932	<	skdat.f1 += dVdr;
933	<
934	<	if (i_is_Dipole)
935	<	skdat.t1 -= cross(uz_i, duduz_i);
936	<	if (j_is_Dipole)
937	<	skdat.t2 -= cross(uz_j, duduz_j);
938	<	}
939	<	}
940	<
941	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1007	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1008		RealType mu1, preVal, chg1, self;
1009
1010		if (!initialized_) initialize();
1011	<
1012	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1011	>
1012	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1013
1014		// logicals
949	–
1015		bool i_is_Charge = data.is_Charge;
1016		bool i_is_Dipole = data.is_Dipole;
1017
1018	<	if (summationMethod_ == REACTION_FIELD) {
1018	>	if (summationMethod_ == esm_REACTION_FIELD) {
1019		if (i_is_Dipole) {
1020		mu1 = data.dipole_moment;
1021		preVal = pre22_ * preRF2_ * mu1 * mu1;
1022	<	scdat.pot -= 0.5 * preVal;
1022	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1023
1024		// The self-correction term adds into the reaction field vector
1025	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1025	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1026		Vector3d ei = preVal * uz_i;
1027
1028		// This looks very wrong. A vector crossed with itself is zero.
1029	<	scdat.t -= cross(uz_i, ei);
1029	>	*(sdat.t) -= cross(uz_i, ei);
1030		}
1031	<	} else if (summationMethod_ == SHIFTED_FORCE \|\| summationMethod_ == SHIFTED_POTENTIAL) {
1031	>	} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1032		if (i_is_Charge) {
1033		chg1 = data.charge;
1034		if (screeningMethod_ == DAMPED) {
1035	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1035	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1036		} else {
1037	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1037	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1038		}
1039	<	scdat.pot += self;
1039	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1040		}
1041		}
1042		}
1043	+
1044	+	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1045	+	// This seems to work moderately well as a default. There's no
1046	+	// inherent scale for 1/r interactions that we can standardize.
1047	+	// 12 angstroms seems to be a reasonably good guess for most
1048	+	// cases.
1049	+	return 12.0;
1050	+	}
1051		}

Diff Legend

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents): Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs. Revision 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC

Diff Legend

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1504 by gezelter, Sat Oct 2 20:41:53 2010 UTC vs.
Revision 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC