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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC vs.
Revision 1665 by gezelter, Tue Nov 22 20:38:56 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).
39	>	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	>	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		#include <stdio.h>
#	Line 47 \| Line 48
48		#include "utils/simError.h"
49		#include "types/NonBondedInteractionType.hpp"
50		#include "types/DirectionalAtomType.hpp"
51	+	#include "io/Globals.hpp"
52
51	–
53		namespace OpenMD {
54
55		Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
56	<	forceField_(NULL), info_(NULL) {}
56	>	forceField_(NULL), info_(NULL),
57	>	haveCutoffRadius_(false),
58	>	haveDampingAlpha_(false),
59	>	haveDielectric_(false),
60	>	haveElectroSpline_(false)
61	>	{}
62
63		void Electrostatic::initialize() {
64	+
65	+	Globals* simParams_ = info_->getSimParams();
66
67		summationMap_["HARD"] = esm_HARD;
68	+	summationMap_["NONE"] = esm_HARD;
69		summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
70		summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
71		summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
#	Line 95 \| Line 104 \| namespace OpenMD {
104		screeningMethod_ = UNDAMPED;
105		dielectric_ = 1.0;
106		one_third_ = 1.0 / 3.0;
98	–	haveCutoffRadius_ = false;
99	–	haveDampingAlpha_ = false;
100	–	haveDielectric_ = false;
101	–	haveElectroSpline_ = false;
107
108		// check the summation method:
109		if (simParams_->haveElectrostaticSummationMethod()) {
#	Line 111 \| Line 116 \| namespace OpenMD {
116		} else {
117		// throw error
118		sprintf( painCave.errMsg,
119	<	"SimInfo error: Unknown electrostaticSummationMethod.\n"
119	>	"Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
120		"\t(Input file specified %s .)\n"
121	<	"\telectrostaticSummationMethod must be one of: \"none\",\n"
121	>	"\telectrostaticSummationMethod must be one of: \"hard\",\n"
122		"\t\"shifted_potential\", \"shifted_force\", or \n"
123		"\t\"reaction_field\".\n", myMethod.c_str() );
124		painCave.isFatal = 1;
#	Line 246 \| Line 251 \| namespace OpenMD {
251		preRF2_ = 2.0 * preRF_;
252		}
253
254	<	RealType dx = cutoffRadius_ / RealType(np_ - 1);
254	>	// Add a 2 angstrom safety window to deal with cutoffGroups that
255	>	// have charged atoms longer than the cutoffRadius away from each
256	>	// other. Splining may not be the best choice here. Direct calls
257	>	// to erfc might be preferrable.
258	>
259	>	RealType dx = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
260		RealType rval;
261		vector<RealType> rvals;
262		vector<RealType> yvals;
#	Line 405 \| Line 415 \| namespace OpenMD {
415		return;
416		}
417
418	<	void Electrostatic::setElectrostaticCutoffRadius( RealType theECR,
419	<	RealType theRSW ) {
420	<	cutoffRadius_ = theECR;
411	<	rrf_ = cutoffRadius_;
412	<	rt_ = theRSW;
418	>	void Electrostatic::setCutoffRadius( RealType rCut ) {
419	>	cutoffRadius_ = rCut;
420	>	rrf_ = cutoffRadius_;
421		haveCutoffRadius_ = true;
422		}
423	+
424	+	void Electrostatic::setSwitchingRadius( RealType rSwitch ) {
425	+	rt_ = rSwitch;
426	+	}
427		void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
428		summationMethod_ = esm;
429		}
#	Line 427 \| Line 439 \| namespace OpenMD {
439		haveDielectric_ = true;
440		}
441
442	<	void Electrostatic::calcForce(InteractionData idat) {
442	>	void Electrostatic::calcForce(InteractionData &idat) {
443
444		// utility variables. Should clean these up and use the Vector3d and
445		// Mat3x3d to replace as many as we can in future versions:
#	Line 441 \| Line 453 \| namespace OpenMD {
453		RealType ct_i, ct_j, ct_ij, a1;
454		RealType riji, ri, ri2, ri3, ri4;
455		RealType pref, vterm, epot, dudr;
456	+	RealType vpair(0.0);
457		RealType scale, sc2;
458		RealType pot_term, preVal, rfVal;
459		RealType c2ri, c3ri, c4rij, cti3, ctj3, ctidotj;
460		RealType preSw, preSwSc;
461		RealType c1, c2, c3, c4;
462	<	RealType erfcVal, derfcVal;
462	>	RealType erfcVal(1.0), derfcVal(0.0);
463		RealType BigR;
464
465		Vector3d Q_i, Q_j;
#	Line 457 \| Line 470 \| namespace OpenMD {
470		Vector3d rhatdot2, rhatc4;
471		Vector3d dVdr;
472
473	+	// variables for indirect (reaction field) interactions for excluded pairs:
474	+	RealType indirect_Pot(0.0);
475	+	RealType indirect_vpair(0.0);
476	+	Vector3d indirect_dVdr(V3Zero);
477	+	Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero);
478	+
479		pair<RealType, RealType> res;
480
481		if (!initialized_) initialize();
482
483	<	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atype1];
484	<	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atype2];
483	>	ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first];
484	>	ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second];
485
486		// some variables we'll need independent of electrostatic type:
487
488	<	riji = 1.0 / idat.rij;
489	<	Vector3d rhat = idat.d * riji;
488	>	riji = 1.0 / *(idat.rij) ;
489	>	Vector3d rhat = (idat.d) riji;
490
491		// logicals
492
#	Line 481 \| Line 500 \| namespace OpenMD {
500		bool j_is_SplitDipole = data2.is_SplitDipole;
501		bool j_is_Quadrupole = data2.is_Quadrupole;
502
503	<	if (i_is_Charge)
503	>	if (i_is_Charge) {
504		q_i = data1.charge;
505	+	if (idat.excluded) {
506	+	*(idat.skippedCharge2) += q_i;
507	+	}
508	+	}
509
510		if (i_is_Dipole) {
511		mu_i = data1.dipole_moment;
512	<	uz_i = idat.eFrame1.getColumn(2);
512	>	uz_i = idat.eFrame1->getColumn(2);
513
514		ct_i = dot(uz_i, rhat);
515
#	Line 502 \| Line 525 \| namespace OpenMD {
525		qyy_i = Q_i.y();
526		qzz_i = Q_i.z();
527
528	<	ux_i = idat.eFrame1.getColumn(0);
529	<	uy_i = idat.eFrame1.getColumn(1);
530	<	uz_i = idat.eFrame1.getColumn(2);
528	>	ux_i = idat.eFrame1->getColumn(0);
529	>	uy_i = idat.eFrame1->getColumn(1);
530	>	uz_i = idat.eFrame1->getColumn(2);
531
532		cx_i = dot(ux_i, rhat);
533		cy_i = dot(uy_i, rhat);
#	Line 515 \| Line 538 \| namespace OpenMD {
538		duduz_i = V3Zero;
539		}
540
541	<	if (j_is_Charge)
541	>	if (j_is_Charge) {
542		q_j = data2.charge;
543	+	if (idat.excluded) {
544	+	*(idat.skippedCharge1) += q_j;
545	+	}
546	+	}
547
548	+
549		if (j_is_Dipole) {
550		mu_j = data2.dipole_moment;
551	<	uz_j = idat.eFrame2.getColumn(2);
551	>	uz_j = idat.eFrame2->getColumn(2);
552
553		ct_j = dot(uz_j, rhat);
554
#	Line 536 \| Line 564 \| namespace OpenMD {
564		qyy_j = Q_j.y();
565		qzz_j = Q_j.z();
566
567	<	ux_j = idat.eFrame2.getColumn(0);
568	<	uy_j = idat.eFrame2.getColumn(1);
569	<	uz_j = idat.eFrame2.getColumn(2);
567	>	ux_j = idat.eFrame2->getColumn(0);
568	>	uy_j = idat.eFrame2->getColumn(1);
569	>	uz_j = idat.eFrame2->getColumn(2);
570
571		cx_j = dot(ux_j, rhat);
572		cy_j = dot(uy_j, rhat);
#	Line 557 \| Line 585 \| namespace OpenMD {
585		if (j_is_Charge) {
586		if (screeningMethod_ == DAMPED) {
587		// assemble the damping variables
588	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
589	<	erfcVal = res.first;
590	<	derfcVal = res.second;
588	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
589	>	//erfcVal = res.first;
590	>	//derfcVal = res.second;
591	>
592	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
593	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
594	>
595		c1 = erfcVal * riji;
596		c2 = (-derfcVal + c1) * riji;
597		} else {
#	Line 567 \| Line 599 \| namespace OpenMD {
599		c2 = c1 * riji;
600		}
601
602	<	preVal = idat.electroMult * pre11_ * q_i * q_j;
602	>	preVal = (idat.electroMult) pre11_ * q_i * q_j;
603
604		if (summationMethod_ == esm_SHIFTED_POTENTIAL) {
605		vterm = preVal * (c1 - c1c_);
606	<	dudr = -idat.sw * preVal * c2;
606	>	dudr = - (idat.sw) preVal * c2;
607
608		} else if (summationMethod_ == esm_SHIFTED_FORCE) {
609	<	vterm = preVal * ( c1 - c1c_ + c2c_*(idat.rij - cutoffRadius_) );
610	<	dudr = idat.sw * preVal * (c2c_ - c2);
609	>	vterm = preVal * ( c1 - c1c_ + c2c_( (idat.rij) - cutoffRadius_) );
610	>	dudr = (idat.sw) preVal * (c2c_ - c2);
611
612		} else if (summationMethod_ == esm_REACTION_FIELD) {
613	<	rfVal = idat.electroMult * preRF_ * idat.rij * idat.rij;
613	>	rfVal = preRF_ * (idat.rij) *(idat.rij);
614	>
615		vterm = preVal * ( riji + rfVal );
616	<	dudr = idat.sw * preVal * ( 2.0 * rfVal - riji ) * riji;
616	>	dudr = (idat.sw) preVal * ( 2.0 * rfVal - riji ) * riji;
617	>
618	>	// if this is an excluded pair, there are still indirect
619	>	// interactions via the reaction field we must worry about:
620
621	+	if (idat.excluded) {
622	+	indirect_vpair += preVal * rfVal;
623	+	indirect_Pot += (idat.sw) preVal * rfVal;
624	+	indirect_dVdr += (idat.sw) preVal * 2.0 * rfVal * riji * rhat;
625	+	}
626	+
627		} else {
586	–	vterm = preVal * riji * erfcVal;
628
629	<	dudr = - idat.sw * preVal * c2;
629	>	vterm = preVal * riji * erfcVal;
630	>	dudr = - (idat.sw) preVal * c2;
631
632		}
591	–
592	–	idat.vpair += vterm;
593	–	epot += idat.sw * vterm;
633
634	<	dVdr += dudr * rhat;
634	>	vpair += vterm;
635	>	epot += (idat.sw) vterm;
636	>	dVdr += dudr * rhat;
637		}
638
639		if (j_is_Dipole) {
640		// pref is used by all the possible methods
641	<	pref = idat.electroMult * pre12_ * q_i * mu_j;
642	<	preSw = idat.sw * pref;
641	>	pref = (idat.electroMult) pre12_ * q_i * mu_j;
642	>	preSw = (idat.sw) pref;
643
644		if (summationMethod_ == esm_REACTION_FIELD) {
645		ri2 = riji * riji;
646		ri3 = ri2 * riji;
647
648	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * idat.rij );
649	<	idat.vpair += vterm;
650	<	epot += idat.sw * vterm;
648	>	vterm = - pref * ct_j * ( ri2 - preRF2_ * *(idat.rij) );
649	>	vpair += vterm;
650	>	epot += (idat.sw) vterm;
651
652		dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
653	<	duduz_j += -preSw * rhat * (ri2 - preRF2_ * idat.rij);
653	>	duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
654
655	+	// Even if we excluded this pair from direct interactions,
656	+	// we still have the reaction-field-mediated charge-dipole
657	+	// interaction:
658	+
659	+	if (idat.excluded) {
660	+	indirect_vpair += pref * ct_j * preRF2_ * *(idat.rij);
661	+	indirect_Pot += preSw * ct_j * preRF2_ * *(idat.rij);
662	+	indirect_dVdr += preSw * preRF2_ * uz_j;
663	+	indirect_duduz_j += preSw * rhat * preRF2_ * *(idat.rij);
664	+	}
665	+
666		} else {
667		// determine the inverse r used if we have split dipoles
668		if (j_is_SplitDipole) {
669	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
669	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
670		ri = 1.0 / BigR;
671	<	scale = idat.rij * ri;
671	>	scale = (idat.rij) ri;
672		} else {
673		ri = riji;
674		scale = 1.0;
#	Line 626 \| Line 678 \| namespace OpenMD {
678
679		if (screeningMethod_ == DAMPED) {
680		// assemble the damping variables
681	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
682	<	erfcVal = res.first;
683	<	derfcVal = res.second;
681	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
682	>	//erfcVal = res.first;
683	>	//derfcVal = res.second;
684	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
685	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
686		c1 = erfcVal * ri;
687		c2 = (-derfcVal + c1) * ri;
688		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 643 \| Line 697 \| namespace OpenMD {
697		// calculate the potential
698		pot_term = scale * c2;
699		vterm = -pref * ct_j * pot_term;
700	<	idat.vpair += vterm;
701	<	epot += idat.sw * vterm;
700	>	vpair += vterm;
701	>	epot += (idat.sw) vterm;
702
703		// calculate derivatives for forces and torques
704
#	Line 659 \| Line 713 \| namespace OpenMD {
713		cx2 = cx_j * cx_j;
714		cy2 = cy_j * cy_j;
715		cz2 = cz_j * cz_j;
716	<	pref = idat.electroMult * pre14_ * q_i * one_third_;
716	>	pref = (idat.electroMult) pre14_ * q_i * one_third_;
717
718		if (screeningMethod_ == DAMPED) {
719		// assemble the damping variables
720	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
721	<	erfcVal = res.first;
722	<	derfcVal = res.second;
720	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
721	>	//erfcVal = res.first;
722	>	//derfcVal = res.second;
723	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
724	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
725		c1 = erfcVal * riji;
726		c2 = (-derfcVal + c1) * riji;
727		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 678 \| Line 734 \| namespace OpenMD {
734		}
735
736		// precompute variables for convenience
737	<	preSw = idat.sw * pref;
737	>	preSw = (idat.sw) pref;
738		c2ri = c2 * riji;
739		c3ri = c3 * riji;
740	<	c4rij = c4 * idat.rij;
740	>	c4rij = c4 * *(idat.rij) ;
741		rhatdot2 = 2.0 * rhat * c3;
742		rhatc4 = rhat * c4rij;
743
#	Line 690 \| Line 746 \| namespace OpenMD {
746		qyy_j * (cy2*c3 - c2ri) +
747		qzz_j * (cz2*c3 - c2ri) );
748		vterm = pref * pot_term;
749	<	idat.vpair += vterm;
750	<	epot += idat.sw * vterm;
749	>	vpair += vterm;
750	>	epot += (idat.sw) vterm;
751
752		// calculate derivatives for the forces and torques
753
#	Line 709 \| Line 765 \| namespace OpenMD {
765
766		if (j_is_Charge) {
767		// variables used by all the methods
768	<	pref = idat.electroMult * pre12_ * q_j * mu_i;
769	<	preSw = idat.sw * pref;
768	>	pref = (idat.electroMult) pre12_ * q_j * mu_i;
769	>	preSw = (idat.sw) pref;
770
771		if (summationMethod_ == esm_REACTION_FIELD) {
772
773		ri2 = riji * riji;
774		ri3 = ri2 * riji;
775
776	<	vterm = pref * ct_i * ( ri2 - preRF2_ * idat.rij );
777	<	idat.vpair += vterm;
778	<	epot += idat.sw * vterm;
776	>	vterm = pref * ct_i * ( ri2 - preRF2_ * *(idat.rij) );
777	>	vpair += vterm;
778	>	epot += (idat.sw) vterm;
779
780		dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
781
782	<	duduz_i += preSw * rhat * (ri2 - preRF2_ * idat.rij);
782	>	duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) );
783	>
784	>	// Even if we excluded this pair from direct interactions,
785	>	// we still have the reaction-field-mediated charge-dipole
786	>	// interaction:
787	>
788	>	if (idat.excluded) {
789	>	indirect_vpair += -pref * ct_i * preRF2_ * *(idat.rij);
790	>	indirect_Pot += -preSw * ct_i * preRF2_ * *(idat.rij);
791	>	indirect_dVdr += -preSw * preRF2_ * uz_i;
792	>	indirect_duduz_i += -preSw * rhat * preRF2_ * *(idat.rij);
793	>	}
794
795		} else {
796
797		// determine inverse r if we are using split dipoles
798		if (i_is_SplitDipole) {
799	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
799	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
800		ri = 1.0 / BigR;
801	<	scale = idat.rij * ri;
801	>	scale = (idat.rij) ri;
802		} else {
803		ri = riji;
804		scale = 1.0;
#	Line 741 \| Line 808 \| namespace OpenMD {
808
809		if (screeningMethod_ == DAMPED) {
810		// assemble the damping variables
811	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
812	<	erfcVal = res.first;
813	<	derfcVal = res.second;
811	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
812	>	//erfcVal = res.first;
813	>	//derfcVal = res.second;
814	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
815	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
816		c1 = erfcVal * ri;
817		c2 = (-derfcVal + c1) * ri;
818		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 758 \| Line 827 \| namespace OpenMD {
827		// calculate the potential
828		pot_term = c2 * scale;
829		vterm = pref * ct_i * pot_term;
830	<	idat.vpair += vterm;
831	<	epot += idat.sw * vterm;
830	>	vpair += vterm;
831	>	epot += (idat.sw) vterm;
832
833		// calculate derivatives for the forces and torques
834		dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3);
#	Line 771 \| Line 840 \| namespace OpenMD {
840		// variables used by all methods
841		ct_ij = dot(uz_i, uz_j);
842
843	<	pref = idat.electroMult * pre22_ * mu_i * mu_j;
844	<	preSw = idat.sw * pref;
843	>	pref = (idat.electroMult) pre22_ * mu_i * mu_j;
844	>	preSw = (idat.sw) pref;
845
846		if (summationMethod_ == esm_REACTION_FIELD) {
847		ri2 = riji * riji;
#	Line 781 \| Line 850 \| namespace OpenMD {
850
851		vterm = pref * ( ri3 * (ct_ij - 3.0 * ct_i * ct_j) -
852		preRF2_ * ct_ij );
853	<	idat.vpair += vterm;
854	<	epot += idat.sw * vterm;
853	>	vpair += vterm;
854	>	epot += (idat.sw) vterm;
855
856		a1 = 5.0 * ct_i * ct_j - ct_ij;
857
#	Line 791 \| Line 860 \| namespace OpenMD {
860		duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j);
861		duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i);
862
863	+	if (idat.excluded) {
864	+	indirect_vpair += - pref * preRF2_ * ct_ij;
865	+	indirect_Pot += - preSw * preRF2_ * ct_ij;
866	+	indirect_duduz_i += -preSw * preRF2_ * uz_j;
867	+	indirect_duduz_j += -preSw * preRF2_ * uz_i;
868	+	}
869	+
870		} else {
871
872		if (i_is_SplitDipole) {
873		if (j_is_SplitDipole) {
874	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i + 0.25 * d_j * d_j);
874	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i + 0.25 * d_j * d_j);
875		} else {
876	<	BigR = sqrt(idat.r2 + 0.25 * d_i * d_i);
876	>	BigR = sqrt( (idat.r2) + 0.25 d_i * d_i);
877		}
878		ri = 1.0 / BigR;
879	<	scale = idat.rij * ri;
879	>	scale = (idat.rij) ri;
880		} else {
881		if (j_is_SplitDipole) {
882	<	BigR = sqrt(idat.r2 + 0.25 * d_j * d_j);
882	>	BigR = sqrt( (idat.r2) + 0.25 d_j * d_j);
883		ri = 1.0 / BigR;
884	<	scale = idat.rij * ri;
884	>	scale = (idat.rij) ri;
885		} else {
886		ri = riji;
887		scale = 1.0;
#	Line 813 \| Line 889 \| namespace OpenMD {
889		}
890		if (screeningMethod_ == DAMPED) {
891		// assemble damping variables
892	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
893	<	erfcVal = res.first;
894	<	derfcVal = res.second;
892	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
893	>	//erfcVal = res.first;
894	>	//derfcVal = res.second;
895	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
896	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
897		c1 = erfcVal * ri;
898		c2 = (-derfcVal + c1) * ri;
899		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * ri;
#	Line 835 \| Line 913 \| namespace OpenMD {
913		preSwSc = preSw * scale;
914		c2ri = c2 * ri;
915		c3ri = c3 * ri;
916	<	c4rij = c4 * idat.rij;
916	>	c4rij = c4 * *(idat.rij) ;
917
918		// calculate the potential
919		pot_term = (ct_ij * c2ri - ctidotj * c3);
920		vterm = pref * pot_term;
921	<	idat.vpair += vterm;
922	<	epot += idat.sw * vterm;
921	>	vpair += vterm;
922	>	epot += (idat.sw) vterm;
923
924		// calculate derivatives for the forces and torques
925		dVdr += preSwSc * ( ctidotj * rhat * c4rij -
#	Line 860 \| Line 938 \| namespace OpenMD {
938		cy2 = cy_i * cy_i;
939		cz2 = cz_i * cz_i;
940
941	<	pref = idat.electroMult * pre14_ * q_j * one_third_;
941	>	pref = (idat.electroMult) pre14_ * q_j * one_third_;
942
943		if (screeningMethod_ == DAMPED) {
944		// assemble the damping variables
945	<	res = erfcSpline_->getValueAndDerivativeAt(idat.rij);
946	<	erfcVal = res.first;
947	<	derfcVal = res.second;
945	>	//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) );
946	>	//erfcVal = res.first;
947	>	//derfcVal = res.second;
948	>	erfcVal = erfc(dampingAlpha_ * *(idat.rij));
949	>	derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2));
950		c1 = erfcVal * riji;
951		c2 = (-derfcVal + c1) * riji;
952		c3 = -2.0 * derfcVal * alpha2_ + 3.0 * c2 * riji;
#	Line 879 \| Line 959 \| namespace OpenMD {
959		}
960
961		// precompute some variables for convenience
962	<	preSw = idat.sw * pref;
962	>	preSw = (idat.sw) pref;
963		c2ri = c2 * riji;
964		c3ri = c3 * riji;
965	<	c4rij = c4 * idat.rij;
965	>	c4rij = c4 * *(idat.rij) ;
966		rhatdot2 = 2.0 * rhat * c3;
967		rhatc4 = rhat * c4rij;
968
#	Line 892 \| Line 972 \| namespace OpenMD {
972		qzz_i * (cz2 * c3 - c2ri) );
973
974		vterm = pref * pot_term;
975	<	idat.vpair += vterm;
976	<	epot += idat.sw * vterm;
975	>	vpair += vterm;
976	>	epot += (idat.sw) vterm;
977
978		// calculate the derivatives for the forces and torques
979
#	Line 907 \| Line 987 \| namespace OpenMD {
987		}
988		}
989
910	–	idat.pot += epot;
911	–	idat.f1 += dVdr;
990
991	<	if (i_is_Dipole \|\| i_is_Quadrupole)
992	<	idat.t1 -= cross(uz_i, duduz_i);
993	<	if (i_is_Quadrupole) {
994	<	idat.t1 -= cross(ux_i, dudux_i);
995	<	idat.t1 -= cross(uy_i, duduy_i);
996	<	}
991	>	if (!idat.excluded) {
992	>	*(idat.vpair) += vpair;
993	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += epot;
994	>	*(idat.f1) += dVdr;
995	>
996	>	if (i_is_Dipole \|\| i_is_Quadrupole)
997	>	*(idat.t1) -= cross(uz_i, duduz_i);
998	>	if (i_is_Quadrupole) {
999	>	*(idat.t1) -= cross(ux_i, dudux_i);
1000	>	*(idat.t1) -= cross(uy_i, duduy_i);
1001	>	}
1002	>
1003	>	if (j_is_Dipole \|\| j_is_Quadrupole)
1004	>	*(idat.t2) -= cross(uz_j, duduz_j);
1005	>	if (j_is_Quadrupole) {
1006	>	*(idat.t2) -= cross(uz_j, dudux_j);
1007	>	*(idat.t2) -= cross(uz_j, duduy_j);
1008	>	}
1009
1010	<	if (j_is_Dipole \|\| j_is_Quadrupole)
1011	<	idat.t2 -= cross(uz_j, duduz_j);
1012	<	if (j_is_Quadrupole) {
1013	<	idat.t2 -= cross(uz_j, dudux_j);
1014	<	idat.t2 -= cross(uz_j, duduy_j);
1010	>	} else {
1011	>
1012	>	// only accumulate the forces and torques resulting from the
1013	>	// indirect reaction field terms.
1014	>
1015	>	*(idat.vpair) += indirect_vpair;
1016	>	(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot;
1017	>	*(idat.f1) += indirect_dVdr;
1018	>
1019	>	if (i_is_Dipole)
1020	>	*(idat.t1) -= cross(uz_i, indirect_duduz_i);
1021	>	if (j_is_Dipole)
1022	>	*(idat.t2) -= cross(uz_j, indirect_duduz_j);
1023		}
1024
1025	+
1026		return;
1027		}
929	–
930	–	void Electrostatic::calcSkipCorrection(SkipCorrectionData skdat) {
931	–
932	–	if (!initialized_) initialize();
1028
1029	<	ElectrostaticAtomData data1 = ElectrostaticMap[skdat.atype1];
935	<	ElectrostaticAtomData data2 = ElectrostaticMap[skdat.atype2];
936	<
937	<	// logicals
938	<
939	<	bool i_is_Charge = data1.is_Charge;
940	<	bool i_is_Dipole = data1.is_Dipole;
941	<
942	<	bool j_is_Charge = data2.is_Charge;
943	<	bool j_is_Dipole = data2.is_Dipole;
944	<
945	<	RealType q_i, q_j;
946	<
947	<	// The skippedCharge computation is needed by the real-space cutoff methods
948	<	// (i.e. shifted force and shifted potential)
949	<
950	<	if (i_is_Charge) {
951	<	q_i = data1.charge;
952	<	skdat.skippedCharge2 += q_i;
953	<	}
954	<
955	<	if (j_is_Charge) {
956	<	q_j = data2.charge;
957	<	skdat.skippedCharge1 += q_j;
958	<	}
959	<
960	<	// the rest of this function should only be necessary for reaction field.
961	<
962	<	if (summationMethod_ == esm_REACTION_FIELD) {
963	<	RealType riji, ri2, ri3;
964	<	RealType q_i, mu_i, ct_i;
965	<	RealType q_j, mu_j, ct_j;
966	<	RealType preVal, rfVal, vterm, dudr, pref, myPot;
967	<	Vector3d dVdr, uz_i, uz_j, duduz_i, duduz_j, rhat;
968	<
969	<	// some variables we'll need independent of electrostatic type:
970	<
971	<	riji = 1.0 / skdat.rij;
972	<	rhat = skdat.d * riji;
973	<
974	<	if (i_is_Dipole) {
975	<	mu_i = data1.dipole_moment;
976	<	uz_i = skdat.eFrame1.getColumn(2);
977	<	ct_i = dot(uz_i, rhat);
978	<	duduz_i = V3Zero;
979	<	}
980	<
981	<	if (j_is_Dipole) {
982	<	mu_j = data2.dipole_moment;
983	<	uz_j = skdat.eFrame2.getColumn(2);
984	<	ct_j = dot(uz_j, rhat);
985	<	duduz_j = V3Zero;
986	<	}
987	<
988	<	if (i_is_Charge) {
989	<	if (j_is_Charge) {
990	<	preVal = skdat.electroMult * pre11_ * q_i * q_j;
991	<	rfVal = preRF_ * skdat.rij * skdat.rij;
992	<	vterm = preVal * rfVal;
993	<	myPot += skdat.sw * vterm;
994	<	dudr = skdat.sw * preVal * 2.0 * rfVal * riji;
995	<	dVdr += dudr * rhat;
996	<	}
997	<
998	<	if (j_is_Dipole) {
999	<	ri2 = riji * riji;
1000	<	ri3 = ri2 * riji;
1001	<	pref = skdat.electroMult * pre12_ * q_i * mu_j;
1002	<	vterm = - pref * ct_j * ( ri2 - preRF2_ * skdat.rij );
1003	<	myPot += skdat.sw * vterm;
1004	<	dVdr += -skdat.sw * pref * ( ri3 * ( uz_j - 3.0 * ct_j * rhat) - preRF2_ * uz_j);
1005	<	duduz_j += -skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1006	<	}
1007	<	}
1008	<	if (i_is_Dipole) {
1009	<	if (j_is_Charge) {
1010	<	ri2 = riji * riji;
1011	<	ri3 = ri2 * riji;
1012	<	pref = skdat.electroMult * pre12_ * q_j * mu_i;
1013	<	vterm = - pref * ct_i * ( ri2 - preRF2_ * skdat.rij );
1014	<	myPot += skdat.sw * vterm;
1015	<	dVdr += skdat.sw * pref * ( ri3 * ( uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i);
1016	<	duduz_i += skdat.sw * pref * rhat * (ri2 - preRF2_ * skdat.rij);
1017	<	}
1018	<	}
1019	<
1020	<	// accumulate the forces and torques resulting from the self term
1021	<	skdat.pot += myPot;
1022	<	skdat.f1 += dVdr;
1023	<
1024	<	if (i_is_Dipole)
1025	<	skdat.t1 -= cross(uz_i, duduz_i);
1026	<	if (j_is_Dipole)
1027	<	skdat.t2 -= cross(uz_j, duduz_j);
1028	<	}
1029	<	}
1030	<
1031	<	void Electrostatic::calcSelfCorrection(SelfCorrectionData scdat) {
1029	>	void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1030		RealType mu1, preVal, chg1, self;
1031
1032		if (!initialized_) initialize();
1033	<
1034	<	ElectrostaticAtomData data = ElectrostaticMap[scdat.atype];
1033	>
1034	>	ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1035
1036		// logicals
1039	–
1037		bool i_is_Charge = data.is_Charge;
1038		bool i_is_Dipole = data.is_Dipole;
1039
#	Line 1044 \| Line 1041 \| namespace OpenMD {
1041		if (i_is_Dipole) {
1042		mu1 = data.dipole_moment;
1043		preVal = pre22_ * preRF2_ * mu1 * mu1;
1044	<	scdat.pot -= 0.5 * preVal;
1044	>	((sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 preVal;
1045
1046		// The self-correction term adds into the reaction field vector
1047	<	Vector3d uz_i = scdat.eFrame.getColumn(2);
1047	>	Vector3d uz_i = sdat.eFrame->getColumn(2);
1048		Vector3d ei = preVal * uz_i;
1049
1050		// This looks very wrong. A vector crossed with itself is zero.
1051	<	scdat.t -= cross(uz_i, ei);
1051	>	*(sdat.t) -= cross(uz_i, ei);
1052		}
1053		} else if (summationMethod_ == esm_SHIFTED_FORCE \|\| summationMethod_ == esm_SHIFTED_POTENTIAL) {
1054		if (i_is_Charge) {
1055		chg1 = data.charge;
1056		if (screeningMethod_ == DAMPED) {
1057	<	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1057	>	self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1058		} else {
1059	<	self = - 0.5 * rcuti_ * chg1 * (chg1 + scdat.skippedCharge) * pre11_;
1059	>	self = - 0.5 * rcuti_ * chg1 * (chg1 + (sdat.skippedCharge)) pre11_;
1060		}
1061	<	scdat.pot += self;
1061	>	(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1062		}
1063		}
1064		}
1065
1066	<	RealType Electrostatic::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
1066	>	RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType, AtomType> atypes) {
1067		// This seems to work moderately well as a default. There's no
1068		// inherent scale for 1/r interactions that we can standardize.
1069		// 12 angstroms seems to be a reasonably good guess for most

Diff Legend

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

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents): Revision 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC vs. Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

Diff Legend

Comparing branches/development/src/nonbonded/Electrostatic.cpp (file contents):
Revision 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC vs.
Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC