ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 2722
Committed: Thu Apr 20 18:24:24 2006 UTC (18 years, 4 months ago) by gezelter
File size: 52454 byte(s)
Log Message: 
Complete rewrite of spline code and everything that uses it.

File Contents

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