ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 696
Committed: Mon Oct 24 14:06:36 2005 UTC (19 years, 7 months ago) by chrisfen
File size: 52352 byte(s)
Log Message: 
added charge-dipole reaction field - don't know if it works...

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 696 elseif (summationMethod .eq. REACTION_FIELD) then
730     ri2 = ri * ri
731     ri3 = ri2 * ri
732    
733     pref = pre12 * q_i * mu_j
734     vterm = - pref * ct_j * ( ri2 - preRF2*rij )
735     vpair = vpair + vterm
736     epot = epot + sw*vterm
737    
738     !! this has a + sign in the () because the rij vector is
739     !! r_j - r_i and the charge-dipole potential takes the origin
740     !! as the point dipole, which is atom j in this case.
741    
742     dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
743     preRF2*uz_j(1) )
744     dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
745     preRF2*uz_j(2) )
746     dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
747     preRF2*uz_j(3) )
748     duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
749     duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
750     duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
751    
752 chrisfen 597 else
753     if (j_is_SplitDipole) then
754     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
755     ri = 1.0_dp / BigR
756     scale = rij * ri
757     else
758     ri = riji
759     scale = 1.0_dp
760     endif
761    
762     ri2 = ri * ri
763     ri3 = ri2 * ri
764     sc2 = scale * scale
765 chrisfen 626
766     pref = pre12 * q_i * mu_j
767 chrisfen 597 vterm = - pref * ct_j * ri2 * scale
768 chrisfen 626 vpair = vpair + vterm
769     epot = epot + sw*vterm
770 chrisfen 597
771     !! this has a + sign in the () because the rij vector is
772     !! r_j - r_i and the charge-dipole potential takes the origin
773     !! as the point dipole, which is atom j in this case.
774    
775 chrisfen 626 dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
776     dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
777     dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
778 chrisfen 597
779 chrisfen 626 duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
780     duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
781     duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
782 gezelter 411
783 chrisfen 597 endif
784 gezelter 411 endif
785 gezelter 421
786 gezelter 439 if (j_is_Quadrupole) then
787     ri2 = riji * riji
788     ri3 = ri2 * riji
789 gezelter 440 ri4 = ri2 * ri2
790 gezelter 439 cx2 = cx_j * cx_j
791     cy2 = cy_j * cy_j
792     cz2 = cz_j * cz_j
793    
794 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
795 chrisfen 626 pref = pre14 * q_i / 3.0_dp
796 chrisfen 597 vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
797     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
798     qzz_j * (3.0_dp*cz2 - 1.0_dp) )
799     vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
800     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
801     qzz_j * (3.0_dp*cz2 - 1.0_dp) )
802 chrisfen 626 vpair = vpair + ( vterm1 - vterm2 )
803     epot = epot + sw*( vterm1 - vterm2 )
804 chrisfen 597
805     dudx = dudx - (5.0_dp * &
806 chrisfen 626 (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
807 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
808     qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
809     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
810     qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
811     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
812     qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
813     dudy = dudy - (5.0_dp * &
814 chrisfen 626 (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
815 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
816     qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
817     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
818     qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
819     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
820     qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
821     dudz = dudz - (5.0_dp * &
822 chrisfen 626 (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
823 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
824     qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
825     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
826     qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
827     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
828     qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
829    
830 chrisfen 626 dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
831 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
832 chrisfen 626 dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
833 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
834 chrisfen 626 dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
835 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
836    
837 chrisfen 626 duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
838 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
839 chrisfen 626 duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
840 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
841 chrisfen 626 duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
842 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
843    
844 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
845 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
846 chrisfen 626 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
847 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
848 chrisfen 626 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
849 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
850    
851     else
852 chrisfen 626 pref = pre14 * q_i / 3.0_dp
853 chrisfen 597 vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
854     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
855     qzz_j * (3.0_dp*cz2 - 1.0_dp))
856 chrisfen 626 vpair = vpair + vterm
857     epot = epot + sw*vterm
858 chrisfen 597
859 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
860 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
861     qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
862     qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
863 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
864 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
865     qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
866     qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
867 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
868 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
869     qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
870     qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
871    
872 chrisfen 626 dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
873     dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
874     dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
875 chrisfen 597
876 chrisfen 626 duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
877     duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
878     duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
879 chrisfen 597
880 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
881     duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
882     duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
883 chrisfen 597
884     endif
885 gezelter 439 endif
886 gezelter 411 endif
887 gezelter 507
888 gezelter 411 if (i_is_Dipole) then
889 gezelter 507
890 gezelter 411 if (j_is_Charge) then
891 chrisfen 626
892     pref = pre12 * q_j * mu_i
893    
894 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
895 chrisfen 597 ri2 = riji * riji
896     ri3 = ri2 * riji
897 gezelter 507
898 chrisfen 626 pref = pre12 * q_j * mu_i
899 chrisfen 597 vterm = pref * ct_i * (ri2 - rcuti2)
900 chrisfen 626 vpair = vpair + vterm
901     epot = epot + sw*vterm
902 chrisfen 597
903     !! this has a + sign in the () because the rij vector is
904     !! r_j - r_i and the charge-dipole potential takes the origin
905     !! as the point dipole, which is atom j in this case.
906    
907 chrisfen 626 dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
908 chrisfen 597 - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
909 chrisfen 626 dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
910 chrisfen 597 - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
911 chrisfen 626 dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
912 chrisfen 597 - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
913    
914 chrisfen 626 duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
915     duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
916     duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
917 gezelter 411
918 chrisfen 696 elseif (summationMethod .eq. REACTION_FIELD) then
919     ri2 = ri * ri
920     ri3 = ri2 * ri
921    
922     pref = pre12 * q_j * mu_i
923     vterm = pref * ct_i * ( ri2 - preRF*rij )
924     vpair = vpair + vterm
925     epot = epot + sw*vterm
926    
927     dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0*ct_i*xhat - &
928     preRF*uz_i(1) )
929     dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0*ct_i*yhat - &
930     preRF*uz_i(2) )
931     dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0*ct_i*zhat - &
932     preRF*uz_i(3) )
933    
934     duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF*rij )
935     duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF*rij )
936     duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF*rij )
937    
938 chrisfen 597 else
939     if (i_is_SplitDipole) then
940 gezelter 421 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
941     ri = 1.0_dp / BigR
942 chrisfen 597 scale = rij * ri
943     else
944 gezelter 421 ri = riji
945     scale = 1.0_dp
946     endif
947 chrisfen 597
948     ri2 = ri * ri
949     ri3 = ri2 * ri
950     sc2 = scale * scale
951 chrisfen 626
952     pref = pre12 * q_j * mu_i
953 chrisfen 597 vterm = pref * ct_i * ri2 * scale
954 chrisfen 626 vpair = vpair + vterm
955     epot = epot + sw*vterm
956 chrisfen 597
957 chrisfen 626 dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
958     dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
959     dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
960 chrisfen 597
961 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
962     duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
963     duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
964 gezelter 421 endif
965 chrisfen 597 endif
966 chrisfen 626
967 chrisfen 597 if (j_is_Dipole) then
968 gezelter 421
969 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
970 chrisfen 597 ri2 = riji * riji
971     ri3 = ri2 * riji
972     ri4 = ri2 * ri2
973 gezelter 507
974 chrisfen 626 pref = pre22 * mu_i * mu_j
975 chrisfen 597 vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
976 chrisfen 626 vpair = vpair + vterm
977     epot = epot + sw*vterm
978 chrisfen 597
979     a1 = 5.0d0 * ct_i * ct_j - ct_ij
980    
981 chrisfen 626 dudx = dudx + sw*pref*3.0d0*ri4 &
982     * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
983     - sw*pref*3.0d0*rcuti4 &
984     * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
985     dudy = dudy + sw*pref*3.0d0*ri4 &
986     * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
987     - sw*pref*3.0d0*rcuti4 &
988     * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
989     dudz = dudz + sw*pref*3.0d0*ri4 &
990     * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
991     - sw*pref*3.0d0*rcuti4 &
992     * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
993 chrisfen 597
994 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
995 chrisfen 597 - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
996 chrisfen 626 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
997 chrisfen 597 - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
998 chrisfen 626 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
999 chrisfen 597 - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
1000 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1001 chrisfen 597 - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
1002 chrisfen 626 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1003 chrisfen 597 - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
1004 chrisfen 626 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1005 chrisfen 597 - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
1006 chrisfen 626
1007 chrisfen 695 elseif (summationMethod .eq. REACTION_FIELD) then
1008     ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1009    
1010     ri2 = riji * riji
1011     ri3 = ri2 * riji
1012     ri4 = ri2 * ri2
1013    
1014     pref = pre22 * mu_i * mu_j
1015    
1016     vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1017     preRF2*ct_ij )
1018     vpair = vpair + vterm
1019     epot = epot + sw*vterm
1020    
1021     a1 = 5.0d0 * ct_i * ct_j - ct_ij
1022    
1023     dudx = dudx + sw*pref*3.0d0*ri4 &
1024     * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1025     dudy = dudy + sw*pref*3.0d0*ri4 &
1026     * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1027     dudz = dudz + sw*pref*3.0d0*ri4 &
1028     * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1029    
1030     duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1031     - preRF2*uz_j(1))
1032     duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1033     - preRF2*uz_j(2))
1034     duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1035     - preRF2*uz_j(3))
1036     duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1037     - preRF2*uz_i(1))
1038     duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1039     - preRF2*uz_i(2))
1040     duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1041     - preRF2*uz_i(3))
1042    
1043 chrisfen 597 else
1044     if (i_is_SplitDipole) then
1045     if (j_is_SplitDipole) then
1046     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
1047     else
1048     BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
1049     endif
1050     ri = 1.0_dp / BigR
1051     scale = rij * ri
1052     else
1053     if (j_is_SplitDipole) then
1054     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
1055     ri = 1.0_dp / BigR
1056     scale = rij * ri
1057     else
1058     ri = riji
1059     scale = 1.0_dp
1060     endif
1061     endif
1062    
1063     ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1064    
1065     ri2 = ri * ri
1066     ri3 = ri2 * ri
1067     ri4 = ri2 * ri2
1068     sc2 = scale * scale
1069    
1070 chrisfen 626 pref = pre22 * mu_i * mu_j
1071 chrisfen 597 vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1072 chrisfen 626 vpair = vpair + vterm
1073     epot = epot + sw*vterm
1074 chrisfen 597
1075     a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1076    
1077 chrisfen 626 dudx = dudx + sw*pref*3.0d0*ri4*scale &
1078     *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1079     dudy = dudy + sw*pref*3.0d0*ri4*scale &
1080     *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1081     dudz = dudz + sw*pref*3.0d0*ri4*scale &
1082     *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1083 chrisfen 597
1084 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1085     *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1086     duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1087     *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1088     duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1089     *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1090 chrisfen 597
1091 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1092     *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1093     duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1094     *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1095     duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1096     *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1097 chrisfen 597 endif
1098 gezelter 411 endif
1099     endif
1100 gezelter 439
1101     if (i_is_Quadrupole) then
1102     if (j_is_Charge) then
1103 gezelter 507
1104 gezelter 439 ri2 = riji * riji
1105     ri3 = ri2 * riji
1106 gezelter 440 ri4 = ri2 * ri2
1107 gezelter 439 cx2 = cx_i * cx_i
1108     cy2 = cy_i * cy_i
1109     cz2 = cz_i * cz_i
1110 gezelter 507
1111 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
1112 chrisfen 626 pref = pre14 * q_j / 3.0_dp
1113 chrisfen 597 vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1114     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1115     qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1116     vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1117     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1118     qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1119 chrisfen 626 vpair = vpair + ( vterm1 - vterm2 )
1120     epot = epot + sw*( vterm1 - vterm2 )
1121 chrisfen 597
1122 chrisfen 626 dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1123     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1124 chrisfen 597 qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1125     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1126     qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1127     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1128     qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1129 chrisfen 626 dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1130     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1131 chrisfen 597 qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1132     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1133     qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1134     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1135     qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1136 chrisfen 626 dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1137     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1138 chrisfen 597 qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1139     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1140     qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1141     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1142     qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1143    
1144 chrisfen 626 dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1145 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1146 chrisfen 626 dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1147 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1148 chrisfen 626 dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1149 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1150    
1151 chrisfen 626 duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1152 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1153 chrisfen 626 duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1154 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1155 chrisfen 626 duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1156 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1157    
1158 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1159 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1160 chrisfen 626 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1161 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1162 chrisfen 626 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1163 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1164 gezelter 507
1165 chrisfen 597 else
1166 chrisfen 626 pref = pre14 * q_j / 3.0_dp
1167 chrisfen 597 vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1168     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1169     qzz_i * (3.0_dp*cz2 - 1.0_dp))
1170 chrisfen 626 vpair = vpair + vterm
1171     epot = epot + sw*vterm
1172 chrisfen 597
1173 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1174 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1175     qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1176     qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1177 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1178 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1179     qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1180     qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1181 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1182 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1183     qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1184     qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1185    
1186 chrisfen 626 dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1187     dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1188     dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1189 chrisfen 597
1190 chrisfen 626 duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1191     duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1192     duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1193 chrisfen 597
1194 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1195     duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1196     duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1197 chrisfen 597 endif
1198 gezelter 439 endif
1199     endif
1200 gezelter 507
1201    
1202 gezelter 411 if (do_pot) then
1203     #ifdef IS_MPI
1204 chuckv 656 pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1205     pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1206 gezelter 411 #else
1207     pot = pot + epot
1208     #endif
1209     endif
1210 gezelter 507
1211 gezelter 411 #ifdef IS_MPI
1212     f_Row(1,atom1) = f_Row(1,atom1) + dudx
1213     f_Row(2,atom1) = f_Row(2,atom1) + dudy
1214     f_Row(3,atom1) = f_Row(3,atom1) + dudz
1215 gezelter 507
1216 gezelter 411 f_Col(1,atom2) = f_Col(1,atom2) - dudx
1217     f_Col(2,atom2) = f_Col(2,atom2) - dudy
1218     f_Col(3,atom2) = f_Col(3,atom2) - dudz
1219 gezelter 507
1220 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1221 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1222     t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1223     t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1224 gezelter 411 endif
1225 gezelter 439 if (i_is_Quadrupole) then
1226     t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1227     t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1228     t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1229 gezelter 411
1230 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1231     t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1232     t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1233     endif
1234    
1235 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1236 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1237     t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1238     t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1239 gezelter 411 endif
1240 gezelter 439 if (j_is_Quadrupole) then
1241     t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1242     t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1243     t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1244 gezelter 411
1245 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1246     t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1247     t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1248     endif
1249    
1250 gezelter 411 #else
1251     f(1,atom1) = f(1,atom1) + dudx
1252     f(2,atom1) = f(2,atom1) + dudy
1253     f(3,atom1) = f(3,atom1) + dudz
1254 gezelter 507
1255 gezelter 411 f(1,atom2) = f(1,atom2) - dudx
1256     f(2,atom2) = f(2,atom2) - dudy
1257     f(3,atom2) = f(3,atom2) - dudz
1258 gezelter 507
1259 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1260 gezelter 439 t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1261     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1262     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1263 gezelter 411 endif
1264 gezelter 439 if (i_is_Quadrupole) then
1265     t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1266     t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1267     t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1268    
1269     t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1270     t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1271     t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1272     endif
1273    
1274 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1275 gezelter 439 t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1276     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1277     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1278 gezelter 411 endif
1279 gezelter 439 if (j_is_Quadrupole) then
1280     t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1281     t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1282     t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1283    
1284     t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1285     t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1286     t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1287     endif
1288    
1289 gezelter 411 #endif
1290 gezelter 507
1291 gezelter 411 #ifdef IS_MPI
1292     id1 = AtomRowToGlobal(atom1)
1293     id2 = AtomColToGlobal(atom2)
1294     #else
1295     id1 = atom1
1296     id2 = atom2
1297     #endif
1298    
1299     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1300 gezelter 507
1301 gezelter 411 fpair(1) = fpair(1) + dudx
1302     fpair(2) = fpair(2) + dudy
1303     fpair(3) = fpair(3) + dudz
1304    
1305     endif
1306    
1307     return
1308     end subroutine doElectrostaticPair
1309 chuckv 492
1310     subroutine destroyElectrostaticTypes()
1311    
1312 gezelter 507 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1313    
1314 chuckv 492 end subroutine destroyElectrostaticTypes
1315    
1316 chrisfen 695 subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1317     logical, intent(in) :: do_pot
1318 chrisfen 682 integer, intent(in) :: atom1
1319 chrisfen 695 integer :: atid1
1320 chrisfen 682 real(kind=dp), dimension(9,nLocal) :: eFrame
1321 chrisfen 695 real(kind=dp), dimension(3,nLocal) :: t
1322     real(kind=dp) :: mu1
1323     real(kind=dp) :: preVal, epot, rfpot
1324     real(kind=dp) :: eix, eiy, eiz
1325 chrisfen 682
1326 chrisfen 695 ! this is a local only array, so we use the local atom type id's:
1327     atid1 = atid(atom1)
1328    
1329     if (ElectrostaticMap(atid1)%is_Dipole) then
1330     mu1 = getDipoleMoment(atid1)
1331    
1332     preVal = pre22 * preRF2 * mu1*mu1
1333     rfpot = rfpot - 0.5d0*preVal
1334 chrisfen 682
1335 chrisfen 695 ! The self-correction term adds into the reaction field vector
1336    
1337     eix = preVal * eFrame(3,atom1)
1338     eiy = preVal * eFrame(6,atom1)
1339     eiz = preVal * eFrame(9,atom1)
1340 chrisfen 682
1341 chrisfen 695 ! once again, this is self-self, so only the local arrays are needed
1342     ! even for MPI jobs:
1343    
1344     t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1345     eFrame(9,atom1)*eiy
1346     t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1347     eFrame(3,atom1)*eiz
1348     t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1349     eFrame(6,atom1)*eix
1350 chrisfen 682
1351     endif
1352 chrisfen 695
1353 chrisfen 682 return
1354 chrisfen 695 end subroutine rf_self_self
1355 chrisfen 682
1356 gezelter 411 end module electrostatic_module