ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 939
Committed: Thu Apr 20 18:24:24 2006 UTC (19 years, 6 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 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 gezelter 938 real(kind=dp), parameter :: pre11 = 332.0637778d0
69 gezelter 434 !! Charge-Dipole, assuming charges are measured in electrons, and
70     !! dipoles are measured in debyes
71 gezelter 938 real(kind=dp), parameter :: pre12 = 69.13373d0
72 gezelter 434 !! Dipole-Dipole, assuming dipoles are measured in debyes
73 gezelter 938 real(kind=dp), parameter :: pre22 = 14.39325d0
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 gezelter 938 real(kind=dp), parameter :: pre14 = 69.13373d0
78 gezelter 411
79 gezelter 939 real(kind=dp), parameter :: zero = 0.0d0
80    
81 chrisfen 712 !! variables to handle different summation methods for long-range
82     !! electrostatics:
83 gezelter 602 integer, save :: summationMethod = NONE
84 chrisfen 710 integer, save :: screeningMethod = UNDAMPED
85 chrisfen 603 logical, save :: summationMethodChecked = .false.
86 gezelter 602 real(kind=DP), save :: defaultCutoff = 0.0_DP
87 chrisfen 682 real(kind=DP), save :: defaultCutoff2 = 0.0_DP
88 gezelter 602 logical, save :: haveDefaultCutoff = .false.
89     real(kind=DP), save :: dampingAlpha = 0.0_DP
90 chrisfen 716 real(kind=DP), save :: alpha2 = 0.0_DP
91 gezelter 602 logical, save :: haveDampingAlpha = .false.
92 chrisfen 682 real(kind=DP), save :: dielectric = 1.0_DP
93 gezelter 602 logical, save :: haveDielectric = .false.
94     real(kind=DP), save :: constEXP = 0.0_DP
95 chrisfen 682 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 695 real(kind=dp), save :: preRF2 = 0.0_DP
106 chrisfen 716 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 849 real(kind=dp), save :: f3 = 0.0_DP
110     real(kind=dp), save :: f4 = 0.0_DP
111 chrisfen 716 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 849 real(kind=dp), save :: f3c = 0.0_DP
115     real(kind=dp), save :: f4c = 0.0_DP
116 chrisfen 716
117 gezelter 809 #if defined(__IFC) || defined(__PGI)
118 chuckv 632 ! error function for ifc version > 7.
119 chuckv 631 double precision, external :: derfc
120 chuckv 632 #endif
121    
122 chrisfen 853 public :: setElectrostaticSummationMethod
123 chrisfen 710 public :: setScreeningMethod
124 gezelter 602 public :: setElectrostaticCutoffRadius
125 chrisfen 710 public :: setDampingAlpha
126 gezelter 602 public :: setReactionFieldDielectric
127 chrisfen 937 public :: buildElectroSplines
128 gezelter 411 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 492 public :: destroyElectrostaticTypes
137 chrisfen 703 public :: self_self
138 chrisfen 700 public :: rf_self_excludes
139 gezelter 411
140 chrisfen 937
141 gezelter 411 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 532 logical :: is_Tap = .false.
148 gezelter 411 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 938 logical, save :: hasElectrostaticMap
157    
158 gezelter 411 contains
159    
160 chrisfen 853 subroutine setElectrostaticSummationMethod(the_ESM)
161 gezelter 602 integer, intent(in) :: the_ESM
162    
163     if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
164 chrisfen 853 call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
165 gezelter 602 endif
166    
167 chrisfen 610 summationMethod = the_ESM
168 chrisfen 626
169 chrisfen 853 end subroutine setElectrostaticSummationMethod
170 gezelter 602
171 chrisfen 710 subroutine setScreeningMethod(the_SM)
172     integer, intent(in) :: the_SM
173     screeningMethod = the_SM
174     end subroutine setScreeningMethod
175    
176 chrisfen 682 subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
177 gezelter 602 real(kind=dp), intent(in) :: thisRcut
178 chrisfen 682 real(kind=dp), intent(in) :: thisRsw
179 gezelter 602 defaultCutoff = thisRcut
180 chrisfen 849 defaultCutoff2 = defaultCutoff*defaultCutoff
181 chrisfen 682 rrf = defaultCutoff
182     rt = thisRsw
183 gezelter 602 haveDefaultCutoff = .true.
184     end subroutine setElectrostaticCutoffRadius
185    
186 chrisfen 710 subroutine setDampingAlpha(thisAlpha)
187 gezelter 602 real(kind=dp), intent(in) :: thisAlpha
188     dampingAlpha = thisAlpha
189 chrisfen 716 alpha2 = dampingAlpha*dampingAlpha
190 gezelter 602 haveDampingAlpha = .true.
191 chrisfen 710 end subroutine setDampingAlpha
192 gezelter 602
193     subroutine setReactionFieldDielectric(thisDielectric)
194     real(kind=dp), intent(in) :: thisDielectric
195     dielectric = thisDielectric
196     haveDielectric = .true.
197     end subroutine setReactionFieldDielectric
198    
199 chrisfen 937 subroutine buildElectroSplines()
200     end subroutine buildElectroSplines
201    
202 gezelter 411 subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
203 chrisfen 532 is_SplitDipole, is_Quadrupole, is_Tap, status)
204 gezelter 507
205 gezelter 411 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 532 logical, intent(in) :: is_Tap
211 gezelter 411 integer, intent(out) :: status
212     integer :: nAtypes, myATID, i, j
213    
214     status = 0
215     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
216 gezelter 507
217 gezelter 411 !! 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 507
222 gezelter 411 nAtypes = getSize(atypes)
223 gezelter 507
224 gezelter 411 if (nAtypes == 0) then
225     status = -1
226     return
227     end if
228 gezelter 507
229 gezelter 938 allocate(ElectrostaticMap(nAtypes))
230 gezelter 507
231 gezelter 411 end if
232    
233     if (myATID .gt. size(ElectrostaticMap)) then
234     status = -1
235     return
236     endif
237 gezelter 507
238 gezelter 411 ! 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 532 ElectrostaticMap(myATID)%is_Tap = is_Tap
246 gezelter 507
247 gezelter 938 hasElectrostaticMap = .true.
248    
249 gezelter 411 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 938 if (.not.hasElectrostaticMap) then
261 gezelter 411 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 507 endif
277 gezelter 411
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 938 if (.not.hasElectrostaticMap) then
291 gezelter 411 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 938 if (.not.hasElectrostaticMap) then
321 gezelter 411 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 938 if (.not.hasElectrostaticMap) then
351 gezelter 411 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 507
368 gezelter 411 do i = 1, 3
369 gezelter 507 ElectrostaticMap(myATID)%quadrupole_moments(i) = &
370     quadrupole_moments(i)
371     enddo
372 gezelter 411
373     end subroutine setQuadrupoleMoments
374    
375 gezelter 507
376 gezelter 411 function getCharge(atid) result (c)
377     integer, intent(in) :: atid
378     integer :: localError
379     real(kind=dp) :: c
380 gezelter 507
381 gezelter 938 if (.not.hasElectrostaticMap) then
382 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
383     return
384     end if
385 gezelter 507
386 gezelter 411 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 507
391 gezelter 411 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 507
399 gezelter 938 if (.not.hasElectrostaticMap) then
400 gezelter 411 call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
401     return
402     end if
403 gezelter 507
404 gezelter 411 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 507
409 gezelter 411 dm = ElectrostaticMap(atid)%dipole_moment
410     end function getDipoleMoment
411    
412 gezelter 602 subroutine checkSummationMethod()
413    
414 chrisfen 607 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 710 if (screeningMethod .eq. DAMPED) then
424 chrisfen 703 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 603
432 chrisfen 849 constEXP = exp(-alpha2*defaultCutoff2)
433 chrisfen 703 invRootPi = 0.56418958354775628695d0
434 chrisfen 716 alphaPi = 2.0d0*dampingAlpha*invRootPi
435     f0c = derfc(dampingAlpha*defaultCutoff)
436     f1c = alphaPi*defaultCutoff*constEXP + f0c
437 chrisfen 849 f2c = alphaPi*2.0d0*alpha2*constEXP
438     f3c = alphaPi*2.0d0*alpha2*constEXP*defaultCutoff2*defaultCutoff
439 gezelter 602 endif
440    
441 chrisfen 603 if (summationMethod .eq. REACTION_FIELD) then
442 chrisfen 703 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 603 endif
450 chrisfen 703
451 chrisfen 603 endif
452    
453     summationMethodChecked = .true.
454 gezelter 602 end subroutine checkSummationMethod
455    
456 chrisfen 712
457 gezelter 762 subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, rcut, sw, &
458 chrisfen 712 vpair, fpair, pot, eFrame, f, t, do_pot)
459 gezelter 507
460 chrisfen 703 logical, intent(in) :: do_pot
461 gezelter 507
462 gezelter 411 integer, intent(in) :: atom1, atom2
463     integer :: localError
464    
465 gezelter 762 real(kind=dp), intent(in) :: rij, r2, sw, rcut
466 gezelter 411 real(kind=dp), intent(in), dimension(3) :: d
467     real(kind=dp), intent(inout) :: vpair
468 chrisfen 703 real(kind=dp), intent(inout), dimension(3) :: fpair
469 gezelter 411
470 chrisfen 626 real( kind = dp ) :: pot
471 gezelter 411 real( kind = dp ), dimension(9,nLocal) :: eFrame
472     real( kind = dp ), dimension(3,nLocal) :: f
473 chrisfen 710 real( kind = dp ), dimension(3,nLocal) :: felec
474 gezelter 411 real( kind = dp ), dimension(3,nLocal) :: t
475 gezelter 507
476 gezelter 439 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 411
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 532 logical :: i_is_Tap, j_is_Tap
484 gezelter 411 integer :: me1, me2, id1, id2
485     real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
486 gezelter 439 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 719 real (kind=dp) :: ct_i, ct_j, ct_ij, a0, a1
492 gezelter 421 real (kind=dp) :: riji, ri, ri2, ri3, ri4
493 chrisfen 597 real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
494 gezelter 421 real (kind=dp) :: xhat, yhat, zhat
495 gezelter 411 real (kind=dp) :: dudx, dudy, dudz
496 chrisfen 626 real (kind=dp) :: scale, sc2, bigR
497 chrisfen 716 real (kind=dp) :: varEXP
498 chrisfen 719 real (kind=dp) :: pot_term
499 chrisfen 695 real (kind=dp) :: preVal, rfVal
500 chrisfen 849 real (kind=dp) :: f13, f134
501 gezelter 411
502 gezelter 602 if (.not.summationMethodChecked) then
503     call checkSummationMethod()
504     endif
505    
506 gezelter 411 #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 644
518 gezelter 421 xhat = d(1) * riji
519     yhat = d(2) * riji
520     zhat = d(3) * riji
521 gezelter 411
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 532 i_is_Tap = ElectrostaticMap(me1)%is_Tap
528 gezelter 411
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 532 j_is_Tap = ElectrostaticMap(me2)%is_Tap
534 gezelter 411
535     if (i_is_Charge) then
536     q_i = ElectrostaticMap(me1)%charge
537     endif
538 gezelter 507
539 gezelter 411 if (i_is_Dipole) then
540     mu_i = ElectrostaticMap(me1)%dipole_moment
541     #ifdef IS_MPI
542 gezelter 439 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 411 #else
546 gezelter 439 uz_i(1) = eFrame(3,atom1)
547     uz_i(2) = eFrame(6,atom1)
548     uz_i(3) = eFrame(9,atom1)
549 gezelter 411 #endif
550 gezelter 439 ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
551 gezelter 411
552     if (i_is_SplitDipole) then
553     d_i = ElectrostaticMap(me1)%split_dipole_distance
554     endif
555 gezelter 939 duduz_i = zero
556 gezelter 411 endif
557    
558 gezelter 439 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 939 dudux_i = zero
587     duduy_i = zero
588     duduz_i = zero
589 gezelter 439 endif
590    
591 gezelter 411 if (j_is_Charge) then
592     q_j = ElectrostaticMap(me2)%charge
593     endif
594 gezelter 507
595 gezelter 411 if (j_is_Dipole) then
596     mu_j = ElectrostaticMap(me2)%dipole_moment
597     #ifdef IS_MPI
598 gezelter 439 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 411 #else
602 gezelter 439 uz_j(1) = eFrame(3,atom2)
603     uz_j(2) = eFrame(6,atom2)
604     uz_j(3) = eFrame(9,atom2)
605 gezelter 411 #endif
606 chrisfen 465 ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
607 gezelter 411
608     if (j_is_SplitDipole) then
609     d_j = ElectrostaticMap(me2)%split_dipole_distance
610     endif
611 gezelter 939 duduz_j = zero
612 gezelter 411 endif
613    
614 gezelter 439 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 939 dudux_j = zero
643     duduy_j = zero
644     duduz_j = zero
645 gezelter 439 endif
646 chrisfen 554
647 gezelter 939 epot = zero
648     dudx = zero
649     dudy = zero
650     dudz = zero
651 gezelter 411
652     if (i_is_Charge) then
653    
654     if (j_is_Charge) then
655 chrisfen 739 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 507
661 chrisfen 739 preVal = pre11 * q_i * q_j
662    
663 chrisfen 710 if (summationMethod .eq. SHIFTED_POTENTIAL) then
664 chrisfen 739 vterm = preVal * (riji*f0 - rcuti*f0c)
665 chrisfen 597
666 chrisfen 739 dudr = -sw * preVal * riji * riji * f1
667    
668 chrisfen 710 elseif (summationMethod .eq. SHIFTED_FORCE) then
669 chrisfen 739 vterm = preVal * ( riji*f0 - rcuti*f0c + &
670 chrisfen 716 f1c*rcuti2*(rij-defaultCutoff) )
671    
672 chrisfen 739 dudr = -sw*preVal * (riji*riji*f1 - rcuti2*f1c)
673    
674 chrisfen 695 elseif (summationMethod .eq. REACTION_FIELD) then
675     rfVal = preRF*rij*rij
676     vterm = preVal * ( riji + rfVal )
677 chrisfen 700
678 chrisfen 695 dudr = sw * preVal * ( 2.0d0*rfVal - riji )*riji
679 chrisfen 739
680 chrisfen 597 else
681 chrisfen 739 vterm = preVal * riji*f0
682 chrisfen 597
683 chrisfen 739 dudr = - sw * preVal * riji*riji*f1
684    
685 chrisfen 597 endif
686    
687 chrisfen 739 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 411 endif
695    
696     if (j_is_Dipole) then
697 chrisfen 849 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 411
704 chrisfen 626 pref = pre12 * q_i * mu_j
705 gezelter 411
706 chrisfen 710 if (summationMethod .eq. REACTION_FIELD) then
707 chrisfen 700 ri2 = riji * riji
708     ri3 = ri2 * riji
709 chrisfen 696
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 597 else
729     if (j_is_SplitDipole) then
730 gezelter 938 BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
731     ri = 1.0d0 / BigR
732 chrisfen 597 scale = rij * ri
733     else
734     ri = riji
735 gezelter 938 scale = 1.0d0
736 chrisfen 597 endif
737    
738     ri2 = ri * ri
739     ri3 = ri2 * ri
740     sc2 = scale * scale
741 chrisfen 626
742 chrisfen 849 pot_term = ri2 * scale * f1
743     vterm = - pref * ct_j * pot_term
744 chrisfen 626 vpair = vpair + vterm
745     epot = epot + sw*vterm
746 chrisfen 597
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 849 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 411
762 chrisfen 597 endif
763 gezelter 411 endif
764 gezelter 421
765 gezelter 439 if (j_is_Quadrupole) then
766 chrisfen 849 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 439 ri2 = riji * riji
776     ri3 = ri2 * riji
777 gezelter 440 ri4 = ri2 * ri2
778 gezelter 439 cx2 = cx_j * cx_j
779     cy2 = cy_j * cy_j
780     cz2 = cz_j * cz_j
781    
782 gezelter 938 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 849 vterm = pref * (pot_term*f1 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f2)
787 chrisfen 740 vpair = vpair + vterm
788     epot = epot + sw*vterm
789 chrisfen 597
790 chrisfen 849 dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
791     sw*pref*ri4 * ( &
792 gezelter 938 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 849 + (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 938 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 849 + (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 938 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 849 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
808 chrisfen 740
809 gezelter 938 dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*xhat) &
810 chrisfen 849 * (3.0d0*f1 + f3) )
811 gezelter 938 dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*yhat) &
812 chrisfen 849 * (3.0d0*f1 + f3) )
813 gezelter 938 dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*zhat) &
814 chrisfen 849 * (3.0d0*f1 + f3) )
815 chrisfen 740
816 gezelter 938 duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*xhat) &
817 chrisfen 849 * (3.0d0*f1 + f3) )
818 gezelter 938 duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*yhat) &
819 chrisfen 849 * (3.0d0*f1 + f3) )
820 gezelter 938 duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*zhat) &
821 chrisfen 849 * (3.0d0*f1 + f3) )
822 chrisfen 740
823 gezelter 938 duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*xhat) &
824 chrisfen 849 * (3.0d0*f1 + f3) )
825 gezelter 938 duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*yhat) &
826 chrisfen 849 * (3.0d0*f1 + f3) )
827 gezelter 938 duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*zhat) &
828 chrisfen 849 * (3.0d0*f1 + f3) )
829    
830 gezelter 439 endif
831 gezelter 411 endif
832 chrisfen 740
833 gezelter 411 if (i_is_Dipole) then
834 gezelter 507
835 gezelter 411 if (j_is_Charge) then
836 chrisfen 849 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 626
843 chrisfen 740 pref = pre12 * q_j * mu_i
844    
845 chrisfen 719 if (summationMethod .eq. SHIFTED_POTENTIAL) then
846 chrisfen 597 ri2 = riji * riji
847     ri3 = ri2 * riji
848 chrisfen 719
849 chrisfen 849 pot_term = ri2*f1 - rcuti2*f1c
850 chrisfen 719 vterm = pref * ct_i * pot_term
851     vpair = vpair + vterm
852     epot = epot + sw*vterm
853    
854 chrisfen 849 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 719
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 507
862 chrisfen 719 elseif (summationMethod .eq. SHIFTED_FORCE) then
863     ri2 = riji * riji
864     ri3 = ri2 * riji
865    
866 chrisfen 849 !! 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 719 vterm = pref * ct_i * pot_term
870     vpair = vpair + vterm
871     epot = epot + sw*vterm
872    
873 chrisfen 849 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 719
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 849
884 chrisfen 719 elseif (summationMethod .eq. REACTION_FIELD) then
885     ri2 = riji * riji
886     ri3 = ri2 * riji
887    
888 chrisfen 700 vterm = pref * ct_i * ( ri2 - preRF2*rij )
889 chrisfen 696 vpair = vpair + vterm
890     epot = epot + sw*vterm
891    
892 chrisfen 700 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 696
899 chrisfen 700 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 696
903 chrisfen 597 else
904     if (i_is_SplitDipole) then
905 gezelter 938 BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
906     ri = 1.0d0 / BigR
907 chrisfen 597 scale = rij * ri
908     else
909 gezelter 421 ri = riji
910 gezelter 938 scale = 1.0d0
911 gezelter 421 endif
912 chrisfen 597
913     ri2 = ri * ri
914     ri3 = ri2 * ri
915     sc2 = scale * scale
916 chrisfen 626
917 chrisfen 849 pot_term = ri2 * f1 * scale
918     vterm = pref * ct_i * pot_term
919 chrisfen 626 vpair = vpair + vterm
920     epot = epot + sw*vterm
921 chrisfen 597
922 chrisfen 849 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 597
929 chrisfen 849 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 421 endif
933 chrisfen 597 endif
934 chrisfen 626
935 chrisfen 597 if (j_is_Dipole) then
936 chrisfen 849 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 719 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 421
953 chrisfen 710 if (summationMethod .eq. REACTION_FIELD) then
954 chrisfen 695 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 597 else
982     if (i_is_SplitDipole) then
983     if (j_is_SplitDipole) then
984 gezelter 938 BigR = sqrt(r2 + 0.25d0 * d_i * d_i + 0.25d0 * d_j * d_j)
985 chrisfen 597 else
986 gezelter 938 BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
987 chrisfen 597 endif
988 gezelter 938 ri = 1.0d0 / BigR
989 chrisfen 597 scale = rij * ri
990     else
991     if (j_is_SplitDipole) then
992 gezelter 938 BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
993     ri = 1.0d0 / BigR
994 chrisfen 597 scale = rij * ri
995     else
996     ri = riji
997 gezelter 938 scale = 1.0d0
998 chrisfen 597 endif
999     endif
1000    
1001     sc2 = scale * scale
1002 chrisfen 719
1003 chrisfen 849 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 626 vpair = vpair + vterm
1006     epot = epot + sw*vterm
1007 chrisfen 597
1008 chrisfen 849 f13 = f1+f3
1009     f134 = f13 + f4
1010 chrisfen 597
1011 chrisfen 849 !!$ 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 597
1037 chrisfen 849 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 597
1044 chrisfen 849 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 597 endif
1051 gezelter 411 endif
1052     endif
1053 gezelter 439
1054     if (i_is_Quadrupole) then
1055     if (j_is_Charge) then
1056 chrisfen 849 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 439 ri2 = riji * riji
1066     ri3 = ri2 * riji
1067 gezelter 440 ri4 = ri2 * ri2
1068 gezelter 439 cx2 = cx_i * cx_i
1069     cy2 = cy_i * cy_i
1070     cz2 = cz_i * cz_i
1071 gezelter 507
1072 gezelter 938 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 849 vterm = pref * (pot_term*f1 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f2)
1077 chrisfen 740 vpair = vpair + vterm
1078     epot = epot + sw*vterm
1079    
1080 chrisfen 849 dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
1081     sw*pref*ri4 * ( &
1082 gezelter 938 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 849 + (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 938 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 849 + (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 938 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 849 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1098 chrisfen 740
1099 gezelter 938 dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*xhat) &
1100 chrisfen 849 * (3.0d0*f1 + f3) )
1101 gezelter 938 dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*yhat) &
1102 chrisfen 849 * (3.0d0*f1 + f3) )
1103 gezelter 938 dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*zhat) &
1104 chrisfen 849 * (3.0d0*f1 + f3) )
1105 chrisfen 740
1106 gezelter 938 duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*xhat) &
1107 chrisfen 849 * (3.0d0*f1 + f3) )
1108 gezelter 938 duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*yhat) &
1109 chrisfen 849 * (3.0d0*f1 + f3) )
1110 gezelter 938 duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*zhat) &
1111 chrisfen 849 * (3.0d0*f1 + f3) )
1112 chrisfen 740
1113 gezelter 938 duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*xhat) &
1114 chrisfen 849 * (3.0d0*f1 + f3) )
1115 gezelter 938 duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*yhat) &
1116 chrisfen 849 * (3.0d0*f1 + f3) )
1117 gezelter 938 duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*zhat) &
1118 chrisfen 849 * (3.0d0*f1 + f3) )
1119 chrisfen 740
1120 gezelter 439 endif
1121     endif
1122 gezelter 507
1123    
1124 gezelter 411 if (do_pot) then
1125     #ifdef IS_MPI
1126 chuckv 656 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 411 #else
1129     pot = pot + epot
1130     #endif
1131     endif
1132 gezelter 507
1133 gezelter 411 #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 507
1138 gezelter 411 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 507
1142 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1143 gezelter 439 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 411 endif
1147 gezelter 439 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 411
1152 gezelter 439 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 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1158 gezelter 439 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 411 endif
1162 gezelter 439 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 411
1167 gezelter 439 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 411 #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 507
1177 gezelter 411 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 507
1181 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1182 gezelter 439 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 411 endif
1186 gezelter 439 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 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1197 gezelter 439 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 411 endif
1201 gezelter 439 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 411 #endif
1212 gezelter 507
1213 gezelter 411 #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 507
1223 gezelter 411 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 492
1232     subroutine destroyElectrostaticTypes()
1233    
1234 gezelter 507 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1235    
1236 chuckv 492 end subroutine destroyElectrostaticTypes
1237    
1238 chrisfen 703 subroutine self_self(atom1, eFrame, mypot, t, do_pot)
1239 chrisfen 695 logical, intent(in) :: do_pot
1240 chrisfen 682 integer, intent(in) :: atom1
1241 chrisfen 695 integer :: atid1
1242 chrisfen 682 real(kind=dp), dimension(9,nLocal) :: eFrame
1243 chrisfen 695 real(kind=dp), dimension(3,nLocal) :: t
1244 chrisfen 703 real(kind=dp) :: mu1, c1
1245     real(kind=dp) :: preVal, epot, mypot
1246 chrisfen 695 real(kind=dp) :: eix, eiy, eiz
1247 chrisfen 682
1248 chrisfen 695 ! this is a local only array, so we use the local atom type id's:
1249     atid1 = atid(atom1)
1250 chrisfen 703
1251     if (.not.summationMethodChecked) then
1252     call checkSummationMethod()
1253     endif
1254 chrisfen 695
1255 chrisfen 703 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 743 elseif ( (summationMethod .eq. SHIFTED_FORCE) .or. &
1281     (summationMethod .eq. SHIFTED_POTENTIAL) ) then
1282 chrisfen 717 if (ElectrostaticMap(atid1)%is_Charge) then
1283     c1 = getCharge(atid1)
1284    
1285     if (screeningMethod .eq. DAMPED) then
1286 gezelter 938 mypot = mypot - (f0c * rcuti * 0.5d0 + &
1287 chrisfen 717 dampingAlpha*invRootPi) * c1 * c1
1288    
1289     else
1290 gezelter 938 mypot = mypot - (rcuti * 0.5d0 * c1 * c1)
1291 chrisfen 717
1292     endif
1293     endif
1294     endif
1295 chrisfen 695
1296 chrisfen 682 return
1297 chrisfen 703 end subroutine self_self
1298 chrisfen 682
1299 chrisfen 703 subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, myPot, &
1300 chrisfen 700 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 703 real(kind=dp) :: pref, preVal, rfVal, myPot
1324 chrisfen 700 real(kind=dp) :: dudx, dudy, dudz, dudr
1325    
1326 chrisfen 703 if (.not.summationMethodChecked) then
1327     call checkSummationMethod()
1328 chrisfen 700 endif
1329    
1330 gezelter 939 dudx = zero
1331     dudy = zero
1332     dudz = zero
1333 chrisfen 700
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 703 myPot = myPot + sw*vterm
1357    
1358 chrisfen 700 dudr = sw*preVal * 2.0d0*rfVal*riji
1359 chrisfen 703
1360 chrisfen 700 dudx = dudx + dudr * xhat
1361     dudy = dudy + dudr * yhat
1362     dudz = dudz + dudr * zhat
1363 chrisfen 703
1364 chrisfen 700 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 703
1372 chrisfen 700 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 703 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 700 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 703
1390 chrisfen 700 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 703
1398 chrisfen 700 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 703 myPot = myPot + sw*vterm
1404 chrisfen 700
1405 chrisfen 703 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 700
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 703
1418    
1419     ! accumulate the forces and torques resulting from the self term
1420 chrisfen 700 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 411 end module electrostatic_module