[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 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC

#	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
57		void Electrostatic::initialize() {
58	+
59	+	Globals* simParams_ = info_->getSimParams();
60	+
61	+	summationMap_["HARD"] = esm_HARD;
62	+	summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
63	+	summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
64	+	summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
65	+	summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD;
66	+	summationMap_["EWALD_FULL"] = esm_EWALD_FULL;
67	+	summationMap_["EWALD_PME"] = esm_EWALD_PME;
68	+	summationMap_["EWALD_SPME"] = esm_EWALD_SPME;
69	+	screeningMap_["DAMPED"] = DAMPED;
70	+	screeningMap_["UNDAMPED"] = UNDAMPED;
71	+
72		// these prefactors convert the multipole interactions into kcal / mol
73		// all were computed assuming distances are measured in angstroms
74		// Charge-Charge, assuming charges are measured in electrons
#	Line 79 \| Line 93 \| namespace OpenMD {
93
94		// variables to handle different summation methods for long-range
95		// electrostatics:
96	<	summationMethod_ = NONE;
96	>	summationMethod_ = esm_HARD;
97		screeningMethod_ = UNDAMPED;
98		dielectric_ = 1.0;
99		one_third_ = 1.0 / 3.0;
100	<	haveDefaultCutoff_ = false;
100	>	haveCutoffRadius_ = false;
101		haveDampingAlpha_ = false;
102		haveDielectric_ = false;
103		haveElectroSpline_ = false;
104
105	+	// check the summation method:
106	+	if (simParams_->haveElectrostaticSummationMethod()) {
107	+	string myMethod = simParams_->getElectrostaticSummationMethod();
108	+	toUpper(myMethod);
109	+	map<string, ElectrostaticSummationMethod>::iterator i;
110	+	i = summationMap_.find(myMethod);
111	+	if ( i != summationMap_.end() ) {
112	+	summationMethod_ = (*i).second;
113	+	} else {
114	+	// throw error
115	+	sprintf( painCave.errMsg,
116	+	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
117	+	"\t(Input file specified %s .)\n"
118	+	"\telectrostaticSummationMethod must be one of: \"none\",\n"
119	+	"\t\"shifted_potential\", \"shifted_force\", or \n"
120	+	"\t\"reaction_field\".\n", myMethod.c_str() );
121	+	painCave.isFatal = 1;
122	+	simError();
123	+	}
124	+	} else {
125	+	// set ElectrostaticSummationMethod to the cutoffMethod:
126	+	if (simParams_->haveCutoffMethod()){
127	+	string myMethod = simParams_->getCutoffMethod();
128	+	toUpper(myMethod);
129	+	map<string, ElectrostaticSummationMethod>::iterator i;
130	+	i = summationMap_.find(myMethod);
131	+	if ( i != summationMap_.end() ) {
132	+	summationMethod_ = (*i).second;
133	+	}
134	+	}
135	+	}
136	+
137	+	if (summationMethod_ == esm_REACTION_FIELD) {
138	+	if (!simParams_->haveDielectric()) {
139	+	// throw warning
140	+	sprintf( painCave.errMsg,
141	+	"SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
142	+	"\tA default value of %f will be used for the dielectric.\n", dielectric_);
143	+	painCave.isFatal = 0;
144	+	painCave.severity = OPENMD_INFO;
145	+	simError();
146	+	} else {
147	+	dielectric_ = simParams_->getDielectric();
148	+	}
149	+	haveDielectric_ = true;
150	+	}
151	+
152	+	if (simParams_->haveElectrostaticScreeningMethod()) {
153	+	string myScreen = simParams_->getElectrostaticScreeningMethod();
154	+	toUpper(myScreen);
155	+	map<string, ElectrostaticScreeningMethod>::iterator i;
156	+	i = screeningMap_.find(myScreen);
157	+	if ( i != screeningMap_.end()) {
158	+	screeningMethod_ = (*i).second;
159	+	} else {
160	+	sprintf( painCave.errMsg,
161	+	"SimInfo error: Unknown electrostaticScreeningMethod.\n"
162	+	"\t(Input file specified %s .)\n"
163	+	"\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
164	+	"or \"damped\".\n", myScreen.c_str() );
165	+	painCave.isFatal = 1;
166	+	simError();
167	+	}
168	+	}
169	+
170	+	// check to make sure a cutoff value has been set:
171	+	if (!haveCutoffRadius_) {
172	+	sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
173	+	"Cutoff value!\n");
174	+	painCave.severity = OPENMD_ERROR;
175	+	painCave.isFatal = 1;
176	+	simError();
177	+	}
178	+
179	+	if (screeningMethod_ == DAMPED) {
180	+	if (!simParams_->haveDampingAlpha()) {
181	+	// first set a cutoff dependent alpha value
182	+	// we assume alpha depends linearly with rcut from 0 to 20.5 ang
183	+	dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
184	+	if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
185	+
186	+	// throw warning
187	+	sprintf( painCave.errMsg,
188	+	"Electrostatic::initialize: dampingAlpha was not specified in the input file.\n"
189	+	"\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n",
190	+	dampingAlpha_, cutoffRadius_);
191	+	painCave.severity = OPENMD_INFO;
192	+	painCave.isFatal = 0;
193	+	simError();
194	+	} else {
195	+	dampingAlpha_ = simParams_->getDampingAlpha();
196	+	}
197	+	haveDampingAlpha_ = true;
198	+	}
199	+
200		// find all of the Electrostatic atom Types:
201		ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
202		ForceField::AtomTypeContainer::MapTypeIterator i;
203		AtomType* at;
204	<
204	>
205		for (at = atomTypes->beginType(i); at != NULL;
206		at = atomTypes->nextType(i)) {
207
208		if (at->isElectrostatic())
209		addType(at);
101	–	}
102	–
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();
210		}
211	+
212
213	<	defaultCutoff2_ = defaultCutoff_ * defaultCutoff_;
214	<	rcuti_ = 1.0 / defaultCutoff_;
213	>	cutoffRadius2_ = cutoffRadius_ * cutoffRadius_;
214	>	rcuti_ = 1.0 / cutoffRadius_;
215		rcuti2_ = rcuti_ * rcuti_;
216		rcuti3_ = rcuti2_ * rcuti_;
217		rcuti4_ = rcuti2_ * rcuti2_;
218
219		if (screeningMethod_ == DAMPED) {
220	<	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	<
220	>
221		alpha2_ = dampingAlpha_ * dampingAlpha_;
222		alpha4_ = alpha2_ * alpha2_;
223		alpha6_ = alpha4_ * alpha2_;
224		alpha8_ = alpha4_ * alpha4_;
225
226	<	constEXP_ = exp(-alpha2_ * defaultCutoff2_);
226	>	constEXP_ = exp(-alpha2_ * cutoffRadius2_);
227		invRootPi_ = 0.56418958354775628695;
228		alphaPi_ = 2.0 * dampingAlpha_ * invRootPi_;
229
230	<	c1c_ = erfc(dampingAlpha_ * defaultCutoff_) * rcuti_;
230	>	c1c_ = erfc(dampingAlpha_ * cutoffRadius_) * rcuti_;
231		c2c_ = alphaPi_ * constEXP_ * rcuti_ + c1c_ * rcuti_;
232		c3c_ = 2.0 * alphaPi_ * alpha2_ + 3.0 * c2c_ * rcuti_;
233		c4c_ = 4.0 * alphaPi_ * alpha4_ + 5.0 * c3c_ * rcuti2_;
#	Line 148 \| Line 242 \| namespace OpenMD {
242		c6c_ = 9.0 * c5c_ * rcuti2_;
243		}
244
245	<	if (summationMethod_ == REACTION_FIELD) {
246	<	if (haveDielectric_) {
247	<	preRF_ = (dielectric_ - 1.0) /
248	<	((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	<	}
245	>	if (summationMethod_ == esm_REACTION_FIELD) {
246	>	preRF_ = (dielectric_ - 1.0) /
247	>	((2.0 * dielectric_ + 1.0) * cutoffRadius2_ * cutoffRadius_);
248	>	preRF2_ = 2.0 * preRF_;
249		}
250	<
251	<	RealType dx = defaultCutoff_ / RealType(np_ - 1);
250	>
251	>	RealType dx = cutoffRadius_ / RealType(np_ - 1);
252		RealType rval;
253		vector<RealType> rvals;
254		vector<RealType> yvals;
#	Line 283 \| Line 369 \| namespace OpenMD {
369		simError();
370		}
371
372	+	// Quadrupoles in OpenMD are set as the diagonal elements
373	+	// of the diagonalized traceless quadrupole moment tensor.
374	+	// The column vectors of the unitary matrix that diagonalizes
375	+	// the quadrupole moment tensor become the eFrame (or the
376	+	// electrostatic version of the body-fixed frame.
377	+
378		Vector3dGenericData* v3dData = dynamic_cast<Vector3dGenericData*>(data);
379		if (v3dData == NULL) {
380		sprintf( painCave.errMsg,
#	Line 315 \| Line 407 \| namespace OpenMD {
407		return;
408		}
409
410	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
411	<	RealType theRSW ) {
412	<	defaultCutoff_ = theECR;
413	<	rrf_ = defaultCutoff_;
322	<	rt_ = theRSW;
323	<	haveDefaultCutoff_ = true;
410	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
411	>	cutoffRadius_ = rCut;
412	>	rrf_ = cutoffRadius_;
413	>	haveCutoffRadius_ = true;
414		}
415	+
416	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
417	+	rt_ = rSwitch;
418	+	}
419		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
420		summationMethod_ = esm;
421		}
#	Line 337 \| Line 431 \| namespace OpenMD {
431		haveDielectric_ = true;
432		}
433
434	<	void Electrostatic::calcForce(InteractionData idat) {
434	>	void Electrostatic::calcForce(InteractionData &idat) {
435
436		// utility variables. Should clean these up and use the Vector3d and
437		// Mat3x3d to replace as many as we can in future versions:
#	Line 371 \| Line 465 \| namespace OpenMD {
465
466		if (!initialized_) initialize();
467
468	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
469	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
468	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
469	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
470
471		// some variables we'll need independent of electrostatic type:
472
473	<	riji = 1.0 / idat.rij;
474	<	Vector3d rhat = idat.d * riji;
473	>	riji = 1.0 / *(idat.rij) ;
474	>	Vector3d rhat = (idat.d) riji;
475
476		// logicals
477
#	Line 396 \| Line 490 \| namespace OpenMD {
490
491		if (i_is_Dipole) {
492		mu_i = data1.dipole_moment;
493	<	uz_i = idat.eFrame1.getColumn(2);
493	>	uz_i = idat.eFrame1->getColumn(2);
494
495		ct_i = dot(uz_i, rhat);
496
#	Line 412 \| Line 506 \| namespace OpenMD {
506		qyy_i = Q_i.y();
507		qzz_i = Q_i.z();
508
509	<	ux_i = idat.eFrame1.getColumn(0);
510	<	uy_i = idat.eFrame1.getColumn(1);
511	<	uz_i = idat.eFrame1.getColumn(2);
509	>	ux_i = idat.eFrame1->getColumn(0);
510	>	uy_i = idat.eFrame1->getColumn(1);
511	>	uz_i = idat.eFrame1->getColumn(2);
512
513		cx_i = dot(ux_i, rhat);
514		cy_i = dot(uy_i, rhat);
#	Line 430 \| Line 524 \| namespace OpenMD {
524
525		if (j_is_Dipole) {
526		mu_j = data2.dipole_moment;
527	<	uz_j = idat.eFrame2.getColumn(2);
527	>	uz_j = idat.eFrame2->getColumn(2);
528
529		ct_j = dot(uz_j, rhat);
530
#	Line 446 \| Line 540 \| namespace OpenMD {
540		qyy_j = Q_j.y();
541		qzz_j = Q_j.z();
542
543	<	ux_j = idat.eFrame2.getColumn(0);
544	<	uy_j = idat.eFrame2.getColumn(1);
545	<	uz_j = idat.eFrame2.getColumn(2);
543	>	ux_j = idat.eFrame2->getColumn(0);
544	>	uy_j = idat.eFrame2->getColumn(1);
545	>	uz_j = idat.eFrame2->getColumn(2);
546
547		cx_j = dot(ux_j, rhat);
548		cy_j = dot(uy_j, rhat);
#	Line 467 \| Line 561 \| namespace OpenMD {
561		if (j_is_Charge) {
562		if (screeningMethod_ == DAMPED) {
563		// assemble the damping variables
564	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
564	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
565		erfcVal = res.first;
566		derfcVal = res.second;
567		c1 = erfcVal * riji;
#	Line 477 \| Line 571 \| namespace OpenMD {
571		c2 = c1 * riji;
572		}
573
574	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
574	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
575
576	<	if (summationMethod_ == SHIFTED_POTENTIAL) {
576	>	if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
577		vterm = preVal * (c1 - c1c_);
578	<	dudr = -idat.sw * preVal * c2;
578	>	dudr = - (idat.sw) preVal * c2;
579
580	<	} else if (summationMethod_ == SHIFTED_FORCE) {
581	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - defaultCutoff_) );
582	<	dudr = idat.sw * preVal * (c2c_ - c2);
580	>	} else if (summationMethod_ == esm_SHIFTED_FORCE) {
581	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
582	>	dudr = (idat.sw) preVal * (c2c_ - c2);
583
584	<	} else if (summationMethod_ == REACTION_FIELD) {
585	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
584	>	} else if (summationMethod_ == esm_REACTION_FIELD) {
585	>	rfVal = (idat.electroMult) preRF_ * (idat.rij) *(idat.rij) ;
586		vterm = preVal * ( riji + rfVal );
587	<	dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
587	>	dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
588
589		} else {
590		vterm = preVal * riji * erfcVal;
591
592	<	dudr = - idat.sw * preVal * c2;
592	>	dudr = - (idat.sw) preVal * c2;
593
594		}
595
596	<	idat.vpair += vterm;
597	<	epot += idat.sw * vterm;
596	>	*(idat.vpair) += vterm;
597	>	epot += (idat.sw) vterm;
598
599		dVdr += dudr * rhat;
600		}
601
602		if (j_is_Dipole) {
603		// pref is used by all the possible methods
604	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
605	<	preSw = idat.sw * pref;
604	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
605	>	preSw = (idat.sw) pref;
606
607	<	if (summationMethod_ == REACTION_FIELD) {
607	>	if (summationMethod_ == esm_REACTION_FIELD) {
608		ri2 = riji * riji;
609		ri3 = ri2 * riji;
610
611	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
612	<	idat.vpair += vterm;
613	<	epot += idat.sw * vterm;
611	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
612	>	*(idat.vpair) += vterm;
613	>	epot += (idat.sw) vterm;
614
615		dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
616	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
616	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
617
618		} else {
619		// determine the inverse r used if we have split dipoles
620		if (j_is_SplitDipole) {
621	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
621	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
622		ri = 1.0 / BigR;
623	<	scale = idat.rij * ri;
623	>	scale = (idat.rij) ri;
624		} else {
625		ri = riji;
626		scale = 1.0;
#	Line 536 \| Line 630 \| namespace OpenMD {
630
631		if (screeningMethod_ == DAMPED) {
632		// assemble the damping variables
633	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
633	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
634		erfcVal = res.first;
635		derfcVal = res.second;
636		c1 = erfcVal * ri;
#	Line 553 \| Line 647 \| namespace OpenMD {
647		// calculate the potential
648		pot_term = scale * c2;
649		vterm = -pref * ct_j * pot_term;
650	<	idat.vpair += vterm;
651	<	epot += idat.sw * vterm;
650	>	*(idat.vpair) += vterm;
651	>	epot += (idat.sw) vterm;
652
653		// calculate derivatives for forces and torques
654
#	Line 569 \| Line 663 \| namespace OpenMD {
663		cx2 = cx_j * cx_j;
664		cy2 = cy_j * cy_j;
665		cz2 = cz_j * cz_j;
666	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
666	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
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 * riji;
#	Line 588 \| Line 682 \| namespace OpenMD {
682		}
683
684		// precompute variables for convenience
685	<	preSw = idat.sw * pref;
685	>	preSw = (idat.sw) pref;
686		c2ri = c2 * riji;
687		c3ri = c3 * riji;
688	<	c4rij = c4 * idat.rij;
688	>	c4rij = c4 * *(idat.rij) ;
689		rhatdot2 = 2.0 * rhat * c3;
690		rhatc4 = rhat * c4rij;
691
#	Line 600 \| Line 694 \| namespace OpenMD {
694		qyy_j * (cy2*c3 - c2ri) +
695		qzz_j * (cz2*c3 - c2ri) );
696		vterm = pref * pot_term;
697	<	idat.vpair += vterm;
698	<	epot += idat.sw * vterm;
697	>	*(idat.vpair) += vterm;
698	>	epot += (idat.sw) vterm;
699
700		// calculate derivatives for the forces and torques
701
#	Line 619 \| Line 713 \| namespace OpenMD {
713
714		if (j_is_Charge) {
715		// variables used by all the methods
716	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
717	<	preSw = idat.sw * pref;
716	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
717	>	preSw = (idat.sw) pref;
718
719	<	if (summationMethod_ == REACTION_FIELD) {
719	>	if (summationMethod_ == esm_REACTION_FIELD) {
720
721		ri2 = riji * riji;
722		ri3 = ri2 * riji;
723
724	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
725	<	idat.vpair += vterm;
726	<	epot += idat.sw * vterm;
724	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
725	>	*(idat.vpair) += vterm;
726	>	epot += (idat.sw) vterm;
727
728		dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
729
730	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
730	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
731
732		} else {
733
734		// determine inverse r if we are using split dipoles
735		if (i_is_SplitDipole) {
736	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
736	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
737		ri = 1.0 / BigR;
738	<	scale = idat.rij * ri;
738	>	scale = (idat.rij) ri;
739		} else {
740		ri = riji;
741		scale = 1.0;
#	Line 651 \| Line 745 \| namespace OpenMD {
745
746		if (screeningMethod_ == DAMPED) {
747		// assemble the damping variables
748	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
748	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
749		erfcVal = res.first;
750		derfcVal = res.second;
751		c1 = erfcVal * ri;
#	Line 668 \| Line 762 \| namespace OpenMD {
762		// calculate the potential
763		pot_term = c2 * scale;
764		vterm = pref * ct_i * pot_term;
765	<	idat.vpair += vterm;
766	<	epot += idat.sw * vterm;
765	>	*(idat.vpair) += vterm;
766	>	epot += (idat.sw) vterm;
767
768		// calculate derivatives for the forces and torques
769		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
#	Line 681 \| Line 775 \| namespace OpenMD {
775		// variables used by all methods
776		ct_ij = dot(uz_i, uz_j);
777
778	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
779	<	preSw = idat.sw * pref;
778	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
779	>	preSw = (idat.sw) pref;
780
781	<	if (summationMethod_ == REACTION_FIELD) {
781	>	if (summationMethod_ == esm_REACTION_FIELD) {
782		ri2 = riji * riji;
783		ri3 = ri2 * riji;
784		ri4 = ri2 * ri2;
785
786		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
787		preRF2_ * ct_ij );
788	<	idat.vpair += vterm;
789	<	epot += idat.sw * vterm;
788	>	*(idat.vpair) += vterm;
789	>	epot += (idat.sw) vterm;
790
791		a1 = 5.0 * ct_i * ct_j - ct_ij;
792
#	Line 705 \| Line 799 \| namespace OpenMD {
799
800		if (i_is_SplitDipole) {
801		if (j_is_SplitDipole) {
802	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
802	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
803		} else {
804	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
804	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
805		}
806		ri = 1.0 / BigR;
807	<	scale = idat.rij * ri;
807	>	scale = (idat.rij) ri;
808		} else {
809		if (j_is_SplitDipole) {
810	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
810	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
811		ri = 1.0 / BigR;
812	<	scale = idat.rij * ri;
812	>	scale = (idat.rij) ri;
813		} else {
814		ri = riji;
815		scale = 1.0;
#	Line 723 \| Line 817 \| namespace OpenMD {
817		}
818		if (screeningMethod_ == DAMPED) {
819		// assemble damping variables
820	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
820	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
821		erfcVal = res.first;
822		derfcVal = res.second;
823		c1 = erfcVal * ri;
#	Line 745 \| Line 839 \| namespace OpenMD {
839		preSwSc = preSw * scale;
840		c2ri = c2 * ri;
841		c3ri = c3 * ri;
842	<	c4rij = c4 * idat.rij;
842	>	c4rij = c4 * *(idat.rij) ;
843
844		// calculate the potential
845		pot_term = (ct_ij * c2ri - ctidotj * c3);
846		vterm = pref * pot_term;
847	<	idat.vpair += vterm;
848	<	epot += idat.sw * vterm;
847	>	*(idat.vpair) += vterm;
848	>	epot += (idat.sw) vterm;
849
850		// calculate derivatives for the forces and torques
851		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 770 \| Line 864 \| namespace OpenMD {
864		cy2 = cy_i * cy_i;
865		cz2 = cz_i * cz_i;
866
867	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
867	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
868
869		if (screeningMethod_ == DAMPED) {
870		// assemble the damping variables
871	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
871	>	res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
872		erfcVal = res.first;
873		derfcVal = res.second;
874		c1 = erfcVal * riji;
#	Line 789 \| Line 883 \| namespace OpenMD {
883		}
884
885		// precompute some variables for convenience
886	<	preSw = idat.sw * pref;
886	>	preSw = (idat.sw) pref;
887		c2ri = c2 * riji;
888		c3ri = c3 * riji;
889	<	c4rij = c4 * idat.rij;
889	>	c4rij = c4 * *(idat.rij) ;
890		rhatdot2 = 2.0 * rhat * c3;
891		rhatc4 = rhat * c4rij;
892
#	Line 802 \| Line 896 \| namespace OpenMD {
896		qzz_i * (cz2 * c3 - c2ri) );
897
898		vterm = pref * pot_term;
899	<	idat.vpair += vterm;
900	<	epot += idat.sw * vterm;
899	>	*(idat.vpair) += vterm;
900	>	epot += (idat.sw) vterm;
901
902		// calculate the derivatives for the forces and torques
903
#	Line 817 \| Line 911 \| namespace OpenMD {
911		}
912		}
913
914	<	idat.pot += epot;
915	<	idat.f1 += dVdr;
914	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
915	>	*(idat.f1) += dVdr;
916
917		if (i_is_Dipole \|\| i_is_Quadrupole)
918	<	idat.t1 -= cross(uz_i, duduz_i);
918	>	*(idat.t1) -= cross(uz_i, duduz_i);
919		if (i_is_Quadrupole) {
920	<	idat.t1 -= cross(ux_i, dudux_i);
921	<	idat.t1 -= cross(uy_i, duduy_i);
920	>	*(idat.t1) -= cross(ux_i, dudux_i);
921	>	*(idat.t1) -= cross(uy_i, duduy_i);
922		}
923	<
924	<	if (j_is_Dipole \|\| j_is_Quadrupole)
925	<	idat.t2 -= cross(uz_j, duduz_j);
923	>
924	>	if (j_is_Dipole \|\| j_is_Quadrupole)
925	>	*(idat.t2) -= cross(uz_j, duduz_j);
926		if (j_is_Quadrupole) {
927	<	idat.t2 -= cross(uz_j, dudux_j);
928	<	idat.t2 -= cross(uz_j, duduy_j);
927	>	*(idat.t2) -= cross(uz_j, dudux_j);
928	>	*(idat.t2) -= cross(uz_j, duduy_j);
929		}
930
931		return;
932		}
933
934	<	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
934	>	void Electrostatic::calcSkipCorrection(InteractionData &idat) {
935
936		if (!initialized_) initialize();
937
938	<	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
939	<	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
938	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
939	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
940
941		// logicals
942
#	Line 859 \| Line 953 \| namespace OpenMD {
953
954		if (i_is_Charge) {
955		q_i = data1.charge;
956	<	skdat.skippedCharge2 += q_i;
956	>	*(idat.skippedCharge2) += q_i;
957		}
958
959		if (j_is_Charge) {
960		q_j = data2.charge;
961	<	skdat.skippedCharge1 += q_j;
961	>	*(idat.skippedCharge1) += q_j;
962		}
963
964		// the rest of this function should only be necessary for reaction field.
965
966	<	if (summationMethod_ == REACTION_FIELD) {
966	>	if (summationMethod_ == esm_REACTION_FIELD) {
967		RealType riji, ri2, ri3;
968	<	RealType q_i, mu_i, ct_i;
969	<	RealType q_j, mu_j, ct_j;
970	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
968	>	RealType mu_i, ct_i;
969	>	RealType mu_j, ct_j;
970	>	RealType preVal, rfVal, vterm, dudr, pref, myPot(0.0);
971		Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
972
973		// some variables we'll need independent of electrostatic type:
974
975	<	riji = 1.0 / skdat.rij;
976	<	rhat = skdat.d * riji;
975	>	riji = 1.0 / *(idat.rij) ;
976	>	rhat = (idat.d) riji;
977
978		if (i_is_Dipole) {
979		mu_i = data1.dipole_moment;
980	<	uz_i = skdat.eFrame1.getColumn(2);
980	>	uz_i = idat.eFrame1->getColumn(2);
981		ct_i = dot(uz_i, rhat);
982		duduz_i = V3Zero;
983		}
984
985		if (j_is_Dipole) {
986		mu_j = data2.dipole_moment;
987	<	uz_j = skdat.eFrame2.getColumn(2);
987	>	uz_j = idat.eFrame2->getColumn(2);
988		ct_j = dot(uz_j, rhat);
989		duduz_j = V3Zero;
990		}
991
992		if (i_is_Charge) {
993		if (j_is_Charge) {
994	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
995	<	rfVal = preRF_ * skdat.rij * skdat.rij;
994	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
995	>	rfVal = preRF_ * (idat.rij) *(idat.rij) ;
996		vterm = preVal * rfVal;
997	<	myPot += skdat.sw * vterm;
998	<	dudr = skdat.sw * preVal * 2.0 * rfVal * riji;
997	>	myPot += (idat.sw) vterm;
998	>	dudr = (idat.sw) preVal * 2.0 * rfVal * riji;
999		dVdr += dudr * rhat;
1000		}
1001
1002		if (j_is_Dipole) {
1003		ri2 = riji * riji;
1004		ri3 = ri2 * riji;
1005	<	pref = skdat.electroMult * pre12_ * q_i * mu_j;
1006	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1007	<	myPot += skdat.sw * vterm;
1008	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1005	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
1006	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
1007	>	myPot += (idat.sw) vterm;
1008	>	dVdr += - (idat.sw) pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1009	>	duduz_j += - (idat.sw) pref * rhat * (ri2 - preRF2_ * *(idat.rij) );
1010		}
1011		}
1012		if (i_is_Dipole) {
1013		if (j_is_Charge) {
1014		ri2 = riji * riji;
1015		ri3 = ri2 * riji;
1016	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
1017	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1018	<	myPot += skdat.sw * vterm;
1019	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1020	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1016	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
1017	>	vterm = - pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
1018	>	myPot += (idat.sw) vterm;
1019	>	dVdr += (idat.sw) pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1020	>	duduz_i += (idat.sw) pref * rhat * (ri2 - preRF2_ * *(idat.rij));
1021		}
1022		}
1023
1024		// accumulate the forces and torques resulting from the self term
1025	<	skdat.pot += myPot;
1026	<	skdat.f1 += dVdr;
1025	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += myPot;
1026	>	*(idat.f1) += dVdr;
1027
1028		if (i_is_Dipole)
1029	<	skdat.t1 -= cross(uz_i, duduz_i);
1029	>	*(idat.t1) -= cross(uz_i, duduz_i);
1030		if (j_is_Dipole)
1031	<	skdat.t2 -= cross(uz_j, duduz_j);
1031	>	*(idat.t2) -= cross(uz_j, duduz_j);
1032		}
1033		}
1034
1035	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1035	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1036		RealType mu1, preVal, chg1, self;
1037
1038		if (!initialized_) initialize();
1039
1040	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1040	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1041
1042		// logicals
1043
1044		bool i_is_Charge = data.is_Charge;
1045		bool i_is_Dipole = data.is_Dipole;
1046
1047	<	if (summationMethod_ == REACTION_FIELD) {
1047	>	if (summationMethod_ == esm_REACTION_FIELD) {
1048		if (i_is_Dipole) {
1049		mu1 = data.dipole_moment;
1050		preVal = pre22_ * preRF2_ * mu1 * mu1;
1051	<	scdat.pot -= 0.5 * preVal;
1051	>	sdat.pot[2] -= 0.5 * preVal;
1052
1053		// The self-correction term adds into the reaction field vector
1054	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1054	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1055		Vector3d ei = preVal * uz_i;
1056
1057		// This looks very wrong. A vector crossed with itself is zero.
1058	<	scdat.t -= cross(uz_i, ei);
1058	>	*(sdat.t) -= cross(uz_i, ei);
1059		}
1060	<	} else if (summationMethod_ == SHIFTED_FORCE \|\| summationMethod_ == SHIFTED_POTENTIAL) {
1060	>	} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1061		if (i_is_Charge) {
1062		chg1 = data.charge;
1063		if (screeningMethod_ == DAMPED) {
1064	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1064	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1065		} else {
1066	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1066	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1067		}
1068	<	scdat.pot += self;
1068	>	sdat.pot[ELECTROSTATIC_FAMILY] += self;
1069		}
1070		}
1071		}
1072	+
1073	+	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1074	+	// This seems to work moderately well as a default. There's no
1075	+	// inherent scale for 1/r interactions that we can standardize.
1076	+	// 12 angstroms seems to be a reasonably good guess for most
1077	+	// cases.
1078	+	return 12.0;
1079	+	}
1080		}

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 1584 by gezelter, Fri Jun 17 20:16: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 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC