ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2118 by gezelter, Fri Mar 11 15:53:18 2005 UTC vs.
Revision 2390 by chrisfen, Wed Oct 19 19:24:40 2005 UTC

#	Line 40 | Line 40 | module electrostatic_module
40		!!
41
42		module electrostatic_module
43	<
43	>
44		use force_globals
45		use definitions
46		use atype_module
#	Line 54 | Line 54 | module electrostatic_module
54
55		PRIVATE
56
57	+
58	+	#define __FORTRAN90
59	+	#include "UseTheForce/DarkSide/fInteractionMap.h"
60	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61	+
62	+
63		!! these prefactors convert the multipole interactions into kcal / mol
64		!! all were computed assuming distances are measured in angstroms
65		!! Charge-Charge, assuming charges are measured in electrons
#	Line 68 | Line 74 | module electrostatic_module
74		!! This unit is also known affectionately as an esu centi-barn.
75		real(kind=dp), parameter :: pre14 = 69.13373_dp
76
77	+	!! variables to handle different summation methods for long-range electrostatics:
78	+	integer, save :: summationMethod = NONE
79	+	logical, save :: summationMethodChecked = .false.
80	+	real(kind=DP), save :: defaultCutoff = 0.0_DP
81	+	real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82	+	logical, save :: haveDefaultCutoff = .false.
83	+	real(kind=DP), save :: dampingAlpha = 0.0_DP
84	+	logical, save :: haveDampingAlpha = .false.
85	+	real(kind=DP), save :: dielectric = 1.0_DP
86	+	logical, save :: haveDielectric = .false.
87	+	real(kind=DP), save :: constERFC = 0.0_DP
88	+	real(kind=DP), save :: constEXP = 0.0_DP
89	+	logical, save :: haveDWAconstants = .false.
90	+	real(kind=dp), save :: rcuti = 0.0_DP
91	+	real(kind=dp), save :: rcuti2 = 0.0_DP
92	+	real(kind=dp), save :: rcuti3 = 0.0_DP
93	+	real(kind=dp), save :: rcuti4 = 0.0_DP
94	+	real(kind=dp), save :: alphaPi = 0.0_DP
95	+	real(kind=dp), save :: invRootPi = 0.0_DP
96	+	real(kind=dp), save :: rrf = 1.0_DP
97	+	real(kind=dp), save :: rt = 1.0_DP
98	+	real(kind=dp), save :: rrfsq = 1.0_DP
99	+	real(kind=dp), save :: preRF = 0.0_DP
100	+	logical, save :: preRFCalculated = .false.
101	+
102	+	#ifdef __IFC
103	+	! error function for ifc version > 7.
104	+	double precision, external :: derfc
105	+	#endif
106	+
107	+	public :: setElectrostaticSummationMethod
108	+	public :: setElectrostaticCutoffRadius
109	+	public :: setDampedWolfAlpha
110	+	public :: setReactionFieldDielectric
111	+	public :: setReactionFieldPrefactor
112		public :: newElectrostaticType
113		public :: setCharge
114		public :: setDipoleMoment
#	Line 76 | Line 117 | module electrostatic_module
117		public :: doElectrostaticPair
118		public :: getCharge
119		public :: getDipoleMoment
120	+	public :: pre22
121	+	public :: destroyElectrostaticTypes
122	+	public :: accumulate_rf
123	+	public :: accumulate_self_rf
124	+	public :: reaction_field_final
125	+	public :: rf_correct_forces
126
127		type :: Electrostatic
128		integer :: c_ident
#	Line 83 | Line 130 | module electrostatic_module
130		logical :: is_Dipole = .false.
131		logical :: is_SplitDipole = .false.
132		logical :: is_Quadrupole = .false.
133	+	logical :: is_Tap = .false.
134		real(kind=DP) :: charge = 0.0_DP
135		real(kind=DP) :: dipole_moment = 0.0_DP
136		real(kind=DP) :: split_dipole_distance = 0.0_DP
#	Line 93 | Line 141 | contains
141
142		contains
143
144	+	subroutine setElectrostaticSummationMethod(the_ESM)
145	+	integer, intent(in) :: the_ESM
146	+
147	+	if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
148	+	call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
149	+	endif
150	+
151	+	summationMethod = the_ESM
152	+
153	+	end subroutine setElectrostaticSummationMethod
154	+
155	+	subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
156	+	real(kind=dp), intent(in) :: thisRcut
157	+	real(kind=dp), intent(in) :: thisRsw
158	+	defaultCutoff = thisRcut
159	+	rrf = defaultCutoff
160	+	rt = thisRsw
161	+	haveDefaultCutoff = .true.
162	+	end subroutine setElectrostaticCutoffRadius
163	+
164	+	subroutine setDampedWolfAlpha(thisAlpha)
165	+	real(kind=dp), intent(in) :: thisAlpha
166	+	dampingAlpha = thisAlpha
167	+	haveDampingAlpha = .true.
168	+	end subroutine setDampedWolfAlpha
169	+
170	+	subroutine setReactionFieldDielectric(thisDielectric)
171	+	real(kind=dp), intent(in) :: thisDielectric
172	+	dielectric = thisDielectric
173	+	haveDielectric = .true.
174	+	end subroutine setReactionFieldDielectric
175	+
176	+	subroutine setReactionFieldPrefactor
177	+	if (haveDefaultCutoff .and. haveDielectric) then
178	+	defaultCutoff2 = defaultCutoff*defaultCutoff
179	+	preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
180	+	((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
181	+	preRFCalculated = .true.
182	+	else if (.not.haveDefaultCutoff) then
183	+	call handleError("setReactionFieldPrefactor", "Default cutoff not set")
184	+	else
185	+	call handleError("setReactionFieldPrefactor", "Dielectric not set")
186	+	endif
187	+	end subroutine setReactionFieldPrefactor
188	+
189		subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
190	<	is_SplitDipole, is_Quadrupole, status)
191	<
190	>	is_SplitDipole, is_Quadrupole, is_Tap, status)
191	>
192		integer, intent(in) :: c_ident
193		logical, intent(in) :: is_Charge
194		logical, intent(in) :: is_Dipole
195		logical, intent(in) :: is_SplitDipole
196		logical, intent(in) :: is_Quadrupole
197	+	logical, intent(in) :: is_Tap
198		integer, intent(out) :: status
199		integer :: nAtypes, myATID, i, j
200
201		status = 0
202		myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
203	<
203	>
204		!! Be simple-minded and assume that we need an ElectrostaticMap that
205		!! is the same size as the total number of atom types
206
207		if (.not.allocated(ElectrostaticMap)) then
208	<
208	>
209		nAtypes = getSize(atypes)
210	<
210	>
211		if (nAtypes == 0) then
212		status = -1
213		return
214		end if
215	<
215	>
216		if (.not. allocated(ElectrostaticMap)) then
217		allocate(ElectrostaticMap(nAtypes))
218		endif
219	<
219	>
220		end if
221
222		if (myATID .gt. size(ElectrostaticMap)) then
223		status = -1
224		return
225		endif
226	<
226	>
227		! set the values for ElectrostaticMap for this atom type:
228
229		ElectrostaticMap(myATID)%c_ident = c_ident
#	Line 137 | Line 231 | contains
231		ElectrostaticMap(myATID)%is_Dipole = is_Dipole
232		ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
233		ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
234	<
234	>	ElectrostaticMap(myATID)%is_Tap = is_Tap
235	>
236		end subroutine newElectrostaticType
237
238		subroutine setCharge(c_ident, charge, status)
#	Line 165 | Line 260 | contains
260		call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
261		status = -1
262		return
263	<	endif
263	>	endif
264
265		ElectrostaticMap(myATID)%charge = charge
266		end subroutine setCharge
#	Line 256 | Line 351 | contains
351		status = -1
352		return
353		endif
354	<
354	>
355		do i = 1, 3
356	<	ElectrostaticMap(myATID)%quadrupole_moments(i) = &
357	<	quadrupole_moments(i)
358	<	enddo
356	>	ElectrostaticMap(myATID)%quadrupole_moments(i) = &
357	>	quadrupole_moments(i)
358	>	enddo
359
360		end subroutine setQuadrupoleMoments
361
362	<
362	>
363		function getCharge(atid) result (c)
364		integer, intent(in) :: atid
365		integer :: localError
366		real(kind=dp) :: c
367	<
367	>
368		if (.not.allocated(ElectrostaticMap)) then
369		call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
370		return
371		end if
372	<
372	>
373		if (.not.ElectrostaticMap(atid)%is_Charge) then
374		call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
375		return
376		endif
377	<
377	>
378		c = ElectrostaticMap(atid)%charge
379		end function getCharge
380
#	Line 287 | Line 382 | contains
382		integer, intent(in) :: atid
383		integer :: localError
384		real(kind=dp) :: dm
385	<
385	>
386		if (.not.allocated(ElectrostaticMap)) then
387		call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
388		return
389		end if
390	<
390	>
391		if (.not.ElectrostaticMap(atid)%is_Dipole) then
392		call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
393		return
394		endif
395	<
395	>
396		dm = ElectrostaticMap(atid)%dipole_moment
397		end function getDipoleMoment
398	+
399	+	subroutine checkSummationMethod()
400	+
401	+	if (.not.haveDefaultCutoff) then
402	+	call handleError("checkSummationMethod", "no Default Cutoff set!")
403	+	endif
404	+
405	+	rcuti = 1.0d0 / defaultCutoff
406	+	rcuti2 = rcuti*rcuti
407	+	rcuti3 = rcuti2*rcuti
408	+	rcuti4 = rcuti2*rcuti2
409	+
410	+	if (summationMethod .eq. DAMPED_WOLF) then
411	+	if (.not.haveDWAconstants) then
412	+
413	+	if (.not.haveDampingAlpha) then
414	+	call handleError("checkSummationMethod", "no Damping Alpha set!")
415	+	endif
416	+
417	+	if (.not.haveDefaultCutoff) then
418	+	call handleError("checkSummationMethod", "no Default Cutoff set!")
419	+	endif
420
421	+	constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
422	+	constERFC = derfc(dampingAlpha*defaultCutoff)
423	+	invRootPi = 0.56418958354775628695d0
424	+	alphaPi = 2*dampingAlpha*invRootPi
425	+
426	+	haveDWAconstants = .true.
427	+	endif
428	+	endif
429	+
430	+	if (summationMethod .eq. REACTION_FIELD) then
431	+	if (.not.haveDielectric) then
432	+	call handleError("checkSummationMethod", "no reaction field Dielectric set!")
433	+	endif
434	+	endif
435	+
436	+	summationMethodChecked = .true.
437	+	end subroutine checkSummationMethod
438	+
439	+
440	+
441		subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
442		vpair, fpair, pot, eFrame, f, t, do_pot)
443	<
443	>
444		logical, intent(in) :: do_pot
445	<
445	>
446		integer, intent(in) :: atom1, atom2
447		integer :: localError
448
#	Line 318 | Line 455 | contains
455		real( kind = dp ), dimension(9,nLocal) :: eFrame
456		real( kind = dp ), dimension(3,nLocal) :: f
457		real( kind = dp ), dimension(3,nLocal) :: t
321	–
322	–	real (kind = dp), dimension(3) :: ul_i
323	–	real (kind = dp), dimension(3) :: ul_j
458
459	+	real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
460	+	real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
461	+	real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
462	+	real (kind = dp), dimension(3) :: dudux_j, duduy_j, duduz_j
463	+
464		logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
465		logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
466	+	logical :: i_is_Tap, j_is_Tap
467		integer :: me1, me2, id1, id2
468		real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
469	+	real (kind=dp) :: qxx_i, qyy_i, qzz_i
470	+	real (kind=dp) :: qxx_j, qyy_j, qzz_j
471	+	real (kind=dp) :: cx_i, cy_i, cz_i
472	+	real (kind=dp) :: cx_j, cy_j, cz_j
473	+	real (kind=dp) :: cx2, cy2, cz2
474		real (kind=dp) :: ct_i, ct_j, ct_ij, a1
475		real (kind=dp) :: riji, ri, ri2, ri3, ri4
476	<	real (kind=dp) :: pref, vterm, epot, dudr
476	>	real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
477		real (kind=dp) :: xhat, yhat, zhat
478		real (kind=dp) :: dudx, dudy, dudz
334	–	real (kind=dp) :: drdxj, drdyj, drdzj
335	–	real (kind=dp) :: duduix, duduiy, duduiz, dudujx, dudujy, dudujz
479		real (kind=dp) :: scale, sc2, bigR
480	+	real (kind=dp) :: varERFC, varEXP
481	+	real (kind=dp) :: limScale
482
483		if (.not.allocated(ElectrostaticMap)) then
484		call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
485		return
486		end if
487
488	+	if (.not.summationMethodChecked) then
489	+	call checkSummationMethod()
490	+
491	+	endif
492	+
493	+
494		#ifdef IS_MPI
495		me1 = atid_Row(atom1)
496		me2 = atid_Col(atom2)
#	Line 351 | Line 502 | contains
502		!! some variables we'll need independent of electrostatic type:
503
504		riji = 1.0d0 / rij
505	<
505	>
506		xhat = d(1) * riji
507		yhat = d(2) * riji
508		zhat = d(3) * riji
509
359	–	drdxj = xhat
360	–	drdyj = yhat
361	–	drdzj = zhat
362	–
510		!! logicals
364	–
511		i_is_Charge = ElectrostaticMap(me1)%is_Charge
512		i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
513		i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
514		i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
515	+	i_is_Tap = ElectrostaticMap(me1)%is_Tap
516
517		j_is_Charge = ElectrostaticMap(me2)%is_Charge
518		j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
519		j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
520		j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
521	+	j_is_Tap = ElectrostaticMap(me2)%is_Tap
522
523		if (i_is_Charge) then
524		q_i = ElectrostaticMap(me1)%charge
525		endif
526	<
526	>
527		if (i_is_Dipole) then
528		mu_i = ElectrostaticMap(me1)%dipole_moment
529		#ifdef IS_MPI
530	<	ul_i(1) = eFrame_Row(3,atom1)
531	<	ul_i(2) = eFrame_Row(6,atom1)
532	<	ul_i(3) = eFrame_Row(9,atom1)
530	>	uz_i(1) = eFrame_Row(3,atom1)
531	>	uz_i(2) = eFrame_Row(6,atom1)
532	>	uz_i(3) = eFrame_Row(9,atom1)
533		#else
534	<	ul_i(1) = eFrame(3,atom1)
535	<	ul_i(2) = eFrame(6,atom1)
536	<	ul_i(3) = eFrame(9,atom1)
534	>	uz_i(1) = eFrame(3,atom1)
535	>	uz_i(2) = eFrame(6,atom1)
536	>	uz_i(3) = eFrame(9,atom1)
537		#endif
538	<	ct_i = ul_i(1)*drdxj + ul_i(2)*drdyj + ul_i(3)*drdzj
538	>	ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
539
540		if (i_is_SplitDipole) then
541		d_i = ElectrostaticMap(me1)%split_dipole_distance
542		endif
543	<
543	>
544		endif
545
546	+	if (i_is_Quadrupole) then
547	+	qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1)
548	+	qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2)
549	+	qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3)
550	+	#ifdef IS_MPI
551	+	ux_i(1) = eFrame_Row(1,atom1)
552	+	ux_i(2) = eFrame_Row(4,atom1)
553	+	ux_i(3) = eFrame_Row(7,atom1)
554	+	uy_i(1) = eFrame_Row(2,atom1)
555	+	uy_i(2) = eFrame_Row(5,atom1)
556	+	uy_i(3) = eFrame_Row(8,atom1)
557	+	uz_i(1) = eFrame_Row(3,atom1)
558	+	uz_i(2) = eFrame_Row(6,atom1)
559	+	uz_i(3) = eFrame_Row(9,atom1)
560	+	#else
561	+	ux_i(1) = eFrame(1,atom1)
562	+	ux_i(2) = eFrame(4,atom1)
563	+	ux_i(3) = eFrame(7,atom1)
564	+	uy_i(1) = eFrame(2,atom1)
565	+	uy_i(2) = eFrame(5,atom1)
566	+	uy_i(3) = eFrame(8,atom1)
567	+	uz_i(1) = eFrame(3,atom1)
568	+	uz_i(2) = eFrame(6,atom1)
569	+	uz_i(3) = eFrame(9,atom1)
570	+	#endif
571	+	cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
572	+	cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
573	+	cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
574	+	endif
575	+
576		if (j_is_Charge) then
577		q_j = ElectrostaticMap(me2)%charge
578		endif
579	<
579	>
580		if (j_is_Dipole) then
581		mu_j = ElectrostaticMap(me2)%dipole_moment
582		#ifdef IS_MPI
583	<	ul_j(1) = eFrame_Col(3,atom2)
584	<	ul_j(2) = eFrame_Col(6,atom2)
585	<	ul_j(3) = eFrame_Col(9,atom2)
583	>	uz_j(1) = eFrame_Col(3,atom2)
584	>	uz_j(2) = eFrame_Col(6,atom2)
585	>	uz_j(3) = eFrame_Col(9,atom2)
586		#else
587	<	ul_j(1) = eFrame(3,atom2)
588	<	ul_j(2) = eFrame(6,atom2)
589	<	ul_j(3) = eFrame(9,atom2)
587	>	uz_j(1) = eFrame(3,atom2)
588	>	uz_j(2) = eFrame(6,atom2)
589	>	uz_j(3) = eFrame(9,atom2)
590		#endif
591	<	ct_j = ul_j(1)*drdxj + ul_j(2)*drdyj + ul_j(3)*drdzj
591	>	ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
592
593		if (j_is_SplitDipole) then
594		d_j = ElectrostaticMap(me2)%split_dipole_distance
595		endif
596		endif
597
598	+	if (j_is_Quadrupole) then
599	+	qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1)
600	+	qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2)
601	+	qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3)
602	+	#ifdef IS_MPI
603	+	ux_j(1) = eFrame_Col(1,atom2)
604	+	ux_j(2) = eFrame_Col(4,atom2)
605	+	ux_j(3) = eFrame_Col(7,atom2)
606	+	uy_j(1) = eFrame_Col(2,atom2)
607	+	uy_j(2) = eFrame_Col(5,atom2)
608	+	uy_j(3) = eFrame_Col(8,atom2)
609	+	uz_j(1) = eFrame_Col(3,atom2)
610	+	uz_j(2) = eFrame_Col(6,atom2)
611	+	uz_j(3) = eFrame_Col(9,atom2)
612	+	#else
613	+	ux_j(1) = eFrame(1,atom2)
614	+	ux_j(2) = eFrame(4,atom2)
615	+	ux_j(3) = eFrame(7,atom2)
616	+	uy_j(1) = eFrame(2,atom2)
617	+	uy_j(2) = eFrame(5,atom2)
618	+	uy_j(3) = eFrame(8,atom2)
619	+	uz_j(1) = eFrame(3,atom2)
620	+	uz_j(2) = eFrame(6,atom2)
621	+	uz_j(3) = eFrame(9,atom2)
622	+	#endif
623	+	cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
624	+	cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
625	+	cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
626	+	endif
627	+
628		epot = 0.0_dp
629		dudx = 0.0_dp
630		dudy = 0.0_dp
631		dudz = 0.0_dp
632
633	<	duduix = 0.0_dp
634	<	duduiy = 0.0_dp
635	<	duduiz = 0.0_dp
633	>	dudux_i = 0.0_dp
634	>	duduy_i = 0.0_dp
635	>	duduz_i = 0.0_dp
636
637	<	dudujx = 0.0_dp
638	<	dudujy = 0.0_dp
639	<	dudujz = 0.0_dp
637	>	dudux_j = 0.0_dp
638	>	duduy_j = 0.0_dp
639	>	duduz_j = 0.0_dp
640
641		if (i_is_Charge) then
642
643		if (j_is_Charge) then
436	–
437	–	vterm = pre11 * q_i * q_j * riji
438	–	vpair = vpair + vterm
439	–	epot = epot + sw*vterm
644
645	<	dudr  = - sw * vterm * riji
645	>	if (summationMethod .eq. UNDAMPED_WOLF) then
646
647	<	dudx = dudx + dudr * drdxj
648	<	dudy = dudy + dudr * drdyj
649	<	dudz = dudz + dudr * drdzj
650	<
651	<	endif
647	>	vterm = pre11 * q_i * q_j * (riji - rcuti)
648	>	vpair = vpair + vterm
649	>	epot = epot + sw*vterm
650	>
651	>	dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
652	>
653	>	dudx = dudx + dudr * d(1)
654	>	dudy = dudy + dudr * d(2)
655	>	dudz = dudz + dudr * d(3)
656
657	<	if (j_is_Dipole) then
657	>	elseif (summationMethod .eq. DAMPED_WOLF) then
658
659	<	if (j_is_SplitDipole) then
660	<	BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
661	<	ri = 1.0_dp / BigR
662	<	scale = rij * ri
663	<	else
456	<	ri = riji
457	<	scale = 1.0_dp
458	<	endif
459	<
460	<	ri2 = ri * ri
461	<	ri3 = ri2 * ri
462	<	sc2 = scale * scale
659	>	varERFC = derfc(dampingAlpha*rij)
660	>	varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
661	>	vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
662	>	vpair = vpair + vterm
663	>	epot = epot + sw*vterm
664
665	<	pref = pre12 * q_i * mu_j
666	<	vterm = pref * ct_j * ri2 * scale
667	<	vpair = vpair + vterm
668	<	epot = epot + sw * vterm
665	>	dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
666	>	+ alphaPi*varEXP) &
667	>	- (constERFC*rcuti2 &
668	>	+ alphaPi*constEXP)) )
669	>
670	>	dudx = dudx + dudr * d(1)
671	>	dudy = dudy + dudr * d(2)
672	>	dudz = dudz + dudr * d(3)
673
469	–	!! this has a + sign in the () because the rij vector is
470	–	!! r_j - r_i and the charge-dipole potential takes the origin
471	–	!! as the point dipole, which is atom j in this case.
472	–
473	–	dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0*ct_j*xhat*sc2)
474	–	dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0*ct_j*yhat*sc2)
475	–	dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0*ct_j*zhat*sc2)
476	–
477	–	dudujx = dudujx - pref * sw * ri2 * xhat * scale
478	–	dudujy = dudujy - pref * sw * ri2 * yhat * scale
479	–	dudujz = dudujz - pref * sw * ri2 * zhat * scale
480	–
481	–	endif
482	–
483	–	endif
484	–
485	–	if (i_is_Dipole) then
486	–
487	–	if (j_is_Charge) then
488	–
489	–	if (i_is_SplitDipole) then
490	–	BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
491	–	ri = 1.0_dp / BigR
492	–	scale = rij * ri
674		else
494	–	ri = riji
495	–	scale = 1.0_dp
496	–	endif
675
676	<	ri2 = ri * ri
677	<	ri3 = ri2 * ri
678	<	sc2 = scale * scale
676	>	vterm = pre11 * q_i * q_j * riji
677	>	vpair = vpair + vterm
678	>	epot = epot + sw*vterm
679
680	<	pref = pre12 * q_j * mu_i
681	<	vterm = pref * ct_i * ri2 * scale
682	<	vpair = vpair + vterm
683	<	epot = epot + sw * vterm
680	>	dudr  = - sw * vterm * riji
681	>
682	>	dudx = dudx + dudr * xhat
683	>	dudy = dudy + dudr * yhat
684	>	dudz = dudz + dudr * zhat
685
686	<	dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * xhat*sc2)
508	<	dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * yhat*sc2)
509	<	dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * zhat*sc2)
686	>	endif
687
511	–	duduix = duduix + pref * sw * ri2 * xhat * scale
512	–	duduiy = duduiy + pref * sw * ri2 * yhat * scale
513	–	duduiz = duduiz + pref * sw * ri2 * zhat * scale
688		endif
689
690		if (j_is_Dipole) then
691
692	<	if (i_is_SplitDipole) then
693	<	if (j_is_SplitDipole) then
694	<	BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
695	<	else
696	<	BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
697	<	endif
698	<	ri = 1.0_dp / BigR
699	<	scale = rij * ri
692	>	pref = pre12 * q_i * mu_j
693	>
694	>	if (summationMethod .eq. UNDAMPED_WOLF) then
695	>	ri2 = riji * riji
696	>	ri3 = ri2 * riji
697	>
698	>	pref = pre12 * q_i * mu_j
699	>	vterm = - pref * ct_j * (ri2 - rcuti2)
700	>	vpair = vpair + vterm
701	>	epot = epot + sw*vterm
702	>
703	>	!! this has a + sign in the () because the rij vector is
704	>	!! r_j - r_i and the charge-dipole potential takes the origin
705	>	!! as the point dipole, which is atom j in this case.
706	>
707	>	dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
708	>	- rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
709	>	dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
710	>	- rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
711	>	dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
712	>	- rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
713	>
714	>	duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
715	>	duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
716	>	duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
717	>
718		else
719		if (j_is_SplitDipole) then
720		BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
721		ri = 1.0_dp / BigR
722	<	scale = rij * ri
723	<	else
722	>	scale = rij * ri
723	>	else
724		ri = riji
725		scale = 1.0_dp
726		endif
727	<	endif
727	>
728	>	ri2 = ri * ri
729	>	ri3 = ri2 * ri
730	>	sc2 = scale * scale
731
732	<	ct_ij = ul_i(1)*ul_j(1) + ul_i(2)*ul_j(2) + ul_i(3)*ul_j(3)
732	>	pref = pre12 * q_i * mu_j
733	>	vterm = - pref * ct_j * ri2 * scale
734	>	vpair = vpair + vterm
735	>	epot = epot + sw*vterm
736	>
737	>	!! this has a + sign in the () because the rij vector is
738	>	!! r_j - r_i and the charge-dipole potential takes the origin
739	>	!! as the point dipole, which is atom j in this case.
740	>
741	>	dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
742	>	dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
743	>	dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
744	>
745	>	duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
746	>	duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
747	>	duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
748
749	<	ri2 = ri * ri
750	<	ri3 = ri2 * ri
749	>	endif
750	>	endif
751	>
752	>	if (j_is_Quadrupole) then
753	>	ri2 = riji * riji
754	>	ri3 = ri2 * riji
755		ri4 = ri2 * ri2
756	<	sc2 = scale * scale
756	>	cx2 = cx_j * cx_j
757	>	cy2 = cy_j * cy_j
758	>	cz2 = cz_j * cz_j
759
760	<	pref = pre22 * mu_i * mu_j
761	<	vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
762	<	vpair = vpair + vterm
763	<	epot = epot + sw * vterm
764	<
765	<	a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
766	<
767	<	dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*ul_j(1)-ct_j*ul_i(1))
768	<	dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*ul_j(2)-ct_j*ul_i(2))
769	<	dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*ul_j(3)-ct_j*ul_i(3))
770	<
771	<	duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*xhat*sc2)
772	<	duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*yhat*sc2)
773	<	duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*zhat*sc2)
774	<
775	<	dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*xhat*sc2)
776	<	dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*yhat*sc2)
777	<	dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*zhat*sc2)
760	>	if (summationMethod .eq. UNDAMPED_WOLF) then
761	>	pref =  pre14 * q_i / 3.0_dp
762	>	vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
763	>	qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
764	>	qzz_j * (3.0_dp*cz2 - 1.0_dp) )
765	>	vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
766	>	qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
767	>	qzz_j * (3.0_dp*cz2 - 1.0_dp) )
768	>	vpair = vpair + ( vterm1 - vterm2 )
769	>	epot = epot + sw*( vterm1 - vterm2 )
770	>
771	>	dudx = dudx - (5.0_dp * &
772	>	(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
773	>	(ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
774	>	qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
775	>	(ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
776	>	qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
777	>	(ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
778	>	qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
779	>	dudy = dudy - (5.0_dp * &
780	>	(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
781	>	(ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
782	>	qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
783	>	(ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
784	>	qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
785	>	(ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
786	>	qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
787	>	dudz = dudz - (5.0_dp * &
788	>	(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
789	>	(ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
790	>	qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
791	>	(ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
792	>	qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
793	>	(ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
794	>	qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
795	>
796	>	dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
797	>	rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
798	>	dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
799	>	rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
800	>	dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
801	>	rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
802	>
803	>	duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
804	>	rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
805	>	duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
806	>	rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
807	>	duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
808	>	rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
809	>
810	>	duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
811	>	rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
812	>	duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
813	>	rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
814	>	duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
815	>	rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
816	>
817	>	else
818	>	pref =  pre14 * q_i / 3.0_dp
819	>	vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
820	>	qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
821	>	qzz_j * (3.0_dp*cz2 - 1.0_dp))
822	>	vpair = vpair + vterm
823	>	epot = epot + sw*vterm
824	>
825	>	dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
826	>	qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
827	>	qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
828	>	qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
829	>	dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
830	>	qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
831	>	qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
832	>	qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
833	>	dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
834	>	qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
835	>	qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
836	>	qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
837	>
838	>	dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
839	>	dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
840	>	dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
841	>
842	>	duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
843	>	duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
844	>	duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
845	>
846	>	duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
847	>	duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
848	>	duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
849	>
850	>	endif
851		endif
852	+	endif
853
854	+	if (i_is_Dipole) then
855	+
856	+	if (j_is_Charge) then
857	+
858	+	pref = pre12 * q_j * mu_i
859	+
860	+	if (summationMethod .eq. UNDAMPED_WOLF) then
861	+	ri2 = riji * riji
862	+	ri3 = ri2 * riji
863	+
864	+	pref = pre12 * q_j * mu_i
865	+	vterm = pref * ct_i * (ri2 - rcuti2)
866	+	vpair = vpair + vterm
867	+	epot = epot + sw*vterm
868	+
869	+	!! this has a + sign in the () because the rij vector is
870	+	!! r_j - r_i and the charge-dipole potential takes the origin
871	+	!! as the point dipole, which is atom j in this case.
872	+
873	+	dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
874	+	- rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
875	+	dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
876	+	- rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
877	+	dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
878	+	- rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
879	+
880	+	duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
881	+	duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
882	+	duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
883	+
884	+	else
885	+	if (i_is_SplitDipole) then
886	+	BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
887	+	ri = 1.0_dp / BigR
888	+	scale = rij * ri
889	+	else
890	+	ri = riji
891	+	scale = 1.0_dp
892	+	endif
893	+
894	+	ri2 = ri * ri
895	+	ri3 = ri2 * ri
896	+	sc2 = scale * scale
897	+
898	+	pref = pre12 * q_j * mu_i
899	+	vterm = pref * ct_i * ri2 * scale
900	+	vpair = vpair + vterm
901	+	epot = epot + sw*vterm
902	+
903	+	dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
904	+	dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
905	+	dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
906	+
907	+	duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
908	+	duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
909	+	duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
910	+	endif
911	+	endif
912	+
913	+	if (j_is_Dipole) then
914	+
915	+	if (summationMethod .eq. UNDAMPED_WOLF) then
916	+	ri2 = riji * riji
917	+	ri3 = ri2 * riji
918	+	ri4 = ri2 * ri2
919	+
920	+	pref = pre22 * mu_i * mu_j
921	+	vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
922	+	vpair = vpair + vterm
923	+	epot = epot + sw*vterm
924	+
925	+	a1 = 5.0d0 * ct_i * ct_j - ct_ij
926	+
927	+	dudx = dudx + sw*pref*3.0d0*ri4 &
928	+	* (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
929	+	- sw*pref*3.0d0*rcuti4 &
930	+	* (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
931	+	dudy = dudy + sw*pref*3.0d0*ri4 &
932	+	* (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
933	+	- sw*pref*3.0d0*rcuti4 &
934	+	* (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
935	+	dudz = dudz + sw*pref*3.0d0*ri4 &
936	+	* (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
937	+	- sw*pref*3.0d0*rcuti4 &
938	+	* (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
939	+
940	+	duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
941	+	- rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
942	+	duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
943	+	- rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
944	+	duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
945	+	- rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
946	+	duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
947	+	- rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
948	+	duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
949	+	- rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
950	+	duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
951	+	- rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
952	+
953	+	else
954	+	if (i_is_SplitDipole) then
955	+	if (j_is_SplitDipole) then
956	+	BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
957	+	else
958	+	BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
959	+	endif
960	+	ri = 1.0_dp / BigR
961	+	scale = rij * ri
962	+	else
963	+	if (j_is_SplitDipole) then
964	+	BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
965	+	ri = 1.0_dp / BigR
966	+	scale = rij * ri
967	+	else
968	+	ri = riji
969	+	scale = 1.0_dp
970	+	endif
971	+	endif
972	+
973	+	ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
974	+
975	+	ri2 = ri * ri
976	+	ri3 = ri2 * ri
977	+	ri4 = ri2 * ri2
978	+	sc2 = scale * scale
979	+
980	+	pref = pre22 * mu_i * mu_j
981	+	vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
982	+	vpair = vpair + vterm
983	+	epot = epot + sw*vterm
984	+
985	+	a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
986	+
987	+	dudx = dudx + sw*pref*3.0d0*ri4*scale &
988	+	*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989	+	dudy = dudy + sw*pref*3.0d0*ri4*scale &
990	+	*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991	+	dudz = dudz + sw*pref*3.0d0*ri4*scale &
992	+	*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
993	+
994	+	duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
995	+	*(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
996	+	duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
997	+	*(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
998	+	duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
999	+	*(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1000	+
1001	+	duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1002	+	*(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1003	+	duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1004	+	*(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1005	+	duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1006	+	*(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1007	+	endif
1008	+	endif
1009		endif
1010	<
1010	>
1011	>	if (i_is_Quadrupole) then
1012	>	if (j_is_Charge) then
1013	>
1014	>	ri2 = riji * riji
1015	>	ri3 = ri2 * riji
1016	>	ri4 = ri2 * ri2
1017	>	cx2 = cx_i * cx_i
1018	>	cy2 = cy_i * cy_i
1019	>	cz2 = cz_i * cz_i
1020	>
1021	>	if (summationMethod .eq. UNDAMPED_WOLF) then
1022	>	pref = pre14 * q_j / 3.0_dp
1023	>	vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1024	>	qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1025	>	qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1026	>	vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1027	>	qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1028	>	qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1029	>	vpair = vpair + ( vterm1 - vterm2 )
1030	>	epot = epot + sw*( vterm1 - vterm2 )
1031	>
1032	>	dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1033	>	sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1034	>	qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1035	>	(ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1036	>	qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1037	>	(ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1038	>	qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1039	>	dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1040	>	sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1041	>	qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1042	>	(ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1043	>	qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1044	>	(ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1045	>	qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1046	>	dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1047	>	sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1048	>	qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1049	>	(ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1050	>	qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1051	>	(ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1052	>	qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1053	>
1054	>	dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1055	>	rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1056	>	dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1057	>	rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1058	>	dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1059	>	rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1060	>
1061	>	duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1062	>	rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1063	>	duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1064	>	rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1065	>	duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1066	>	rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1067	>
1068	>	duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1069	>	rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1070	>	duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1071	>	rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1072	>	duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1073	>	rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1074	>
1075	>	else
1076	>	pref = pre14 * q_j / 3.0_dp
1077	>	vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1078	>	qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1079	>	qzz_i * (3.0_dp*cz2 - 1.0_dp))
1080	>	vpair = vpair + vterm
1081	>	epot = epot + sw*vterm
1082	>
1083	>	dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1084	>	qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1085	>	qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1086	>	qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1087	>	dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1088	>	qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1089	>	qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1090	>	qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1091	>	dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1092	>	qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1093	>	qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1094	>	qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1095	>
1096	>	dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1097	>	dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1098	>	dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1099	>
1100	>	duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1101	>	duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1102	>	duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1103	>
1104	>	duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1105	>	duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1106	>	duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1107	>	endif
1108	>	endif
1109	>	endif
1110	>
1111	>
1112		if (do_pot) then
1113		#ifdef IS_MPI
1114	<	pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1115	<	pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1114	>	pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1115	>	pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1116		#else
1117		pot = pot + epot
1118		#endif
1119		endif
1120	<
1120	>
1121		#ifdef IS_MPI
1122		f_Row(1,atom1) = f_Row(1,atom1) + dudx
1123		f_Row(2,atom1) = f_Row(2,atom1) + dudy
1124		f_Row(3,atom1) = f_Row(3,atom1) + dudz
1125	<
1125	>
1126		f_Col(1,atom2) = f_Col(1,atom2) - dudx
1127		f_Col(2,atom2) = f_Col(2,atom2) - dudy
1128		f_Col(3,atom2) = f_Col(3,atom2) - dudz
1129	<
1129	>
1130		if (i_is_Dipole .or. i_is_Quadrupole) then
1131	<	t_Row(1,atom1) = t_Row(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy
1132	<	t_Row(2,atom1) = t_Row(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz
1133	<	t_Row(3,atom1) = t_Row(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix
1131	>	t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1132	>	t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1133	>	t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1134		endif
1135	+	if (i_is_Quadrupole) then
1136	+	t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1137	+	t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1138	+	t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1139
1140	+	t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1141	+	t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1142	+	t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1143	+	endif
1144	+
1145		if (j_is_Dipole .or. j_is_Quadrupole) then
1146	<	t_Col(1,atom2) = t_Col(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy
1147	<	t_Col(2,atom2) = t_Col(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz
1148	<	t_Col(3,atom2) = t_Col(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx
1146	>	t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1147	>	t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1148	>	t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1149		endif
1150	+	if (j_is_Quadrupole) then
1151	+	t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1152	+	t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1153	+	t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1154
1155	+	t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1156	+	t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1157	+	t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1158	+	endif
1159	+
1160		#else
1161		f(1,atom1) = f(1,atom1) + dudx
1162		f(2,atom1) = f(2,atom1) + dudy
1163		f(3,atom1) = f(3,atom1) + dudz
1164	<
1164	>
1165		f(1,atom2) = f(1,atom2) - dudx
1166		f(2,atom2) = f(2,atom2) - dudy
1167		f(3,atom2) = f(3,atom2) - dudz
1168	<
1168	>
1169		if (i_is_Dipole .or. i_is_Quadrupole) then
1170	<	t(1,atom1) = t(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy
1171	<	t(2,atom1) = t(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz
1172	<	t(3,atom1) = t(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix
1170	>	t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1171	>	t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1172	>	t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1173		endif
1174	<
1174	>	if (i_is_Quadrupole) then
1175	>	t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1176	>	t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1177	>	t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1178	>
1179	>	t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1180	>	t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1181	>	t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1182	>	endif
1183	>
1184		if (j_is_Dipole .or. j_is_Quadrupole) then
1185	<	t(1,atom2) = t(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy
1186	<	t(2,atom2) = t(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz
1187	<	t(3,atom2) = t(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx
1185	>	t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1186	>	t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1187	>	t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1188		endif
1189	+	if (j_is_Quadrupole) then
1190	+	t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1191	+	t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1192	+	t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1193	+
1194	+	t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1195	+	t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1196	+	t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1197	+	endif
1198	+
1199		#endif
1200	<
1200	>
1201		#ifdef IS_MPI
1202		id1 = AtomRowToGlobal(atom1)
1203		id2 = AtomColToGlobal(atom2)
#	Line 624 | Line 1207 | contains
1207		#endif
1208
1209		if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1210	<
1210	>
1211		fpair(1) = fpair(1) + dudx
1212		fpair(2) = fpair(2) + dudy
1213		fpair(3) = fpair(3) + dudz
#	Line 633 | Line 1216 | contains
1216
1217		return
1218		end subroutine doElectrostaticPair
1219	<
1219	>
1220	>	subroutine destroyElectrostaticTypes()
1221	>
1222	>	if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1223	>
1224	>	end subroutine destroyElectrostaticTypes
1225	>
1226	>	subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
1227	>
1228	>	integer, intent(in) :: atom1, atom2
1229	>	real (kind = dp), intent(in) :: rij
1230	>	real (kind = dp), dimension(9,nLocal) :: eFrame
1231	>
1232	>	integer :: me1, me2
1233	>	real (kind = dp), intent(in) :: taper
1234	>	real (kind = dp):: mu1, mu2
1235	>	real (kind = dp), dimension(3) :: ul1
1236	>	real (kind = dp), dimension(3) :: ul2
1237	>
1238	>	integer :: localError
1239	>
1240	>	#ifdef IS_MPI
1241	>	me1 = atid_Row(atom1)
1242	>	ul1(1) = eFrame_Row(3,atom1)
1243	>	ul1(2) = eFrame_Row(6,atom1)
1244	>	ul1(3) = eFrame_Row(9,atom1)
1245	>
1246	>	me2 = atid_Col(atom2)
1247	>	ul2(1) = eFrame_Col(3,atom2)
1248	>	ul2(2) = eFrame_Col(6,atom2)
1249	>	ul2(3) = eFrame_Col(9,atom2)
1250	>	#else
1251	>	me1 = atid(atom1)
1252	>	ul1(1) = eFrame(3,atom1)
1253	>	ul1(2) = eFrame(6,atom1)
1254	>	ul1(3) = eFrame(9,atom1)
1255	>
1256	>	me2 = atid(atom2)
1257	>	ul2(1) = eFrame(3,atom2)
1258	>	ul2(2) = eFrame(6,atom2)
1259	>	ul2(3) = eFrame(9,atom2)
1260	>	#endif
1261	>
1262	>	mu1 = getDipoleMoment(me1)
1263	>	mu2 = getDipoleMoment(me2)
1264	>
1265	>	#ifdef IS_MPI
1266	>	rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
1267	>	rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
1268	>	rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
1269	>
1270	>	rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
1271	>	rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
1272	>	rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
1273	>	#else
1274	>	rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
1275	>	rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
1276	>	rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
1277	>
1278	>	rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
1279	>	rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
1280	>	rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper
1281	>	#endif
1282	>	return
1283	>	end subroutine accumulate_rf
1284	>
1285	>	subroutine accumulate_self_rf(atom1, mu1, eFrame)
1286	>
1287	>	integer, intent(in) :: atom1
1288	>	real(kind=dp), intent(in) :: mu1
1289	>	real(kind=dp), dimension(9,nLocal) :: eFrame
1290	>
1291	>	!! should work for both MPI and non-MPI version since this is not pairwise.
1292	>	rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
1293	>	rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
1294	>	rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
1295	>
1296	>	return
1297	>	end subroutine accumulate_self_rf
1298	>
1299	>	subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
1300	>
1301	>	integer, intent(in) :: a1
1302	>	real (kind=dp), intent(in) :: mu1
1303	>	real (kind=dp), intent(inout) :: rfpot
1304	>	logical, intent(in) :: do_pot
1305	>	real (kind = dp), dimension(9,nLocal) :: eFrame
1306	>	real (kind = dp), dimension(3,nLocal) :: t
1307	>
1308	>	integer :: localError
1309	>
1310	>	if (.not.preRFCalculated) then
1311	>	call setReactionFieldPrefactor()
1312	>	endif
1313	>
1314	>	! compute torques on dipoles:
1315	>	! pre converts from mu in units of debye to kcal/mol
1316	>
1317	>	! The torque contribution is dipole cross reaction_field
1318	>
1319	>	t(1,a1) = t(1,a1) + preRF*mu1*(eFrame(6,a1)*rf(3,a1) - &
1320	>	eFrame(9,a1)*rf(2,a1))
1321	>	t(2,a1) = t(2,a1) + preRF*mu1*(eFrame(9,a1)*rf(1,a1) - &
1322	>	eFrame(3,a1)*rf(3,a1))
1323	>	t(3,a1) = t(3,a1) + preRF*mu1*(eFrame(3,a1)*rf(2,a1) - &
1324	>	eFrame(6,a1)*rf(1,a1))
1325	>
1326	>	! the potential contribution is -1/2 dipole dot reaction_field
1327	>
1328	>	if (do_pot) then
1329	>	rfpot = rfpot - 0.5d0 * preRF * mu1 * &
1330	>	(rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + &
1331	>	rf(3,a1)*eFrame(9,a1))
1332	>	endif
1333	>
1334	>	return
1335	>	end subroutine reaction_field_final
1336	>
1337	>	subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
1338	>
1339	>	integer, intent(in) :: atom1, atom2
1340	>	real(kind=dp), dimension(3), intent(in) :: d
1341	>	real(kind=dp), intent(in) :: rij, taper
1342	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1343	>	real( kind = dp ), dimension(3,nLocal) :: f
1344	>	real( kind = dp ), dimension(3), intent(inout) :: fpair
1345	>
1346	>	real (kind = dp), dimension(3) :: ul1
1347	>	real (kind = dp), dimension(3) :: ul2
1348	>	real (kind = dp) :: dtdr
1349	>	real (kind = dp) :: dudx, dudy, dudz, u1dotu2
1350	>	integer :: me1, me2, id1, id2
1351	>	real (kind = dp) :: mu1, mu2
1352	>
1353	>	integer :: localError
1354	>
1355	>	if (.not.preRFCalculated) then
1356	>	call setReactionFieldPrefactor()
1357	>	endif
1358	>
1359	>	if (rij.le.rrf) then
1360	>
1361	>	if (rij.lt.rt) then
1362	>	dtdr = 0.0d0
1363	>	else
1364	>	!         write(*,*) 'rf correct in taper region'
1365	>	dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
1366	>	endif
1367	>
1368	>	#ifdef IS_MPI
1369	>	me1 = atid_Row(atom1)
1370	>	ul1(1) = eFrame_Row(3,atom1)
1371	>	ul1(2) = eFrame_Row(6,atom1)
1372	>	ul1(3) = eFrame_Row(9,atom1)
1373	>
1374	>	me2 = atid_Col(atom2)
1375	>	ul2(1) = eFrame_Col(3,atom2)
1376	>	ul2(2) = eFrame_Col(6,atom2)
1377	>	ul2(3) = eFrame_Col(9,atom2)
1378	>	#else
1379	>	me1 = atid(atom1)
1380	>	ul1(1) = eFrame(3,atom1)
1381	>	ul1(2) = eFrame(6,atom1)
1382	>	ul1(3) = eFrame(9,atom1)
1383	>
1384	>	me2 = atid(atom2)
1385	>	ul2(1) = eFrame(3,atom2)
1386	>	ul2(2) = eFrame(6,atom2)
1387	>	ul2(3) = eFrame(9,atom2)
1388	>	#endif
1389	>
1390	>	mu1 = getDipoleMoment(me1)
1391	>	mu2 = getDipoleMoment(me2)
1392	>
1393	>	u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
1394	>
1395	>	dudx = - preRF*mu1*mu2*u1dotu2*dtdr*d(1)/rij
1396	>	dudy = - preRF*mu1*mu2*u1dotu2*dtdr*d(2)/rij
1397	>	dudz = - preRF*mu1*mu2*u1dotu2*dtdr*d(3)/rij
1398	>
1399	>	#ifdef IS_MPI
1400	>	f_Row(1,atom1) = f_Row(1,atom1) + dudx
1401	>	f_Row(2,atom1) = f_Row(2,atom1) + dudy
1402	>	f_Row(3,atom1) = f_Row(3,atom1) + dudz
1403	>
1404	>	f_Col(1,atom2) = f_Col(1,atom2) - dudx
1405	>	f_Col(2,atom2) = f_Col(2,atom2) - dudy
1406	>	f_Col(3,atom2) = f_Col(3,atom2) - dudz
1407	>	#else
1408	>	f(1,atom1) = f(1,atom1) + dudx
1409	>	f(2,atom1) = f(2,atom1) + dudy
1410	>	f(3,atom1) = f(3,atom1) + dudz
1411	>
1412	>	f(1,atom2) = f(1,atom2) - dudx
1413	>	f(2,atom2) = f(2,atom2) - dudy
1414	>	f(3,atom2) = f(3,atom2) - dudz
1415	>	#endif
1416	>
1417	>	#ifdef IS_MPI
1418	>	id1 = AtomRowToGlobal(atom1)
1419	>	id2 = AtomColToGlobal(atom2)
1420	>	#else
1421	>	id1 = atom1
1422	>	id2 = atom2
1423	>	#endif
1424	>
1425	>	if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1426	>
1427	>	fpair(1) = fpair(1) + dudx
1428	>	fpair(2) = fpair(2) + dudy
1429	>	fpair(3) = fpair(3) + dudz
1430	>
1431	>	endif
1432	>
1433	>	end if
1434	>	return
1435	>	end subroutine rf_correct_forces
1436	>
1437		end module electrostatic_module

Diff Legend

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