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

Diff Legend

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