ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 986
Committed: Wed Jun 7 22:49:26 2006 UTC (19 years, 4 months ago) by chrisfen
File size: 49419 byte(s)
Log Message: 
damping now works for charge-quadrupole and dipole-dipole

File Contents

# User Rev Content
1 gezelter 411 !!
2     !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     !!
4     !! The University of Notre Dame grants you ("Licensee") a
5     !! non-exclusive, royalty free, license to use, modify and
6     !! redistribute this software in source and binary code form, provided
7     !! that the following conditions are met:
8     !!
9     !! 1. Acknowledgement of the program authors must be made in any
10     !! publication of scientific results based in part on use of the
11     !! program. An acceptable form of acknowledgement is citation of
12     !! the article in which the program was described (Matthew
13     !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     !! Parallel Simulation Engine for Molecular Dynamics,"
16     !! J. Comput. Chem. 26, pp. 252-271 (2005))
17     !!
18     !! 2. Redistributions of source code must retain the above copyright
19     !! notice, this list of conditions and the following disclaimer.
20     !!
21     !! 3. Redistributions in binary form must reproduce the above copyright
22     !! notice, this list of conditions and the following disclaimer in the
23     !! documentation and/or other materials provided with the
24     !! distribution.
25     !!
26     !! This software is provided "AS IS," without a warranty of any
27     !! kind. All express or implied conditions, representations and
28     !! warranties, including any implied warranty of merchantability,
29     !! fitness for a particular purpose or non-infringement, are hereby
30     !! excluded. The University of Notre Dame and its licensors shall not
31     !! be liable for any damages suffered by licensee as a result of
32     !! using, modifying or distributing the software or its
33     !! derivatives. In no event will the University of Notre Dame or its
34     !! licensors be liable for any lost revenue, profit or data, or for
35     !! direct, indirect, special, consequential, incidental or punitive
36     !! damages, however caused and regardless of the theory of liability,
37     !! arising out of the use of or inability to use software, even if the
38     !! University of Notre Dame has been advised of the possibility of
39     !! such damages.
40     !!
41    
42     module electrostatic_module
43 gezelter 507
44 gezelter 411 use force_globals
45     use definitions
46     use atype_module
47     use vector_class
48     use simulation
49     use status
50 chrisfen 937 use interpolation
51 gezelter 411 #ifdef IS_MPI
52     use mpiSimulation
53     #endif
54     implicit none
55    
56     PRIVATE
57    
58 chuckv 656
59 gezelter 602 #define __FORTRAN90
60 chuckv 656 #include "UseTheForce/DarkSide/fInteractionMap.h"
61 gezelter 602 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
62 chrisfen 716 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
63 gezelter 602
64 chuckv 656
65 gezelter 434 !! these prefactors convert the multipole interactions into kcal / mol
66     !! all were computed assuming distances are measured in angstroms
67     !! Charge-Charge, assuming charges are measured in electrons
68 chrisfen 959 real(kind=dp), parameter :: pre11 = 332.0637778_dp
69 gezelter 434 !! Charge-Dipole, assuming charges are measured in electrons, and
70     !! dipoles are measured in debyes
71 chrisfen 959 real(kind=dp), parameter :: pre12 = 69.13373_dp
72 gezelter 434 !! Dipole-Dipole, assuming dipoles are measured in debyes
73 chrisfen 959 real(kind=dp), parameter :: pre22 = 14.39325_dp
74 gezelter 434 !! Charge-Quadrupole, assuming charges are measured in electrons, and
75     !! quadrupoles are measured in 10^-26 esu cm^2
76     !! This unit is also known affectionately as an esu centi-barn.
77 chrisfen 959 real(kind=dp), parameter :: pre14 = 69.13373_dp
78 gezelter 411
79 chrisfen 959 real(kind=dp), parameter :: zero = 0.0_dp
80 chrisfen 941
81     !! number of points for electrostatic splines
82     integer, parameter :: np = 100
83 gezelter 939
84 chrisfen 712 !! variables to handle different summation methods for long-range
85     !! electrostatics:
86 gezelter 602 integer, save :: summationMethod = NONE
87 chrisfen 710 integer, save :: screeningMethod = UNDAMPED
88 chrisfen 603 logical, save :: summationMethodChecked = .false.
89 gezelter 602 real(kind=DP), save :: defaultCutoff = 0.0_DP
90 chrisfen 682 real(kind=DP), save :: defaultCutoff2 = 0.0_DP
91 gezelter 602 logical, save :: haveDefaultCutoff = .false.
92     real(kind=DP), save :: dampingAlpha = 0.0_DP
93 chrisfen 716 real(kind=DP), save :: alpha2 = 0.0_DP
94 gezelter 602 logical, save :: haveDampingAlpha = .false.
95 chrisfen 682 real(kind=DP), save :: dielectric = 1.0_DP
96 gezelter 602 logical, save :: haveDielectric = .false.
97     real(kind=DP), save :: constEXP = 0.0_DP
98 chrisfen 682 real(kind=dp), save :: rcuti = 0.0_DP
99     real(kind=dp), save :: rcuti2 = 0.0_DP
100     real(kind=dp), save :: rcuti3 = 0.0_DP
101     real(kind=dp), save :: rcuti4 = 0.0_DP
102     real(kind=dp), save :: alphaPi = 0.0_DP
103     real(kind=dp), save :: invRootPi = 0.0_DP
104     real(kind=dp), save :: rrf = 1.0_DP
105     real(kind=dp), save :: rt = 1.0_DP
106     real(kind=dp), save :: rrfsq = 1.0_DP
107     real(kind=dp), save :: preRF = 0.0_DP
108 chrisfen 695 real(kind=dp), save :: preRF2 = 0.0_DP
109 chrisfen 716 real(kind=dp), save :: f0 = 1.0_DP
110     real(kind=dp), save :: f1 = 1.0_DP
111     real(kind=dp), save :: f2 = 0.0_DP
112 chrisfen 849 real(kind=dp), save :: f3 = 0.0_DP
113     real(kind=dp), save :: f4 = 0.0_DP
114 chrisfen 716 real(kind=dp), save :: f0c = 1.0_DP
115     real(kind=dp), save :: f1c = 1.0_DP
116     real(kind=dp), save :: f2c = 0.0_DP
117 chrisfen 849 real(kind=dp), save :: f3c = 0.0_DP
118     real(kind=dp), save :: f4c = 0.0_DP
119 chrisfen 941 real(kind=dp), save :: df0 = 0.0_DP
120     type(cubicSpline), save :: f0spline
121     logical, save :: haveElectroSpline = .false.
122 chrisfen 959 real(kind=dp), save :: one_third = 1.0_DP / 3.0_DP
123 chrisfen 716
124 gezelter 809 #if defined(__IFC) || defined(__PGI)
125 chuckv 632 ! error function for ifc version > 7.
126 gezelter 960 real(kind=dp), external :: erfc
127 chuckv 632 #endif
128    
129 chrisfen 853 public :: setElectrostaticSummationMethod
130 chrisfen 710 public :: setScreeningMethod
131 gezelter 602 public :: setElectrostaticCutoffRadius
132 chrisfen 710 public :: setDampingAlpha
133 gezelter 602 public :: setReactionFieldDielectric
134 chrisfen 941 public :: buildElectroSpline
135 gezelter 411 public :: newElectrostaticType
136     public :: setCharge
137     public :: setDipoleMoment
138     public :: setSplitDipoleDistance
139     public :: setQuadrupoleMoments
140     public :: doElectrostaticPair
141     public :: getCharge
142     public :: getDipoleMoment
143 chuckv 492 public :: destroyElectrostaticTypes
144 chrisfen 703 public :: self_self
145 chrisfen 700 public :: rf_self_excludes
146 gezelter 411
147 chrisfen 937
148 gezelter 411 type :: Electrostatic
149     integer :: c_ident
150     logical :: is_Charge = .false.
151     logical :: is_Dipole = .false.
152     logical :: is_SplitDipole = .false.
153     logical :: is_Quadrupole = .false.
154 chrisfen 532 logical :: is_Tap = .false.
155 gezelter 411 real(kind=DP) :: charge = 0.0_DP
156     real(kind=DP) :: dipole_moment = 0.0_DP
157     real(kind=DP) :: split_dipole_distance = 0.0_DP
158     real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP
159     end type Electrostatic
160    
161     type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
162    
163 gezelter 938 logical, save :: hasElectrostaticMap
164    
165 gezelter 411 contains
166    
167 chrisfen 853 subroutine setElectrostaticSummationMethod(the_ESM)
168 gezelter 602 integer, intent(in) :: the_ESM
169    
170     if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
171 chrisfen 853 call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
172 gezelter 602 endif
173    
174 chrisfen 610 summationMethod = the_ESM
175 chrisfen 626
176 chrisfen 853 end subroutine setElectrostaticSummationMethod
177 gezelter 602
178 chrisfen 710 subroutine setScreeningMethod(the_SM)
179     integer, intent(in) :: the_SM
180     screeningMethod = the_SM
181     end subroutine setScreeningMethod
182    
183 chrisfen 682 subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
184 gezelter 602 real(kind=dp), intent(in) :: thisRcut
185 chrisfen 682 real(kind=dp), intent(in) :: thisRsw
186 gezelter 602 defaultCutoff = thisRcut
187 chrisfen 849 defaultCutoff2 = defaultCutoff*defaultCutoff
188 chrisfen 682 rrf = defaultCutoff
189     rt = thisRsw
190 gezelter 602 haveDefaultCutoff = .true.
191     end subroutine setElectrostaticCutoffRadius
192    
193 chrisfen 710 subroutine setDampingAlpha(thisAlpha)
194 gezelter 602 real(kind=dp), intent(in) :: thisAlpha
195     dampingAlpha = thisAlpha
196 chrisfen 716 alpha2 = dampingAlpha*dampingAlpha
197 gezelter 602 haveDampingAlpha = .true.
198 chrisfen 710 end subroutine setDampingAlpha
199 gezelter 602
200     subroutine setReactionFieldDielectric(thisDielectric)
201     real(kind=dp), intent(in) :: thisDielectric
202     dielectric = thisDielectric
203     haveDielectric = .true.
204     end subroutine setReactionFieldDielectric
205    
206 chrisfen 941 subroutine buildElectroSpline()
207     real( kind = dp ), dimension(np) :: xvals, yvals
208     real( kind = dp ) :: dx, rmin, rval
209     integer :: i
210 chrisfen 937
211 chrisfen 959 rmin = 0.0_dp
212 chrisfen 941
213     dx = (defaultCutoff-rmin) / dble(np-1)
214    
215     do i = 1, np
216     rval = rmin + dble(i-1)*dx
217     xvals(i) = rval
218 gezelter 960 yvals(i) = erfc(dampingAlpha*rval)
219 chrisfen 941 enddo
220    
221     call newSpline(f0spline, xvals, yvals, .true.)
222    
223     haveElectroSpline = .true.
224     end subroutine buildElectroSpline
225    
226 gezelter 411 subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
227 chrisfen 532 is_SplitDipole, is_Quadrupole, is_Tap, status)
228 gezelter 507
229 gezelter 411 integer, intent(in) :: c_ident
230     logical, intent(in) :: is_Charge
231     logical, intent(in) :: is_Dipole
232     logical, intent(in) :: is_SplitDipole
233     logical, intent(in) :: is_Quadrupole
234 chrisfen 532 logical, intent(in) :: is_Tap
235 gezelter 411 integer, intent(out) :: status
236     integer :: nAtypes, myATID, i, j
237    
238     status = 0
239     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
240 gezelter 507
241 gezelter 411 !! Be simple-minded and assume that we need an ElectrostaticMap that
242     !! is the same size as the total number of atom types
243    
244     if (.not.allocated(ElectrostaticMap)) then
245 gezelter 507
246 gezelter 411 nAtypes = getSize(atypes)
247 gezelter 507
248 gezelter 411 if (nAtypes == 0) then
249     status = -1
250     return
251     end if
252 gezelter 507
253 gezelter 938 allocate(ElectrostaticMap(nAtypes))
254 gezelter 507
255 gezelter 411 end if
256    
257     if (myATID .gt. size(ElectrostaticMap)) then
258     status = -1
259     return
260     endif
261 gezelter 507
262 gezelter 411 ! set the values for ElectrostaticMap for this atom type:
263    
264     ElectrostaticMap(myATID)%c_ident = c_ident
265     ElectrostaticMap(myATID)%is_Charge = is_Charge
266     ElectrostaticMap(myATID)%is_Dipole = is_Dipole
267     ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
268     ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
269 chrisfen 532 ElectrostaticMap(myATID)%is_Tap = is_Tap
270 gezelter 507
271 gezelter 938 hasElectrostaticMap = .true.
272    
273 gezelter 411 end subroutine newElectrostaticType
274    
275     subroutine setCharge(c_ident, charge, status)
276     integer, intent(in) :: c_ident
277     real(kind=dp), intent(in) :: charge
278     integer, intent(out) :: status
279     integer :: myATID
280    
281     status = 0
282     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
283    
284 gezelter 938 if (.not.hasElectrostaticMap) then
285 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!")
286     status = -1
287     return
288     end if
289    
290     if (myATID .gt. size(ElectrostaticMap)) then
291     call handleError("electrostatic", "ElectrostaticMap was found to be too small during setCharge!")
292     status = -1
293     return
294     endif
295    
296     if (.not.ElectrostaticMap(myATID)%is_Charge) then
297     call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
298     status = -1
299     return
300 gezelter 507 endif
301 gezelter 411
302     ElectrostaticMap(myATID)%charge = charge
303     end subroutine setCharge
304    
305     subroutine setDipoleMoment(c_ident, dipole_moment, status)
306     integer, intent(in) :: c_ident
307     real(kind=dp), intent(in) :: dipole_moment
308     integer, intent(out) :: status
309     integer :: myATID
310    
311     status = 0
312     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
313    
314 gezelter 938 if (.not.hasElectrostaticMap) then
315 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!")
316     status = -1
317     return
318     end if
319    
320     if (myATID .gt. size(ElectrostaticMap)) then
321     call handleError("electrostatic", "ElectrostaticMap was found to be too small during setDipoleMoment!")
322     status = -1
323     return
324     endif
325    
326     if (.not.ElectrostaticMap(myATID)%is_Dipole) then
327     call handleError("electrostatic", "Attempt to setDipoleMoment of an atom type that is not a dipole!")
328     status = -1
329     return
330     endif
331    
332     ElectrostaticMap(myATID)%dipole_moment = dipole_moment
333     end subroutine setDipoleMoment
334    
335     subroutine setSplitDipoleDistance(c_ident, split_dipole_distance, status)
336     integer, intent(in) :: c_ident
337     real(kind=dp), intent(in) :: split_dipole_distance
338     integer, intent(out) :: status
339     integer :: myATID
340    
341     status = 0
342     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
343    
344 gezelter 938 if (.not.hasElectrostaticMap) then
345 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!")
346     status = -1
347     return
348     end if
349    
350     if (myATID .gt. size(ElectrostaticMap)) then
351     call handleError("electrostatic", "ElectrostaticMap was found to be too small during setSplitDipoleDistance!")
352     status = -1
353     return
354     endif
355    
356     if (.not.ElectrostaticMap(myATID)%is_SplitDipole) then
357     call handleError("electrostatic", "Attempt to setSplitDipoleDistance of an atom type that is not a splitDipole!")
358     status = -1
359     return
360     endif
361    
362     ElectrostaticMap(myATID)%split_dipole_distance = split_dipole_distance
363     end subroutine setSplitDipoleDistance
364    
365     subroutine setQuadrupoleMoments(c_ident, quadrupole_moments, status)
366     integer, intent(in) :: c_ident
367     real(kind=dp), intent(in), dimension(3) :: quadrupole_moments
368     integer, intent(out) :: status
369     integer :: myATID, i, j
370    
371     status = 0
372     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
373    
374 gezelter 938 if (.not.hasElectrostaticMap) then
375 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
376     status = -1
377     return
378     end if
379    
380     if (myATID .gt. size(ElectrostaticMap)) then
381     call handleError("electrostatic", "ElectrostaticMap was found to be too small during setQuadrupoleMoments!")
382     status = -1
383     return
384     endif
385    
386     if (.not.ElectrostaticMap(myATID)%is_Quadrupole) then
387     call handleError("electrostatic", "Attempt to setQuadrupoleMoments of an atom type that is not a quadrupole!")
388     status = -1
389     return
390     endif
391 gezelter 507
392 gezelter 411 do i = 1, 3
393 gezelter 507 ElectrostaticMap(myATID)%quadrupole_moments(i) = &
394     quadrupole_moments(i)
395     enddo
396 gezelter 411
397     end subroutine setQuadrupoleMoments
398    
399 gezelter 507
400 gezelter 411 function getCharge(atid) result (c)
401     integer, intent(in) :: atid
402     integer :: localError
403     real(kind=dp) :: c
404 gezelter 507
405 gezelter 938 if (.not.hasElectrostaticMap) then
406 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
407     return
408     end if
409 gezelter 507
410 gezelter 411 if (.not.ElectrostaticMap(atid)%is_Charge) then
411     call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
412     return
413     endif
414 gezelter 507
415 gezelter 411 c = ElectrostaticMap(atid)%charge
416     end function getCharge
417    
418     function getDipoleMoment(atid) result (dm)
419     integer, intent(in) :: atid
420     integer :: localError
421     real(kind=dp) :: dm
422 gezelter 507
423 gezelter 938 if (.not.hasElectrostaticMap) then
424 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
425     return
426     end if
427 gezelter 507
428 gezelter 411 if (.not.ElectrostaticMap(atid)%is_Dipole) then
429     call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
430     return
431     endif
432 gezelter 507
433 gezelter 411 dm = ElectrostaticMap(atid)%dipole_moment
434     end function getDipoleMoment
435    
436 gezelter 602 subroutine checkSummationMethod()
437    
438 chrisfen 607 if (.not.haveDefaultCutoff) then
439     call handleError("checkSummationMethod", "no Default Cutoff set!")
440     endif
441    
442 chrisfen 959 rcuti = 1.0_dp / defaultCutoff
443 chrisfen 607 rcuti2 = rcuti*rcuti
444     rcuti3 = rcuti2*rcuti
445     rcuti4 = rcuti2*rcuti2
446    
447 chrisfen 710 if (screeningMethod .eq. DAMPED) then
448 chrisfen 703 if (.not.haveDampingAlpha) then
449     call handleError("checkSummationMethod", "no Damping Alpha set!")
450     endif
451    
452     if (.not.haveDefaultCutoff) then
453     call handleError("checkSummationMethod", "no Default Cutoff set!")
454     endif
455 chrisfen 603
456 chrisfen 849 constEXP = exp(-alpha2*defaultCutoff2)
457 chrisfen 959 invRootPi = 0.56418958354775628695_dp
458     alphaPi = 2.0_dp*dampingAlpha*invRootPi
459 gezelter 960 f0c = erfc(dampingAlpha*defaultCutoff)
460 chrisfen 716 f1c = alphaPi*defaultCutoff*constEXP + f0c
461 chrisfen 959 f2c = alphaPi*2.0_dp*alpha2*constEXP
462     f3c = alphaPi*2.0_dp*alpha2*constEXP*defaultCutoff2*defaultCutoff
463 gezelter 602 endif
464    
465 chrisfen 603 if (summationMethod .eq. REACTION_FIELD) then
466 chrisfen 703 if (haveDielectric) then
467     defaultCutoff2 = defaultCutoff*defaultCutoff
468 chrisfen 959 preRF = (dielectric-1.0_dp) / &
469     ((2.0_dp*dielectric+1.0_dp)*defaultCutoff2*defaultCutoff)
470     preRF2 = 2.0_dp*preRF
471 chrisfen 703 else
472     call handleError("checkSummationMethod", "Dielectric not set")
473 chrisfen 603 endif
474 chrisfen 703
475 chrisfen 603 endif
476    
477 chrisfen 941 if (.not.haveElectroSpline) then
478     call buildElectroSpline()
479     end if
480    
481 chrisfen 603 summationMethodChecked = .true.
482 gezelter 602 end subroutine checkSummationMethod
483    
484 chrisfen 712
485 gezelter 762 subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, rcut, sw, &
486 chrisfen 712 vpair, fpair, pot, eFrame, f, t, do_pot)
487 gezelter 507
488 chrisfen 703 logical, intent(in) :: do_pot
489 gezelter 507
490 gezelter 411 integer, intent(in) :: atom1, atom2
491     integer :: localError
492    
493 gezelter 762 real(kind=dp), intent(in) :: rij, r2, sw, rcut
494 gezelter 411 real(kind=dp), intent(in), dimension(3) :: d
495     real(kind=dp), intent(inout) :: vpair
496 chrisfen 703 real(kind=dp), intent(inout), dimension(3) :: fpair
497 gezelter 411
498 chrisfen 626 real( kind = dp ) :: pot
499 gezelter 411 real( kind = dp ), dimension(9,nLocal) :: eFrame
500     real( kind = dp ), dimension(3,nLocal) :: f
501 chrisfen 710 real( kind = dp ), dimension(3,nLocal) :: felec
502 gezelter 411 real( kind = dp ), dimension(3,nLocal) :: t
503 gezelter 507
504 gezelter 439 real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
505     real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
506     real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
507     real (kind = dp), dimension(3) :: dudux_j, duduy_j, duduz_j
508 gezelter 411
509     logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
510     logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
511 chrisfen 532 logical :: i_is_Tap, j_is_Tap
512 gezelter 411 integer :: me1, me2, id1, id2
513     real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
514 gezelter 439 real (kind=dp) :: qxx_i, qyy_i, qzz_i
515     real (kind=dp) :: qxx_j, qyy_j, qzz_j
516     real (kind=dp) :: cx_i, cy_i, cz_i
517     real (kind=dp) :: cx_j, cy_j, cz_j
518     real (kind=dp) :: cx2, cy2, cz2
519 chrisfen 719 real (kind=dp) :: ct_i, ct_j, ct_ij, a0, a1
520 gezelter 421 real (kind=dp) :: riji, ri, ri2, ri3, ri4
521 chrisfen 597 real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
522 gezelter 421 real (kind=dp) :: xhat, yhat, zhat
523 gezelter 411 real (kind=dp) :: dudx, dudy, dudz
524 chrisfen 626 real (kind=dp) :: scale, sc2, bigR
525 chrisfen 716 real (kind=dp) :: varEXP
526 chrisfen 719 real (kind=dp) :: pot_term
527 chrisfen 695 real (kind=dp) :: preVal, rfVal
528 chrisfen 959 real (kind=dp) :: cti3, ctj3, ctidotj
529     real (kind=dp) :: ri7damp, ri5damp, prei3, prei4
530     real (kind=dp) :: xhatdot2, yhatdot2, zhatdot2
531     real (kind=dp) :: xhatdot5, yhatdot5, zhatdot5
532 gezelter 411
533 gezelter 602 if (.not.summationMethodChecked) then
534     call checkSummationMethod()
535     endif
536    
537 gezelter 411 #ifdef IS_MPI
538     me1 = atid_Row(atom1)
539     me2 = atid_Col(atom2)
540     #else
541     me1 = atid(atom1)
542     me2 = atid(atom2)
543     #endif
544    
545     !! some variables we'll need independent of electrostatic type:
546    
547 chrisfen 959 riji = 1.0_dp / rij
548 chrisfen 644
549 gezelter 421 xhat = d(1) * riji
550     yhat = d(2) * riji
551     zhat = d(3) * riji
552 gezelter 411
553     !! logicals
554     i_is_Charge = ElectrostaticMap(me1)%is_Charge
555     i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
556     i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
557     i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
558 chrisfen 532 i_is_Tap = ElectrostaticMap(me1)%is_Tap
559 gezelter 411
560     j_is_Charge = ElectrostaticMap(me2)%is_Charge
561     j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
562     j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
563     j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
564 chrisfen 532 j_is_Tap = ElectrostaticMap(me2)%is_Tap
565 gezelter 411
566     if (i_is_Charge) then
567     q_i = ElectrostaticMap(me1)%charge
568     endif
569 gezelter 507
570 gezelter 411 if (i_is_Dipole) then
571     mu_i = ElectrostaticMap(me1)%dipole_moment
572     #ifdef IS_MPI
573 gezelter 439 uz_i(1) = eFrame_Row(3,atom1)
574     uz_i(2) = eFrame_Row(6,atom1)
575     uz_i(3) = eFrame_Row(9,atom1)
576 gezelter 411 #else
577 gezelter 439 uz_i(1) = eFrame(3,atom1)
578     uz_i(2) = eFrame(6,atom1)
579     uz_i(3) = eFrame(9,atom1)
580 gezelter 411 #endif
581 gezelter 439 ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
582 gezelter 411
583     if (i_is_SplitDipole) then
584     d_i = ElectrostaticMap(me1)%split_dipole_distance
585     endif
586 gezelter 939 duduz_i = zero
587 gezelter 411 endif
588    
589 gezelter 439 if (i_is_Quadrupole) then
590     qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1)
591     qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2)
592     qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3)
593     #ifdef IS_MPI
594     ux_i(1) = eFrame_Row(1,atom1)
595     ux_i(2) = eFrame_Row(4,atom1)
596     ux_i(3) = eFrame_Row(7,atom1)
597     uy_i(1) = eFrame_Row(2,atom1)
598     uy_i(2) = eFrame_Row(5,atom1)
599     uy_i(3) = eFrame_Row(8,atom1)
600     uz_i(1) = eFrame_Row(3,atom1)
601     uz_i(2) = eFrame_Row(6,atom1)
602     uz_i(3) = eFrame_Row(9,atom1)
603     #else
604     ux_i(1) = eFrame(1,atom1)
605     ux_i(2) = eFrame(4,atom1)
606     ux_i(3) = eFrame(7,atom1)
607     uy_i(1) = eFrame(2,atom1)
608     uy_i(2) = eFrame(5,atom1)
609     uy_i(3) = eFrame(8,atom1)
610     uz_i(1) = eFrame(3,atom1)
611     uz_i(2) = eFrame(6,atom1)
612     uz_i(3) = eFrame(9,atom1)
613     #endif
614     cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
615     cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
616     cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
617 gezelter 939 dudux_i = zero
618     duduy_i = zero
619     duduz_i = zero
620 gezelter 439 endif
621    
622 gezelter 411 if (j_is_Charge) then
623     q_j = ElectrostaticMap(me2)%charge
624     endif
625 gezelter 507
626 gezelter 411 if (j_is_Dipole) then
627     mu_j = ElectrostaticMap(me2)%dipole_moment
628     #ifdef IS_MPI
629 gezelter 439 uz_j(1) = eFrame_Col(3,atom2)
630     uz_j(2) = eFrame_Col(6,atom2)
631     uz_j(3) = eFrame_Col(9,atom2)
632 gezelter 411 #else
633 gezelter 439 uz_j(1) = eFrame(3,atom2)
634     uz_j(2) = eFrame(6,atom2)
635     uz_j(3) = eFrame(9,atom2)
636 gezelter 411 #endif
637 chrisfen 465 ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
638 gezelter 411
639     if (j_is_SplitDipole) then
640     d_j = ElectrostaticMap(me2)%split_dipole_distance
641     endif
642 gezelter 939 duduz_j = zero
643 gezelter 411 endif
644    
645 gezelter 439 if (j_is_Quadrupole) then
646     qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1)
647     qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2)
648     qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3)
649     #ifdef IS_MPI
650     ux_j(1) = eFrame_Col(1,atom2)
651     ux_j(2) = eFrame_Col(4,atom2)
652     ux_j(3) = eFrame_Col(7,atom2)
653     uy_j(1) = eFrame_Col(2,atom2)
654     uy_j(2) = eFrame_Col(5,atom2)
655     uy_j(3) = eFrame_Col(8,atom2)
656     uz_j(1) = eFrame_Col(3,atom2)
657     uz_j(2) = eFrame_Col(6,atom2)
658     uz_j(3) = eFrame_Col(9,atom2)
659     #else
660     ux_j(1) = eFrame(1,atom2)
661     ux_j(2) = eFrame(4,atom2)
662     ux_j(3) = eFrame(7,atom2)
663     uy_j(1) = eFrame(2,atom2)
664     uy_j(2) = eFrame(5,atom2)
665     uy_j(3) = eFrame(8,atom2)
666     uz_j(1) = eFrame(3,atom2)
667     uz_j(2) = eFrame(6,atom2)
668     uz_j(3) = eFrame(9,atom2)
669     #endif
670     cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
671     cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
672     cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
673 gezelter 939 dudux_j = zero
674     duduy_j = zero
675     duduz_j = zero
676 gezelter 439 endif
677 chrisfen 554
678 gezelter 939 epot = zero
679     dudx = zero
680     dudy = zero
681     dudz = zero
682 gezelter 411
683     if (i_is_Charge) then
684    
685     if (j_is_Charge) then
686 chrisfen 739 if (screeningMethod .eq. DAMPED) then
687 chrisfen 959 ! assemble the damping variables
688 chrisfen 941 call lookupUniformSpline1d(f0spline, rij, f0, df0)
689     f1 = -rij * df0 + f0
690 chrisfen 739 endif
691 gezelter 507
692 chrisfen 739 preVal = pre11 * q_i * q_j
693    
694 chrisfen 710 if (summationMethod .eq. SHIFTED_POTENTIAL) then
695 chrisfen 739 vterm = preVal * (riji*f0 - rcuti*f0c)
696 chrisfen 597
697 chrisfen 739 dudr = -sw * preVal * riji * riji * f1
698    
699 chrisfen 710 elseif (summationMethod .eq. SHIFTED_FORCE) then
700 chrisfen 739 vterm = preVal * ( riji*f0 - rcuti*f0c + &
701 chrisfen 716 f1c*rcuti2*(rij-defaultCutoff) )
702    
703 chrisfen 739 dudr = -sw*preVal * (riji*riji*f1 - rcuti2*f1c)
704    
705 chrisfen 695 elseif (summationMethod .eq. REACTION_FIELD) then
706     rfVal = preRF*rij*rij
707     vterm = preVal * ( riji + rfVal )
708 chrisfen 700
709 chrisfen 959 dudr = sw * preVal * ( 2.0_dp*rfVal - riji )*riji
710 chrisfen 739
711 chrisfen 597 else
712 chrisfen 739 vterm = preVal * riji*f0
713 chrisfen 597
714 chrisfen 739 dudr = - sw * preVal * riji*riji*f1
715    
716 chrisfen 597 endif
717    
718 chrisfen 739 vpair = vpair + vterm
719     epot = epot + sw*vterm
720    
721     dudx = dudx + dudr * xhat
722     dudy = dudy + dudr * yhat
723     dudz = dudz + dudr * zhat
724    
725 gezelter 411 endif
726    
727     if (j_is_Dipole) then
728 chrisfen 849 if (screeningMethod .eq. DAMPED) then
729 chrisfen 959 ! assemble the damping variables
730 chrisfen 941 call lookupUniformSpline1d(f0spline, rij, f0, df0)
731     f1 = -rij * df0 + f0
732 chrisfen 959 f3 = -2.0_dp*alpha2*df0*rij*rij*rij
733 chrisfen 849 endif
734 gezelter 411
735 chrisfen 626 pref = pre12 * q_i * mu_j
736 gezelter 411
737 chrisfen 710 if (summationMethod .eq. REACTION_FIELD) then
738 chrisfen 700 ri2 = riji * riji
739     ri3 = ri2 * riji
740 chrisfen 696
741     vterm = - pref * ct_j * ( ri2 - preRF2*rij )
742     vpair = vpair + vterm
743     epot = epot + sw*vterm
744    
745 chrisfen 959 dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0_dp*ct_j*xhat) - &
746 chrisfen 696 preRF2*uz_j(1) )
747 chrisfen 959 dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0_dp*ct_j*yhat) - &
748 chrisfen 696 preRF2*uz_j(2) )
749 chrisfen 959 dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0_dp*ct_j*zhat) - &
750 chrisfen 696 preRF2*uz_j(3) )
751     duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
752     duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
753     duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
754    
755 chrisfen 597 else
756     if (j_is_SplitDipole) then
757 chrisfen 959 BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
758     ri = 1.0_dp / BigR
759 chrisfen 597 scale = rij * ri
760     else
761     ri = riji
762 chrisfen 959 scale = 1.0_dp
763 chrisfen 597 endif
764    
765     ri2 = ri * ri
766     ri3 = ri2 * ri
767     sc2 = scale * scale
768 chrisfen 626
769 chrisfen 849 pot_term = ri2 * scale * f1
770 chrisfen 959 vterm = -pref * ct_j * pot_term
771 chrisfen 626 vpair = vpair + vterm
772     epot = epot + sw*vterm
773 chrisfen 597
774 chrisfen 959 prei3 = sw*pref*ri3
775     ri5damp = 3.0_dp*f1 + f3
776 chrisfen 597
777 chrisfen 959 dudx = dudx - prei3 * ( uz_j(1)*f1 - ct_j*xhat*sc2*ri5damp )
778     dudy = dudy - prei3 * ( uz_j(2)*f1 - ct_j*yhat*sc2*ri5damp )
779     dudz = dudz - prei3 * ( uz_j(3)*f1 - ct_j*zhat*sc2*ri5damp )
780 chrisfen 849
781     duduz_j(1) = duduz_j(1) - sw*pref * pot_term * xhat
782     duduz_j(2) = duduz_j(2) - sw*pref * pot_term * yhat
783     duduz_j(3) = duduz_j(3) - sw*pref * pot_term * zhat
784 gezelter 411
785 chrisfen 597 endif
786 gezelter 411 endif
787 gezelter 421
788 gezelter 439 if (j_is_Quadrupole) then
789 chrisfen 849 if (screeningMethod .eq. DAMPED) then
790 chrisfen 959 ! assemble the damping variables
791 chrisfen 941 call lookupUniformSpline1d(f0spline, rij, f0, df0)
792     f1 = -rij * df0 + f0
793 chrisfen 959 f2 = -2.0_dp*alpha2*df0
794     f3 = f2*r2*rij
795     f4 = 0.4_dp*alpha2*f3*r2
796 chrisfen 849 endif
797 chrisfen 959 ri5damp = f1 + f3*one_third
798 chrisfen 986 ri7damp = ri5damp + f4*one_third
799 chrisfen 849
800 gezelter 439 ri2 = riji * riji
801     ri3 = ri2 * riji
802     cx2 = cx_j * cx_j
803     cy2 = cy_j * cy_j
804     cz2 = cz_j * cz_j
805    
806 chrisfen 959 pref = pre14 * q_i * one_third
807    
808     pot_term = ri3*( qxx_j*(3.0_dp*cx2*ri5damp - f1) + &
809     qyy_j*(3.0_dp*cy2*ri5damp - f1) + &
810     qzz_j*(3.0_dp*cz2*ri5damp - f1) )
811     vterm = pref * pot_term
812 chrisfen 740 vpair = vpair + vterm
813     epot = epot + sw*vterm
814 chrisfen 959
815     ! precompute variables for convenience (and obfuscation unfortunatly)
816     prei3 = 3.0_dp*sw*pref*ri3
817     prei4 = prei3*riji
818     xhatdot2 = xhat*2.0_dp * ri5damp
819     yhatdot2 = yhat*2.0_dp * ri5damp
820     zhatdot2 = zhat*2.0_dp * ri5damp
821     xhatdot5 = xhat*5.0_dp * ri7damp
822     yhatdot5 = yhat*5.0_dp * ri7damp
823     zhatdot5 = zhat*5.0_dp * ri7damp
824    
825     dudx = dudx - prei4 * ( &
826     qxx_j*(cx2*xhatdot5 - (2.0_dp*cx_j*ux_j(1) + xhat)*ri5damp) + &
827     qyy_j*(cy2*xhatdot5 - (2.0_dp*cy_j*uy_j(1) + xhat)*ri5damp) + &
828     qzz_j*(cz2*xhatdot5 - (2.0_dp*cz_j*uz_j(1) + xhat)*ri5damp) )
829     dudy = dudy - prei4 * ( &
830     qxx_j*(cx2*yhatdot5 - (2.0_dp*cx_j*ux_j(2) + yhat)*ri5damp) + &
831     qyy_j*(cy2*yhatdot5 - (2.0_dp*cy_j*uy_j(2) + yhat)*ri5damp) + &
832     qzz_j*(cz2*yhatdot5 - (2.0_dp*cz_j*uz_j(2) + yhat)*ri5damp) )
833     dudz = dudz - prei4 * ( &
834     qxx_j*(cx2*zhatdot5 - (2.0_dp*cx_j*ux_j(3) + zhat)*ri5damp) + &
835     qyy_j*(cy2*zhatdot5 - (2.0_dp*cy_j*uy_j(3) + zhat)*ri5damp) + &
836     qzz_j*(cz2*zhatdot5 - (2.0_dp*cz_j*uz_j(3) + zhat)*ri5damp) )
837 chrisfen 597
838 chrisfen 959 dudux_j(1) = dudux_j(1) + prei3*(qxx_j*cx_j*xhatdot2)
839     dudux_j(2) = dudux_j(2) + prei3*(qxx_j*cx_j*yhatdot2)
840     dudux_j(3) = dudux_j(3) + prei3*(qxx_j*cx_j*zhatdot2)
841 chrisfen 740
842 chrisfen 959 duduy_j(1) = duduy_j(1) + prei3*(qyy_j*cy_j*xhatdot2)
843     duduy_j(2) = duduy_j(2) + prei3*(qyy_j*cy_j*yhatdot2)
844     duduy_j(3) = duduy_j(3) + prei3*(qyy_j*cy_j*zhatdot2)
845 chrisfen 740
846 chrisfen 959 duduz_j(1) = duduz_j(1) + prei3*(qzz_j*cz_j*xhatdot2)
847     duduz_j(2) = duduz_j(2) + prei3*(qzz_j*cz_j*yhatdot2)
848     duduz_j(3) = duduz_j(3) + prei3*(qzz_j*cz_j*zhatdot2)
849    
850 chrisfen 849
851 gezelter 439 endif
852 gezelter 411 endif
853 chrisfen 740
854 gezelter 411 if (i_is_Dipole) then
855 gezelter 507
856 gezelter 411 if (j_is_Charge) then
857 chrisfen 849 if (screeningMethod .eq. DAMPED) then
858 chrisfen 959 ! assemble the damping variables
859 chrisfen 941 call lookupUniformSpline1d(f0spline, rij, f0, df0)
860     f1 = -rij * df0 + f0
861 chrisfen 959 f3 = -2.0_dp*alpha2*df0*r2*rij
862 chrisfen 849 endif
863 chrisfen 626
864 chrisfen 740 pref = pre12 * q_j * mu_i
865    
866 chrisfen 959 if (summationMethod .eq. REACTION_FIELD) then
867 gezelter 507
868 chrisfen 719 ri2 = riji * riji
869     ri3 = ri2 * riji
870    
871 chrisfen 700 vterm = pref * ct_i * ( ri2 - preRF2*rij )
872 chrisfen 696 vpair = vpair + vterm
873     epot = epot + sw*vterm
874    
875 chrisfen 959 dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0_dp*ct_i*xhat) - &
876 chrisfen 700 preRF2*uz_i(1) )
877 chrisfen 959 dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0_dp*ct_i*yhat) - &
878 chrisfen 700 preRF2*uz_i(2) )
879 chrisfen 959 dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0_dp*ct_i*zhat) - &
880 chrisfen 700 preRF2*uz_i(3) )
881 chrisfen 696
882 chrisfen 700 duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
883     duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
884     duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
885 chrisfen 696
886 chrisfen 597 else
887     if (i_is_SplitDipole) then
888 chrisfen 959 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
889     ri = 1.0_dp / BigR
890 chrisfen 597 scale = rij * ri
891     else
892 gezelter 421 ri = riji
893 chrisfen 959 scale = 1.0_dp
894 gezelter 421 endif
895 chrisfen 597
896     ri2 = ri * ri
897     ri3 = ri2 * ri
898     sc2 = scale * scale
899 chrisfen 626
900 chrisfen 849 pot_term = ri2 * f1 * scale
901     vterm = pref * ct_i * pot_term
902 chrisfen 626 vpair = vpair + vterm
903     epot = epot + sw*vterm
904 chrisfen 597
905 chrisfen 959 prei3 = sw*pref*ri3
906     ri5damp = 3.0_dp*f1 + f3
907 chrisfen 597
908 chrisfen 959 dudx = dudx + prei3 * ( uz_i(1)*f1 - ct_i*xhat*sc2*ri5damp )
909     dudy = dudy + prei3 * ( uz_i(2)*f1 - ct_i*yhat*sc2*ri5damp )
910     dudz = dudz + prei3 * ( uz_i(3)*f1 - ct_i*zhat*sc2*ri5damp )
911    
912 chrisfen 849 duduz_i(1) = duduz_i(1) + sw*pref * pot_term * xhat
913     duduz_i(2) = duduz_i(2) + sw*pref * pot_term * yhat
914     duduz_i(3) = duduz_i(3) + sw*pref * pot_term * zhat
915 chrisfen 959
916 gezelter 421 endif
917 chrisfen 597 endif
918 chrisfen 626
919 chrisfen 597 if (j_is_Dipole) then
920 chrisfen 986 if (screeningMethod .eq. DAMPED) then
921     ! assemble the damping variables
922     call lookupUniformSpline1d(f0spline, rij, f0, df0)
923     f1 = -rij * df0 + f0
924     f2 = -2.0_dp*alpha2*df0
925     f3 = f2*r2*rij
926     f4 = 0.4_dp*alpha2*f3*r2
927     endif
928 chrisfen 849
929 chrisfen 719 ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
930    
931     ri2 = riji * riji
932     ri3 = ri2 * riji
933     ri4 = ri2 * ri2
934    
935     pref = pre22 * mu_i * mu_j
936 gezelter 421
937 chrisfen 710 if (summationMethod .eq. REACTION_FIELD) then
938 chrisfen 959 vterm = pref*( ri3*(ct_ij - 3.0_dp * ct_i * ct_j) - &
939 chrisfen 695 preRF2*ct_ij )
940     vpair = vpair + vterm
941     epot = epot + sw*vterm
942    
943 chrisfen 959 a1 = 5.0_dp * ct_i * ct_j - ct_ij
944 chrisfen 695
945 chrisfen 959 dudx = dudx + sw*pref*3.0_dp*ri4 &
946 chrisfen 695 * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
947 chrisfen 959 dudy = dudy + sw*pref*3.0_dp*ri4 &
948 chrisfen 695 * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
949 chrisfen 959 dudz = dudz + sw*pref*3.0_dp*ri4 &
950 chrisfen 695 * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
951    
952 chrisfen 959 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0_dp*ct_j*xhat) &
953 chrisfen 695 - preRF2*uz_j(1))
954 chrisfen 959 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0_dp*ct_j*yhat) &
955 chrisfen 695 - preRF2*uz_j(2))
956 chrisfen 959 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0_dp*ct_j*zhat) &
957 chrisfen 695 - preRF2*uz_j(3))
958 chrisfen 959 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0_dp*ct_i*xhat) &
959 chrisfen 695 - preRF2*uz_i(1))
960 chrisfen 959 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0_dp*ct_i*yhat) &
961 chrisfen 695 - preRF2*uz_i(2))
962 chrisfen 959 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0_dp*ct_i*zhat) &
963 chrisfen 695 - preRF2*uz_i(3))
964    
965 chrisfen 597 else
966     if (i_is_SplitDipole) then
967     if (j_is_SplitDipole) then
968 chrisfen 959 BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
969 chrisfen 597 else
970 chrisfen 959 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
971 chrisfen 597 endif
972 chrisfen 959 ri = 1.0_dp / BigR
973 chrisfen 597 scale = rij * ri
974     else
975     if (j_is_SplitDipole) then
976 chrisfen 959 BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
977     ri = 1.0_dp / BigR
978 chrisfen 597 scale = rij * ri
979     else
980     ri = riji
981 chrisfen 959 scale = 1.0_dp
982 chrisfen 597 endif
983     endif
984 chrisfen 719
985 chrisfen 959 ! precompute variables for convenience (and obfuscation
986     ! unfortunatly)
987 chrisfen 986 sc2 = scale * scale
988     ri5damp = f1 + f3*one_third
989     ri7damp = 5.0_dp*(ri5damp + f4*one_third)
990 chrisfen 959 prei3 = sw*pref*ri3
991     prei4 = 3.0_dp*sw*pref*ri4*scale
992 chrisfen 986 cti3 = 3.0_dp*ct_i*sc2*ri5damp
993     ctj3 = 3.0_dp*ct_j*sc2*ri5damp
994     ctidotj = ct_i * ct_j * sc2
995 chrisfen 959
996 chrisfen 986 ! calculate the potential
997     pot_term = (ct_ij*f1 - 3.0_dp*ctidotj*ri5damp)
998     vterm = pref * ri3 * pot_term
999     vpair = vpair + vterm
1000     epot = epot + sw*vterm
1001 chrisfen 959
1002 chrisfen 986 ! calculate derivatives for the forces and torques
1003     dudx = dudx + prei4 * ( ctidotj*xhat*ri7damp - &
1004     (ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat)*ri5damp )
1005     dudy = dudy + prei4 * ( ctidotj*yhat*ri7damp - &
1006     (ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat)*ri5damp )
1007     dudz = dudz + prei4 * ( ctidotj*zhat*ri7damp - &
1008     (ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat)*ri5damp )
1009    
1010     duduz_i(1) = duduz_i(1) + prei3 * ( uz_j(1)*f1 - ctj3*xhat )
1011     duduz_i(2) = duduz_i(2) + prei3 * ( uz_j(2)*f1 - ctj3*yhat )
1012     duduz_i(3) = duduz_i(3) + prei3 * ( uz_j(3)*f1 - ctj3*zhat )
1013 chrisfen 597
1014 chrisfen 986 duduz_j(1) = duduz_j(1) + prei3 * ( uz_i(1)*f1 - cti3*xhat )
1015     duduz_j(2) = duduz_j(2) + prei3 * ( uz_i(2)*f1 - cti3*yhat )
1016     duduz_j(3) = duduz_j(3) + prei3 * ( uz_i(3)*f1 - cti3*zhat )
1017 chrisfen 849
1018 chrisfen 597 endif
1019 gezelter 411 endif
1020     endif
1021 gezelter 439
1022     if (i_is_Quadrupole) then
1023     if (j_is_Charge) then
1024 chrisfen 849 if (screeningMethod .eq. DAMPED) then
1025 chrisfen 959 ! assemble the damping variables
1026 chrisfen 941 call lookupUniformSpline1d(f0spline, rij, f0, df0)
1027     f1 = -rij * df0 + f0
1028 chrisfen 959 f2 = -2.0_dp*alpha2*df0
1029     f3 = f2*r2*rij
1030     f4 = 0.4_dp*alpha2*f3*r2
1031 chrisfen 849 endif
1032 chrisfen 959 ri5damp = f1 + f3*one_third
1033 chrisfen 986 ri7damp = ri5damp + f4*one_third
1034 chrisfen 849
1035 gezelter 439 ri2 = riji * riji
1036     ri3 = ri2 * riji
1037 gezelter 440 ri4 = ri2 * ri2
1038 gezelter 439 cx2 = cx_i * cx_i
1039     cy2 = cy_i * cy_i
1040     cz2 = cz_i * cz_i
1041 gezelter 507
1042 chrisfen 959 pref = pre14 * q_j * one_third
1043    
1044     pot_term = ri3 * ( qxx_i * (3.0_dp*cx2*ri5damp - f1) + &
1045     qyy_i * (3.0_dp*cy2*ri5damp - f1) + &
1046     qzz_i * (3.0_dp*cz2*ri5damp - f1) )
1047    
1048     vterm = pref * pot_term
1049 chrisfen 740 vpair = vpair + vterm
1050     epot = epot + sw*vterm
1051 chrisfen 959
1052     ! precompute variables for convenience (and obfuscation unfortunatly)
1053     prei3 = 3.0_dp*sw*pref*ri3
1054     prei4 = prei3*riji
1055     xhatdot2 = xhat*2.0_dp * ri5damp
1056     yhatdot2 = yhat*2.0_dp * ri5damp
1057     zhatdot2 = zhat*2.0_dp * ri5damp
1058     xhatdot5 = xhat*5.0_dp * ri7damp
1059     yhatdot5 = yhat*5.0_dp * ri7damp
1060     zhatdot5 = zhat*5.0_dp * ri7damp
1061    
1062     dudx = dudx - prei4 * ( &
1063     qxx_i*(cx2*xhatdot5 - (2.0_dp*cx_i*ux_i(1) + xhat)*ri5damp) + &
1064     qyy_i*(cy2*xhatdot5 - (2.0_dp*cy_i*uy_i(1) + xhat)*ri5damp) + &
1065     qzz_i*(cz2*xhatdot5 - (2.0_dp*cz_i*uz_i(1) + xhat)*ri5damp) )
1066     dudy = dudy - prei4 * ( &
1067     qxx_i*(cx2*yhatdot5 - (2.0_dp*cx_i*ux_i(2) + yhat)*ri5damp) + &
1068     qyy_i*(cy2*yhatdot5 - (2.0_dp*cy_i*uy_i(2) + yhat)*ri5damp) + &
1069     qzz_i*(cz2*yhatdot5 - (2.0_dp*cz_i*uz_i(2) + yhat)*ri5damp) )
1070     dudz = dudz - prei4 * ( &
1071     qxx_i*(cx2*zhatdot5 - (2.0_dp*cx_i*ux_i(3) + zhat)*ri5damp) + &
1072     qyy_i*(cy2*zhatdot5 - (2.0_dp*cy_i*uy_i(3) + zhat)*ri5damp) + &
1073     qzz_i*(cz2*zhatdot5 - (2.0_dp*cz_i*uz_i(3) + zhat)*ri5damp) )
1074 chrisfen 740
1075 chrisfen 959 dudux_i(1) = dudux_i(1) + prei3*(qxx_i*cx_i*xhatdot2)
1076     dudux_i(2) = dudux_i(2) + prei3*(qxx_i*cx_i*yhatdot2)
1077     dudux_i(3) = dudux_i(3) + prei3*(qxx_i*cx_i*zhatdot2)
1078 chrisfen 740
1079 chrisfen 959 duduy_i(1) = duduy_i(1) + prei3*(qyy_i*cy_i*xhatdot2)
1080     duduy_i(2) = duduy_i(2) + prei3*(qyy_i*cy_i*yhatdot2)
1081     duduy_i(3) = duduy_i(3) + prei3*(qyy_i*cy_i*zhatdot2)
1082 chrisfen 740
1083 chrisfen 959 duduz_i(1) = duduz_i(1) + prei3*(qzz_i*cz_i*xhatdot2)
1084     duduz_i(2) = duduz_i(2) + prei3*(qzz_i*cz_i*yhatdot2)
1085     duduz_i(3) = duduz_i(3) + prei3*(qzz_i*cz_i*zhatdot2)
1086 gezelter 439 endif
1087     endif
1088 gezelter 507
1089    
1090 gezelter 411 if (do_pot) then
1091     #ifdef IS_MPI
1092 chrisfen 959 pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5_dp*epot
1093     pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5_dp*epot
1094 gezelter 411 #else
1095     pot = pot + epot
1096     #endif
1097     endif
1098 gezelter 507
1099 gezelter 411 #ifdef IS_MPI
1100     f_Row(1,atom1) = f_Row(1,atom1) + dudx
1101     f_Row(2,atom1) = f_Row(2,atom1) + dudy
1102     f_Row(3,atom1) = f_Row(3,atom1) + dudz
1103 gezelter 507
1104 gezelter 411 f_Col(1,atom2) = f_Col(1,atom2) - dudx
1105     f_Col(2,atom2) = f_Col(2,atom2) - dudy
1106     f_Col(3,atom2) = f_Col(3,atom2) - dudz
1107 gezelter 507
1108 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1109 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1110     t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1111     t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1112 gezelter 411 endif
1113 gezelter 439 if (i_is_Quadrupole) then
1114     t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1115     t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1116     t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1117 gezelter 411
1118 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1119     t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1120     t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1121     endif
1122    
1123 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1124 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1125     t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1126     t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1127 gezelter 411 endif
1128 gezelter 439 if (j_is_Quadrupole) then
1129     t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1130     t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1131     t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1132 gezelter 411
1133 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1134     t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1135     t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1136     endif
1137    
1138 gezelter 411 #else
1139     f(1,atom1) = f(1,atom1) + dudx
1140     f(2,atom1) = f(2,atom1) + dudy
1141     f(3,atom1) = f(3,atom1) + dudz
1142 gezelter 507
1143 gezelter 411 f(1,atom2) = f(1,atom2) - dudx
1144     f(2,atom2) = f(2,atom2) - dudy
1145     f(3,atom2) = f(3,atom2) - dudz
1146 gezelter 507
1147 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1148 gezelter 439 t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1149     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1150     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1151 gezelter 411 endif
1152 gezelter 439 if (i_is_Quadrupole) then
1153     t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1154     t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1155     t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1156    
1157     t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1158     t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1159     t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1160     endif
1161    
1162 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1163 gezelter 439 t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1164     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1165     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1166 gezelter 411 endif
1167 gezelter 439 if (j_is_Quadrupole) then
1168     t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1169     t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1170     t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1171    
1172     t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1173     t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1174     t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1175     endif
1176    
1177 gezelter 411 #endif
1178 gezelter 507
1179 gezelter 411 #ifdef IS_MPI
1180     id1 = AtomRowToGlobal(atom1)
1181     id2 = AtomColToGlobal(atom2)
1182     #else
1183     id1 = atom1
1184     id2 = atom2
1185     #endif
1186    
1187     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1188 gezelter 507
1189 gezelter 411 fpair(1) = fpair(1) + dudx
1190     fpair(2) = fpair(2) + dudy
1191     fpair(3) = fpair(3) + dudz
1192    
1193     endif
1194    
1195     return
1196     end subroutine doElectrostaticPair
1197 chuckv 492
1198     subroutine destroyElectrostaticTypes()
1199    
1200 gezelter 507 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1201    
1202 chuckv 492 end subroutine destroyElectrostaticTypes
1203    
1204 chrisfen 703 subroutine self_self(atom1, eFrame, mypot, t, do_pot)
1205 chrisfen 695 logical, intent(in) :: do_pot
1206 chrisfen 682 integer, intent(in) :: atom1
1207 chrisfen 695 integer :: atid1
1208 chrisfen 682 real(kind=dp), dimension(9,nLocal) :: eFrame
1209 chrisfen 695 real(kind=dp), dimension(3,nLocal) :: t
1210 chrisfen 703 real(kind=dp) :: mu1, c1
1211     real(kind=dp) :: preVal, epot, mypot
1212 chrisfen 695 real(kind=dp) :: eix, eiy, eiz
1213 chrisfen 682
1214 chrisfen 695 ! this is a local only array, so we use the local atom type id's:
1215     atid1 = atid(atom1)
1216 chrisfen 703
1217     if (.not.summationMethodChecked) then
1218     call checkSummationMethod()
1219     endif
1220 chrisfen 695
1221 chrisfen 703 if (summationMethod .eq. REACTION_FIELD) then
1222     if (ElectrostaticMap(atid1)%is_Dipole) then
1223     mu1 = getDipoleMoment(atid1)
1224    
1225     preVal = pre22 * preRF2 * mu1*mu1
1226 chrisfen 959 mypot = mypot - 0.5_dp*preVal
1227 chrisfen 703
1228     ! The self-correction term adds into the reaction field vector
1229    
1230     eix = preVal * eFrame(3,atom1)
1231     eiy = preVal * eFrame(6,atom1)
1232     eiz = preVal * eFrame(9,atom1)
1233    
1234     ! once again, this is self-self, so only the local arrays are needed
1235     ! even for MPI jobs:
1236    
1237     t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1238     eFrame(9,atom1)*eiy
1239     t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1240     eFrame(3,atom1)*eiz
1241     t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1242     eFrame(6,atom1)*eix
1243    
1244     endif
1245    
1246 chrisfen 743 elseif ( (summationMethod .eq. SHIFTED_FORCE) .or. &
1247     (summationMethod .eq. SHIFTED_POTENTIAL) ) then
1248 chrisfen 717 if (ElectrostaticMap(atid1)%is_Charge) then
1249     c1 = getCharge(atid1)
1250    
1251     if (screeningMethod .eq. DAMPED) then
1252 chrisfen 959 mypot = mypot - (f0c * rcuti * 0.5_dp + &
1253 chrisfen 717 dampingAlpha*invRootPi) * c1 * c1
1254    
1255     else
1256 chrisfen 959 mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1257 chrisfen 717
1258     endif
1259     endif
1260     endif
1261 chrisfen 695
1262 chrisfen 682 return
1263 chrisfen 703 end subroutine self_self
1264 chrisfen 682
1265 chrisfen 703 subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, myPot, &
1266 chrisfen 700 f, t, do_pot)
1267     logical, intent(in) :: do_pot
1268     integer, intent(in) :: atom1
1269     integer, intent(in) :: atom2
1270     logical :: i_is_Charge, j_is_Charge
1271     logical :: i_is_Dipole, j_is_Dipole
1272     integer :: atid1
1273     integer :: atid2
1274     real(kind=dp), intent(in) :: rij
1275     real(kind=dp), intent(in) :: sw
1276     real(kind=dp), intent(in), dimension(3) :: d
1277     real(kind=dp), intent(inout) :: vpair
1278     real(kind=dp), dimension(9,nLocal) :: eFrame
1279     real(kind=dp), dimension(3,nLocal) :: f
1280     real(kind=dp), dimension(3,nLocal) :: t
1281     real (kind = dp), dimension(3) :: duduz_i
1282     real (kind = dp), dimension(3) :: duduz_j
1283     real (kind = dp), dimension(3) :: uz_i
1284     real (kind = dp), dimension(3) :: uz_j
1285     real(kind=dp) :: q_i, q_j, mu_i, mu_j
1286     real(kind=dp) :: xhat, yhat, zhat
1287     real(kind=dp) :: ct_i, ct_j
1288     real(kind=dp) :: ri2, ri3, riji, vterm
1289 chrisfen 703 real(kind=dp) :: pref, preVal, rfVal, myPot
1290 chrisfen 700 real(kind=dp) :: dudx, dudy, dudz, dudr
1291    
1292 chrisfen 703 if (.not.summationMethodChecked) then
1293     call checkSummationMethod()
1294 chrisfen 700 endif
1295    
1296 gezelter 939 dudx = zero
1297     dudy = zero
1298     dudz = zero
1299 chrisfen 700
1300 chrisfen 959 riji = 1.0_dp/rij
1301 chrisfen 700
1302     xhat = d(1) * riji
1303     yhat = d(2) * riji
1304     zhat = d(3) * riji
1305    
1306     ! this is a local only array, so we use the local atom type id's:
1307     atid1 = atid(atom1)
1308     atid2 = atid(atom2)
1309     i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1310     j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1311     i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1312     j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1313    
1314     if (i_is_Charge.and.j_is_Charge) then
1315     q_i = ElectrostaticMap(atid1)%charge
1316     q_j = ElectrostaticMap(atid2)%charge
1317    
1318     preVal = pre11 * q_i * q_j
1319     rfVal = preRF*rij*rij
1320     vterm = preVal * rfVal
1321    
1322 chrisfen 703 myPot = myPot + sw*vterm
1323    
1324 chrisfen 959 dudr = sw*preVal * 2.0_dp*rfVal*riji
1325 chrisfen 703
1326 chrisfen 700 dudx = dudx + dudr * xhat
1327     dudy = dudy + dudr * yhat
1328     dudz = dudz + dudr * zhat
1329 chrisfen 703
1330 chrisfen 700 elseif (i_is_Charge.and.j_is_Dipole) then
1331     q_i = ElectrostaticMap(atid1)%charge
1332     mu_j = ElectrostaticMap(atid2)%dipole_moment
1333     uz_j(1) = eFrame(3,atom2)
1334     uz_j(2) = eFrame(6,atom2)
1335     uz_j(3) = eFrame(9,atom2)
1336     ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1337 chrisfen 703
1338 chrisfen 700 ri2 = riji * riji
1339     ri3 = ri2 * riji
1340    
1341     pref = pre12 * q_i * mu_j
1342     vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1343 chrisfen 703 myPot = myPot + sw*vterm
1344    
1345 chrisfen 959 dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0_dp*ct_j*xhat) &
1346 chrisfen 703 - preRF2*uz_j(1) )
1347 chrisfen 959 dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0_dp*ct_j*yhat) &
1348 chrisfen 703 - preRF2*uz_j(2) )
1349 chrisfen 959 dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0_dp*ct_j*zhat) &
1350 chrisfen 703 - preRF2*uz_j(3) )
1351    
1352 chrisfen 700 duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1353     duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1354     duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1355 chrisfen 703
1356 chrisfen 700 elseif (i_is_Dipole.and.j_is_Charge) then
1357     mu_i = ElectrostaticMap(atid1)%dipole_moment
1358     q_j = ElectrostaticMap(atid2)%charge
1359     uz_i(1) = eFrame(3,atom1)
1360     uz_i(2) = eFrame(6,atom1)
1361     uz_i(3) = eFrame(9,atom1)
1362     ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1363 chrisfen 703
1364 chrisfen 700 ri2 = riji * riji
1365     ri3 = ri2 * riji
1366    
1367     pref = pre12 * q_j * mu_i
1368     vterm = pref * ct_i * ( ri2 - preRF2*rij )
1369 chrisfen 703 myPot = myPot + sw*vterm
1370 chrisfen 700
1371 chrisfen 959 dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0_dp*ct_i*xhat) &
1372 chrisfen 703 - preRF2*uz_i(1) )
1373 chrisfen 959 dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0_dp*ct_i*yhat) &
1374 chrisfen 703 - preRF2*uz_i(2) )
1375 chrisfen 959 dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0_dp*ct_i*zhat) &
1376 chrisfen 703 - preRF2*uz_i(3) )
1377 chrisfen 700
1378     duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1379     duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1380     duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1381    
1382     endif
1383 chrisfen 703
1384    
1385     ! accumulate the forces and torques resulting from the self term
1386 chrisfen 700 f(1,atom1) = f(1,atom1) + dudx
1387     f(2,atom1) = f(2,atom1) + dudy
1388     f(3,atom1) = f(3,atom1) + dudz
1389    
1390     f(1,atom2) = f(1,atom2) - dudx
1391     f(2,atom2) = f(2,atom2) - dudy
1392     f(3,atom2) = f(3,atom2) - dudz
1393    
1394     if (i_is_Dipole) then
1395     t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1396     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1397     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1398     elseif (j_is_Dipole) then
1399     t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1400     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1401     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1402     endif
1403    
1404     return
1405     end subroutine rf_self_excludes
1406    
1407 gezelter 411 end module electrostatic_module