ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 695
Committed: Sun Oct 23 21:08:08 2005 UTC (19 years, 7 months ago) by chrisfen
File size: 50256 byte(s)
Log Message: 
streamlined reaction field for dipoles (now a good bit faster) and added reaction field for charges - still need to do charge-dipole RF

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