ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 941
Committed: Fri Apr 21 03:19:52 2006 UTC (19 years, 1 month ago) by chrisfen
File size: 53991 byte(s)
Log Message: 
Electrosplines added...

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