ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 739
Committed: Tue Nov 15 19:04:02 2005 UTC (19 years, 6 months ago) by chrisfen
File size: 58866 byte(s)
Log Message: 
cleaned up the charge-charge interactions a bit...

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 chrisfen 716 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
62 gezelter 602
63 chuckv 656
64 gezelter 434 !! these prefactors convert the multipole interactions into kcal / mol
65     !! all were computed assuming distances are measured in angstroms
66     !! Charge-Charge, assuming charges are measured in electrons
67 gezelter 411 real(kind=dp), parameter :: pre11 = 332.0637778_dp
68 gezelter 434 !! Charge-Dipole, assuming charges are measured in electrons, and
69     !! dipoles are measured in debyes
70     real(kind=dp), parameter :: pre12 = 69.13373_dp
71     !! Dipole-Dipole, assuming dipoles are measured in debyes
72     real(kind=dp), parameter :: pre22 = 14.39325_dp
73     !! Charge-Quadrupole, assuming charges are measured in electrons, and
74     !! quadrupoles are measured in 10^-26 esu cm^2
75     !! This unit is also known affectionately as an esu centi-barn.
76     real(kind=dp), parameter :: pre14 = 69.13373_dp
77 gezelter 411
78 chrisfen 712 !! variables to handle different summation methods for long-range
79     !! electrostatics:
80 gezelter 602 integer, save :: summationMethod = NONE
81 chrisfen 710 integer, save :: screeningMethod = UNDAMPED
82 chrisfen 603 logical, save :: summationMethodChecked = .false.
83 gezelter 602 real(kind=DP), save :: defaultCutoff = 0.0_DP
84 chrisfen 682 real(kind=DP), save :: defaultCutoff2 = 0.0_DP
85 gezelter 602 logical, save :: haveDefaultCutoff = .false.
86     real(kind=DP), save :: dampingAlpha = 0.0_DP
87 chrisfen 716 real(kind=DP), save :: alpha2 = 0.0_DP
88 gezelter 602 logical, save :: haveDampingAlpha = .false.
89 chrisfen 682 real(kind=DP), save :: dielectric = 1.0_DP
90 gezelter 602 logical, save :: haveDielectric = .false.
91     real(kind=DP), save :: constEXP = 0.0_DP
92 chrisfen 682 real(kind=dp), save :: rcuti = 0.0_DP
93     real(kind=dp), save :: rcuti2 = 0.0_DP
94     real(kind=dp), save :: rcuti3 = 0.0_DP
95     real(kind=dp), save :: rcuti4 = 0.0_DP
96     real(kind=dp), save :: alphaPi = 0.0_DP
97     real(kind=dp), save :: invRootPi = 0.0_DP
98     real(kind=dp), save :: rrf = 1.0_DP
99     real(kind=dp), save :: rt = 1.0_DP
100     real(kind=dp), save :: rrfsq = 1.0_DP
101     real(kind=dp), save :: preRF = 0.0_DP
102 chrisfen 695 real(kind=dp), save :: preRF2 = 0.0_DP
103 chrisfen 716 real(kind=dp), save :: f0 = 1.0_DP
104     real(kind=dp), save :: f1 = 1.0_DP
105     real(kind=dp), save :: f2 = 0.0_DP
106     real(kind=dp), save :: f0c = 1.0_DP
107     real(kind=dp), save :: f1c = 1.0_DP
108     real(kind=dp), save :: f2c = 0.0_DP
109    
110 chuckv 632 #ifdef __IFC
111     ! error function for ifc version > 7.
112 chuckv 631 double precision, external :: derfc
113 chuckv 632 #endif
114    
115 gezelter 602 public :: setElectrostaticSummationMethod
116 chrisfen 710 public :: setScreeningMethod
117 gezelter 602 public :: setElectrostaticCutoffRadius
118 chrisfen 710 public :: setDampingAlpha
119 gezelter 602 public :: setReactionFieldDielectric
120 gezelter 411 public :: newElectrostaticType
121     public :: setCharge
122     public :: setDipoleMoment
123     public :: setSplitDipoleDistance
124     public :: setQuadrupoleMoments
125     public :: doElectrostaticPair
126     public :: getCharge
127     public :: getDipoleMoment
128 chuckv 492 public :: destroyElectrostaticTypes
129 chrisfen 703 public :: self_self
130 chrisfen 700 public :: rf_self_excludes
131 gezelter 411
132     type :: Electrostatic
133     integer :: c_ident
134     logical :: is_Charge = .false.
135     logical :: is_Dipole = .false.
136     logical :: is_SplitDipole = .false.
137     logical :: is_Quadrupole = .false.
138 chrisfen 532 logical :: is_Tap = .false.
139 gezelter 411 real(kind=DP) :: charge = 0.0_DP
140     real(kind=DP) :: dipole_moment = 0.0_DP
141     real(kind=DP) :: split_dipole_distance = 0.0_DP
142     real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP
143     end type Electrostatic
144    
145     type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
146    
147     contains
148    
149 gezelter 602 subroutine setElectrostaticSummationMethod(the_ESM)
150     integer, intent(in) :: the_ESM
151    
152     if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
153     call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
154     endif
155    
156 chrisfen 610 summationMethod = the_ESM
157 chrisfen 626
158 gezelter 602 end subroutine setElectrostaticSummationMethod
159    
160 chrisfen 710 subroutine setScreeningMethod(the_SM)
161     integer, intent(in) :: the_SM
162     screeningMethod = the_SM
163     end subroutine setScreeningMethod
164    
165 chrisfen 682 subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
166 gezelter 602 real(kind=dp), intent(in) :: thisRcut
167 chrisfen 682 real(kind=dp), intent(in) :: thisRsw
168 gezelter 602 defaultCutoff = thisRcut
169 chrisfen 682 rrf = defaultCutoff
170     rt = thisRsw
171 gezelter 602 haveDefaultCutoff = .true.
172     end subroutine setElectrostaticCutoffRadius
173    
174 chrisfen 710 subroutine setDampingAlpha(thisAlpha)
175 gezelter 602 real(kind=dp), intent(in) :: thisAlpha
176     dampingAlpha = thisAlpha
177 chrisfen 716 alpha2 = dampingAlpha*dampingAlpha
178 gezelter 602 haveDampingAlpha = .true.
179 chrisfen 710 end subroutine setDampingAlpha
180 gezelter 602
181     subroutine setReactionFieldDielectric(thisDielectric)
182     real(kind=dp), intent(in) :: thisDielectric
183     dielectric = thisDielectric
184     haveDielectric = .true.
185     end subroutine setReactionFieldDielectric
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 chrisfen 710 if (screeningMethod .eq. DAMPED) then
409 chrisfen 703 if (.not.haveDampingAlpha) then
410     call handleError("checkSummationMethod", "no Damping Alpha set!")
411     endif
412    
413     if (.not.haveDefaultCutoff) then
414     call handleError("checkSummationMethod", "no Default Cutoff set!")
415     endif
416 chrisfen 603
417 chrisfen 716 constEXP = exp(-alpha2*defaultCutoff*defaultCutoff)
418 chrisfen 703 invRootPi = 0.56418958354775628695d0
419 chrisfen 716 alphaPi = 2.0d0*dampingAlpha*invRootPi
420     f0c = derfc(dampingAlpha*defaultCutoff)
421     f1c = alphaPi*defaultCutoff*constEXP + f0c
422     f2c = alphaPi*2.0d0*alpha2*constEXP*rcuti2
423    
424 gezelter 602 endif
425    
426 chrisfen 603 if (summationMethod .eq. REACTION_FIELD) then
427 chrisfen 703 if (haveDielectric) then
428     defaultCutoff2 = defaultCutoff*defaultCutoff
429     preRF = (dielectric-1.0d0) / &
430     ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
431     preRF2 = 2.0d0*preRF
432     else
433     call handleError("checkSummationMethod", "Dielectric not set")
434 chrisfen 603 endif
435 chrisfen 703
436 chrisfen 603 endif
437    
438     summationMethodChecked = .true.
439 gezelter 602 end subroutine checkSummationMethod
440    
441 chrisfen 712
442 gezelter 411 subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
443 chrisfen 712 vpair, fpair, pot, eFrame, f, t, do_pot)
444 gezelter 507
445 chrisfen 703 logical, intent(in) :: do_pot
446 gezelter 507
447 gezelter 411 integer, intent(in) :: atom1, atom2
448     integer :: localError
449    
450 chrisfen 607 real(kind=dp), intent(in) :: rij, r2, sw
451 gezelter 411 real(kind=dp), intent(in), dimension(3) :: d
452     real(kind=dp), intent(inout) :: vpair
453 chrisfen 703 real(kind=dp), intent(inout), dimension(3) :: fpair
454 gezelter 411
455 chrisfen 626 real( kind = dp ) :: pot
456 gezelter 411 real( kind = dp ), dimension(9,nLocal) :: eFrame
457     real( kind = dp ), dimension(3,nLocal) :: f
458 chrisfen 710 real( kind = dp ), dimension(3,nLocal) :: felec
459 gezelter 411 real( kind = dp ), dimension(3,nLocal) :: t
460 gezelter 507
461 gezelter 439 real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
462     real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
463     real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
464     real (kind = dp), dimension(3) :: dudux_j, duduy_j, duduz_j
465 gezelter 411
466     logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
467     logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
468 chrisfen 532 logical :: i_is_Tap, j_is_Tap
469 gezelter 411 integer :: me1, me2, id1, id2
470     real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
471 gezelter 439 real (kind=dp) :: qxx_i, qyy_i, qzz_i
472     real (kind=dp) :: qxx_j, qyy_j, qzz_j
473     real (kind=dp) :: cx_i, cy_i, cz_i
474     real (kind=dp) :: cx_j, cy_j, cz_j
475     real (kind=dp) :: cx2, cy2, cz2
476 chrisfen 719 real (kind=dp) :: ct_i, ct_j, ct_ij, a0, a1
477 gezelter 421 real (kind=dp) :: riji, ri, ri2, ri3, ri4
478 chrisfen 597 real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
479 gezelter 421 real (kind=dp) :: xhat, yhat, zhat
480 gezelter 411 real (kind=dp) :: dudx, dudy, dudz
481 chrisfen 626 real (kind=dp) :: scale, sc2, bigR
482 chrisfen 716 real (kind=dp) :: varEXP
483 chrisfen 719 real (kind=dp) :: pot_term
484 chrisfen 695 real (kind=dp) :: preVal, rfVal
485 gezelter 411
486     if (.not.allocated(ElectrostaticMap)) then
487     call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
488     return
489     end if
490    
491 gezelter 602 if (.not.summationMethodChecked) then
492     call checkSummationMethod()
493     endif
494    
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 chrisfen 710 !!$ if (rij .ge. defaultCutoff) then
504     !!$ write(*,*) 'warning: rij = ', rij, ' rcut = ', defaultCutoff, ' sw = ', sw
505     !!$ endif
506    
507 gezelter 411 !! some variables we'll need independent of electrostatic type:
508    
509     riji = 1.0d0 / rij
510 chrisfen 644
511 gezelter 421 xhat = d(1) * riji
512     yhat = d(2) * riji
513     zhat = d(3) * riji
514 gezelter 411
515     !! logicals
516     i_is_Charge = ElectrostaticMap(me1)%is_Charge
517     i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
518     i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
519     i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
520 chrisfen 532 i_is_Tap = ElectrostaticMap(me1)%is_Tap
521 gezelter 411
522     j_is_Charge = ElectrostaticMap(me2)%is_Charge
523     j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
524     j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
525     j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
526 chrisfen 532 j_is_Tap = ElectrostaticMap(me2)%is_Tap
527 gezelter 411
528     if (i_is_Charge) then
529     q_i = ElectrostaticMap(me1)%charge
530     endif
531 gezelter 507
532 gezelter 411 if (i_is_Dipole) then
533     mu_i = ElectrostaticMap(me1)%dipole_moment
534     #ifdef IS_MPI
535 gezelter 439 uz_i(1) = eFrame_Row(3,atom1)
536     uz_i(2) = eFrame_Row(6,atom1)
537     uz_i(3) = eFrame_Row(9,atom1)
538 gezelter 411 #else
539 gezelter 439 uz_i(1) = eFrame(3,atom1)
540     uz_i(2) = eFrame(6,atom1)
541     uz_i(3) = eFrame(9,atom1)
542 gezelter 411 #endif
543 gezelter 439 ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
544 gezelter 411
545     if (i_is_SplitDipole) then
546     d_i = ElectrostaticMap(me1)%split_dipole_distance
547     endif
548 gezelter 507
549 gezelter 411 endif
550    
551 gezelter 439 if (i_is_Quadrupole) then
552     qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1)
553     qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2)
554     qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3)
555     #ifdef IS_MPI
556     ux_i(1) = eFrame_Row(1,atom1)
557     ux_i(2) = eFrame_Row(4,atom1)
558     ux_i(3) = eFrame_Row(7,atom1)
559     uy_i(1) = eFrame_Row(2,atom1)
560     uy_i(2) = eFrame_Row(5,atom1)
561     uy_i(3) = eFrame_Row(8,atom1)
562     uz_i(1) = eFrame_Row(3,atom1)
563     uz_i(2) = eFrame_Row(6,atom1)
564     uz_i(3) = eFrame_Row(9,atom1)
565     #else
566     ux_i(1) = eFrame(1,atom1)
567     ux_i(2) = eFrame(4,atom1)
568     ux_i(3) = eFrame(7,atom1)
569     uy_i(1) = eFrame(2,atom1)
570     uy_i(2) = eFrame(5,atom1)
571     uy_i(3) = eFrame(8,atom1)
572     uz_i(1) = eFrame(3,atom1)
573     uz_i(2) = eFrame(6,atom1)
574     uz_i(3) = eFrame(9,atom1)
575     #endif
576     cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
577     cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
578     cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
579     endif
580    
581 gezelter 411 if (j_is_Charge) then
582     q_j = ElectrostaticMap(me2)%charge
583     endif
584 gezelter 507
585 gezelter 411 if (j_is_Dipole) then
586     mu_j = ElectrostaticMap(me2)%dipole_moment
587     #ifdef IS_MPI
588 gezelter 439 uz_j(1) = eFrame_Col(3,atom2)
589     uz_j(2) = eFrame_Col(6,atom2)
590     uz_j(3) = eFrame_Col(9,atom2)
591 gezelter 411 #else
592 gezelter 439 uz_j(1) = eFrame(3,atom2)
593     uz_j(2) = eFrame(6,atom2)
594     uz_j(3) = eFrame(9,atom2)
595 gezelter 411 #endif
596 chrisfen 465 ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
597 gezelter 411
598     if (j_is_SplitDipole) then
599     d_j = ElectrostaticMap(me2)%split_dipole_distance
600     endif
601     endif
602    
603 gezelter 439 if (j_is_Quadrupole) then
604     qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1)
605     qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2)
606     qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3)
607     #ifdef IS_MPI
608     ux_j(1) = eFrame_Col(1,atom2)
609     ux_j(2) = eFrame_Col(4,atom2)
610     ux_j(3) = eFrame_Col(7,atom2)
611     uy_j(1) = eFrame_Col(2,atom2)
612     uy_j(2) = eFrame_Col(5,atom2)
613     uy_j(3) = eFrame_Col(8,atom2)
614     uz_j(1) = eFrame_Col(3,atom2)
615     uz_j(2) = eFrame_Col(6,atom2)
616     uz_j(3) = eFrame_Col(9,atom2)
617     #else
618     ux_j(1) = eFrame(1,atom2)
619     ux_j(2) = eFrame(4,atom2)
620     ux_j(3) = eFrame(7,atom2)
621     uy_j(1) = eFrame(2,atom2)
622     uy_j(2) = eFrame(5,atom2)
623     uy_j(3) = eFrame(8,atom2)
624     uz_j(1) = eFrame(3,atom2)
625     uz_j(2) = eFrame(6,atom2)
626     uz_j(3) = eFrame(9,atom2)
627     #endif
628     cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
629     cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
630     cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
631     endif
632 chrisfen 554
633 gezelter 411 epot = 0.0_dp
634     dudx = 0.0_dp
635     dudy = 0.0_dp
636     dudz = 0.0_dp
637    
638 gezelter 439 dudux_i = 0.0_dp
639     duduy_i = 0.0_dp
640     duduz_i = 0.0_dp
641 gezelter 411
642 gezelter 439 dudux_j = 0.0_dp
643     duduy_j = 0.0_dp
644     duduz_j = 0.0_dp
645 gezelter 411
646     if (i_is_Charge) then
647    
648     if (j_is_Charge) then
649 chrisfen 739 if (screeningMethod .eq. DAMPED) then
650     f0 = derfc(dampingAlpha*rij)
651     varEXP = exp(-alpha2*rij*rij)
652     f1 = alphaPi*rij*varEXP + f0
653     endif
654 gezelter 507
655 chrisfen 739 preVal = pre11 * q_i * q_j
656    
657 chrisfen 710 if (summationMethod .eq. SHIFTED_POTENTIAL) then
658 chrisfen 739 vterm = preVal * (riji*f0 - rcuti*f0c)
659 chrisfen 597
660 chrisfen 739 dudr = -sw * preVal * riji * riji * f1
661    
662 chrisfen 710 elseif (summationMethod .eq. SHIFTED_FORCE) then
663 chrisfen 739 vterm = preVal * ( riji*f0 - rcuti*f0c + &
664 chrisfen 716 f1c*rcuti2*(rij-defaultCutoff) )
665    
666 chrisfen 739 dudr = -sw*preVal * (riji*riji*f1 - rcuti2*f1c)
667    
668 chrisfen 695 elseif (summationMethod .eq. REACTION_FIELD) then
669     rfVal = preRF*rij*rij
670     vterm = preVal * ( riji + rfVal )
671 chrisfen 700
672 chrisfen 695 dudr = sw * preVal * ( 2.0d0*rfVal - riji )*riji
673 chrisfen 739
674 chrisfen 597 else
675 chrisfen 739 vterm = preVal * riji*f0
676 chrisfen 597
677 chrisfen 739 dudr = - sw * preVal * riji*riji*f1
678    
679 chrisfen 597 endif
680    
681 chrisfen 739 vpair = vpair + vterm
682     epot = epot + sw*vterm
683    
684     dudx = dudx + dudr * xhat
685     dudy = dudy + dudr * yhat
686     dudz = dudz + dudr * zhat
687    
688 gezelter 411 endif
689    
690     if (j_is_Dipole) then
691    
692 chrisfen 626 pref = pre12 * q_i * mu_j
693 gezelter 411
694 chrisfen 710 !!$ if (summationMethod .eq. UNDAMPED_WOLF) then
695     !!$ ri2 = riji * riji
696     !!$ ri3 = ri2 * riji
697     !!$
698     !!$ pref = pre12 * q_i * mu_j
699     !!$ vterm = - pref * ct_j * (ri2 - rcuti2)
700     !!$ vpair = vpair + vterm
701     !!$ epot = epot + sw*vterm
702     !!$
703     !!$ !! this has a + sign in the () because the rij vector is
704     !!$ !! r_j - r_i and the charge-dipole potential takes the origin
705     !!$ !! as the point dipole, which is atom j in this case.
706     !!$
707     !!$ dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
708     !!$ - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
709     !!$ dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
710     !!$ - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
711     !!$ dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
712     !!$ - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
713     !!$
714     !!$ duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
715     !!$ duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
716     !!$ duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
717     !!$
718     !!$ elseif (summationMethod .eq. REACTION_FIELD) then
719 gezelter 507
720 chrisfen 710 if (summationMethod .eq. REACTION_FIELD) then
721 chrisfen 700 ri2 = riji * riji
722     ri3 = ri2 * riji
723 chrisfen 696
724     pref = pre12 * q_i * mu_j
725     vterm = - pref * ct_j * ( ri2 - preRF2*rij )
726     vpair = vpair + vterm
727     epot = epot + sw*vterm
728    
729     !! this has a + sign in the () because the rij vector is
730     !! r_j - r_i and the charge-dipole potential takes the origin
731     !! as the point dipole, which is atom j in this case.
732    
733     dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
734     preRF2*uz_j(1) )
735     dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
736     preRF2*uz_j(2) )
737     dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
738     preRF2*uz_j(3) )
739     duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
740     duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
741     duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
742    
743 chrisfen 597 else
744     if (j_is_SplitDipole) then
745     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
746     ri = 1.0_dp / BigR
747     scale = rij * ri
748     else
749     ri = riji
750     scale = 1.0_dp
751     endif
752    
753     ri2 = ri * ri
754     ri3 = ri2 * ri
755     sc2 = scale * scale
756 chrisfen 626
757     pref = pre12 * q_i * mu_j
758 chrisfen 597 vterm = - pref * ct_j * ri2 * scale
759 chrisfen 626 vpair = vpair + vterm
760     epot = epot + sw*vterm
761 chrisfen 597
762     !! this has a + sign in the () because the rij vector is
763     !! r_j - r_i and the charge-dipole potential takes the origin
764     !! as the point dipole, which is atom j in this case.
765    
766 chrisfen 626 dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
767     dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
768     dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
769 chrisfen 597
770 chrisfen 626 duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
771     duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
772     duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
773 gezelter 411
774 chrisfen 597 endif
775 gezelter 411 endif
776 gezelter 421
777 gezelter 439 if (j_is_Quadrupole) then
778     ri2 = riji * riji
779     ri3 = ri2 * riji
780 gezelter 440 ri4 = ri2 * ri2
781 gezelter 439 cx2 = cx_j * cx_j
782     cy2 = cy_j * cy_j
783     cz2 = cz_j * cz_j
784    
785 chrisfen 710 !!$ if (summationMethod .eq. UNDAMPED_WOLF) then
786     !!$ pref = pre14 * q_i / 3.0_dp
787     !!$ vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
788     !!$ qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
789     !!$ qzz_j * (3.0_dp*cz2 - 1.0_dp) )
790     !!$ vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
791     !!$ qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
792     !!$ qzz_j * (3.0_dp*cz2 - 1.0_dp) )
793     !!$ vpair = vpair + ( vterm1 - vterm2 )
794     !!$ epot = epot + sw*( vterm1 - vterm2 )
795     !!$
796     !!$ dudx = dudx - (5.0_dp * &
797     !!$ (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
798     !!$ (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
799     !!$ qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
800     !!$ (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
801     !!$ qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
802     !!$ (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
803     !!$ qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
804     !!$ dudy = dudy - (5.0_dp * &
805     !!$ (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
806     !!$ (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
807     !!$ qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
808     !!$ (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
809     !!$ qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
810     !!$ (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
811     !!$ qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
812     !!$ dudz = dudz - (5.0_dp * &
813     !!$ (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
814     !!$ (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
815     !!$ qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
816     !!$ (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
817     !!$ qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
818     !!$ (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
819     !!$ qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
820     !!$
821     !!$ dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
822     !!$ rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
823     !!$ dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
824     !!$ rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
825     !!$ dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
826     !!$ rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
827     !!$
828     !!$ duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
829     !!$ rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
830     !!$ duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
831     !!$ rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
832     !!$ duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
833     !!$ rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
834     !!$
835     !!$ duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
836     !!$ rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
837     !!$ duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
838     !!$ rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
839     !!$ duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
840     !!$ rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
841     !!$
842     !!$ else
843 chrisfen 626 pref = pre14 * q_i / 3.0_dp
844 chrisfen 597 vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
845     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
846     qzz_j * (3.0_dp*cz2 - 1.0_dp))
847 chrisfen 626 vpair = vpair + vterm
848     epot = epot + sw*vterm
849 chrisfen 597
850 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
851 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
852     qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
853     qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
854 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
855 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
856     qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
857     qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
858 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
859 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
860     qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
861     qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
862    
863 chrisfen 626 dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
864     dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
865     dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
866 chrisfen 597
867 chrisfen 626 duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
868     duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
869     duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
870 chrisfen 597
871 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
872     duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
873     duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
874 chrisfen 597
875 chrisfen 710 !!$ endif
876 gezelter 439 endif
877 gezelter 411 endif
878 gezelter 507
879 gezelter 411 if (i_is_Dipole) then
880 gezelter 507
881 gezelter 411 if (j_is_Charge) then
882 chrisfen 626
883 chrisfen 719 if (summationMethod .eq. SHIFTED_POTENTIAL) then
884 chrisfen 597 ri2 = riji * riji
885     ri3 = ri2 * riji
886 chrisfen 719
887     pref = pre12 * q_j * mu_i
888     pot_term = ri2 - rcuti2
889     vterm = pref * ct_i * pot_term
890     vpair = vpair + vterm
891     epot = epot + sw*vterm
892    
893     dudx = dudx + sw*pref * ( ri3*(uz_i(1)-3.0d0*ct_i*xhat) )
894     dudy = dudy + sw*pref * ( ri3*(uz_i(2)-3.0d0*ct_i*yhat) )
895     dudz = dudz + sw*pref * ( ri3*(uz_i(3)-3.0d0*ct_i*zhat) )
896    
897     duduz_i(1) = duduz_i(1) + sw*pref * xhat * pot_term
898     duduz_i(2) = duduz_i(2) + sw*pref * yhat * pot_term
899     duduz_i(3) = duduz_i(3) + sw*pref * zhat * pot_term
900 gezelter 507
901 chrisfen 719 elseif (summationMethod .eq. SHIFTED_FORCE) then
902     ri2 = riji * riji
903     ri3 = ri2 * riji
904    
905 chrisfen 626 pref = pre12 * q_j * mu_i
906 chrisfen 719 pot_term = ri2 - rcuti2 + 2.0d0*rcuti3*( rij - defaultCutoff )
907     vterm = pref * ct_i * pot_term
908     vpair = vpair + vterm
909     epot = epot + sw*vterm
910    
911     dudx = dudx + sw*pref * ( (ri3-rcuti3)*(uz_i(1)-3.0d0*ct_i*xhat) )
912     dudy = dudy + sw*pref * ( (ri3-rcuti3)*(uz_i(2)-3.0d0*ct_i*yhat) )
913     dudz = dudz + sw*pref * ( (ri3-rcuti3)*(uz_i(3)-3.0d0*ct_i*zhat) )
914    
915     duduz_i(1) = duduz_i(1) + sw*pref * xhat * pot_term
916     duduz_i(2) = duduz_i(2) + sw*pref * yhat * pot_term
917     duduz_i(3) = duduz_i(3) + sw*pref * zhat * pot_term
918    
919     elseif (summationMethod .eq. REACTION_FIELD) then
920     ri2 = riji * riji
921     ri3 = ri2 * riji
922    
923     pref = pre12 * q_j * mu_i
924 chrisfen 700 vterm = pref * ct_i * ( ri2 - preRF2*rij )
925 chrisfen 696 vpair = vpair + vterm
926     epot = epot + sw*vterm
927    
928 chrisfen 700 dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
929     preRF2*uz_i(1) )
930     dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
931     preRF2*uz_i(2) )
932     dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
933     preRF2*uz_i(3) )
934 chrisfen 696
935 chrisfen 700 duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
936     duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
937     duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
938 chrisfen 696
939 chrisfen 597 else
940     if (i_is_SplitDipole) then
941 gezelter 421 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
942     ri = 1.0_dp / BigR
943 chrisfen 597 scale = rij * ri
944     else
945 gezelter 421 ri = riji
946     scale = 1.0_dp
947     endif
948 chrisfen 597
949     ri2 = ri * ri
950     ri3 = ri2 * ri
951     sc2 = scale * scale
952 chrisfen 626
953     pref = pre12 * q_j * mu_i
954 chrisfen 597 vterm = pref * ct_i * ri2 * scale
955 chrisfen 626 vpair = vpair + vterm
956     epot = epot + sw*vterm
957 chrisfen 597
958 chrisfen 626 dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
959     dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
960     dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
961 chrisfen 597
962 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
963     duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
964     duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
965 gezelter 421 endif
966 chrisfen 597 endif
967 chrisfen 626
968 chrisfen 597 if (j_is_Dipole) then
969 chrisfen 719 ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
970    
971     ri2 = riji * riji
972     ri3 = ri2 * riji
973     ri4 = ri2 * ri2
974    
975     pref = pre22 * mu_i * mu_j
976 gezelter 421
977 chrisfen 719 !!$ if (summationMethod .eq. SHIFTED_POTENTIAL) then
978     !!$ a0 = ct_ij - 3.0d0 * ct_i * ct_j
979     !!$ pot_term = ri3 - rcuti3
980     !!$
981     !!$ vterm = pref*pot_term*a0
982     !!$ vpair = vpair + vterm
983     !!$ epot = epot + sw*vterm
984     !!$
985     !!$ a1 = 5.0d0 * ct_i * ct_j - ct_ij
986     !!$
987     !!$ dudx = dudx + sw*pref*3.0d0*ri4 &
988     !!$ * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989     !!$ dudy = dudy + sw*pref*3.0d0*ri4 &
990     !!$ * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991     !!$ dudz = dudz + sw*pref*3.0d0*ri4 &
992     !!$ * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
993     !!$
994     !!$ duduz_i(1) = duduz_i(1) + sw*pref*( pot_term &
995     !!$ * (uz_j(1) - 3.0d0*ct_j*xhat) )
996     !!$ duduz_i(2) = duduz_i(2) + sw*pref*( pot_term &
997     !!$ * (uz_j(2) - 3.0d0*ct_j*yhat) )
998     !!$ duduz_i(3) = duduz_i(3) + sw*pref*( pot_term &
999     !!$ * (uz_j(3) - 3.0d0*ct_j*zhat) )
1000     !!$ duduz_j(1) = duduz_j(1) + sw*pref*( pot_term &
1001     !!$ * (uz_i(1) - 3.0d0*ct_i*xhat) )
1002     !!$ duduz_j(2) = duduz_j(2) + sw*pref*( pot_term &
1003     !!$ * (uz_i(2) - 3.0d0*ct_i*yhat) )
1004     !!$ duduz_j(3) = duduz_j(3) + sw*pref*( pot_term &
1005     !!$ * (uz_i(3) - 3.0d0*ct_i*zhat) )
1006 chrisfen 703 !!$
1007 chrisfen 719 !!$ elseif (summationMethod .eq. SHIFTED_FORCE) then
1008     !!$ a0 = ct_ij - 3.0d0 * ct_i * ct_j
1009     !!$ pot_term = ri3 - rcuti3 + 3.0d0*rcuti4*( rij - defaultCutoff )
1010     !!$
1011     !!$ vterm = pref*pot_term*a0
1012 chrisfen 703 !!$ vpair = vpair + vterm
1013     !!$ epot = epot + sw*vterm
1014     !!$
1015     !!$ a1 = 5.0d0 * ct_i * ct_j - ct_ij
1016     !!$
1017 chrisfen 719 !!$ dudx = dudx + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1018     !!$ * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1019     !!$ dudy = dudy + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1020     !!$ * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1021     !!$ dudz = dudz + sw*pref*3.0d0*( ri4 - rcuti4 ) &
1022     !!$ * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1023 chrisfen 703 !!$
1024 chrisfen 719 !!$ duduz_i(1) = duduz_i(1) + sw*pref*( pot_term &
1025     !!$ * (uz_j(1) - 3.0d0*ct_j*xhat) )
1026     !!$ duduz_i(2) = duduz_i(2) + sw*pref*( pot_term &
1027     !!$ * (uz_j(2) - 3.0d0*ct_j*yhat) )
1028     !!$ duduz_i(3) = duduz_i(3) + sw*pref*( pot_term &
1029     !!$ * (uz_j(3) - 3.0d0*ct_j*zhat) )
1030     !!$ duduz_j(1) = duduz_j(1) + sw*pref*( pot_term &
1031     !!$ * (uz_i(1) - 3.0d0*ct_i*xhat) )
1032     !!$ duduz_j(2) = duduz_j(2) + sw*pref*( pot_term &
1033     !!$ * (uz_i(2) - 3.0d0*ct_i*yhat) )
1034     !!$ duduz_j(3) = duduz_j(3) + sw*pref*( pot_term &
1035     !!$ * (uz_i(3) - 3.0d0*ct_i*zhat) )
1036     !!$
1037 chrisfen 710 !!$ elseif (summationMethod .eq. REACTION_FIELD) then
1038     if (summationMethod .eq. REACTION_FIELD) then
1039 chrisfen 695 vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1040     preRF2*ct_ij )
1041     vpair = vpair + vterm
1042     epot = epot + sw*vterm
1043    
1044     a1 = 5.0d0 * ct_i * ct_j - ct_ij
1045    
1046     dudx = dudx + sw*pref*3.0d0*ri4 &
1047     * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1048     dudy = dudy + sw*pref*3.0d0*ri4 &
1049     * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1050     dudz = dudz + sw*pref*3.0d0*ri4 &
1051     * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1052    
1053     duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1054     - preRF2*uz_j(1))
1055     duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1056     - preRF2*uz_j(2))
1057     duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1058     - preRF2*uz_j(3))
1059     duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1060     - preRF2*uz_i(1))
1061     duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1062     - preRF2*uz_i(2))
1063     duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1064     - preRF2*uz_i(3))
1065    
1066 chrisfen 597 else
1067     if (i_is_SplitDipole) then
1068     if (j_is_SplitDipole) then
1069     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
1070     else
1071     BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
1072     endif
1073     ri = 1.0_dp / BigR
1074     scale = rij * ri
1075     else
1076     if (j_is_SplitDipole) then
1077     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
1078     ri = 1.0_dp / BigR
1079     scale = rij * ri
1080     else
1081     ri = riji
1082     scale = 1.0_dp
1083     endif
1084     endif
1085    
1086     sc2 = scale * scale
1087 chrisfen 719
1088 chrisfen 597 vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1089 chrisfen 626 vpair = vpair + vterm
1090     epot = epot + sw*vterm
1091 chrisfen 597
1092     a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1093    
1094 chrisfen 626 dudx = dudx + sw*pref*3.0d0*ri4*scale &
1095     *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1096     dudy = dudy + sw*pref*3.0d0*ri4*scale &
1097     *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1098     dudz = dudz + sw*pref*3.0d0*ri4*scale &
1099     *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1100 chrisfen 597
1101 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1102     *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1103     duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1104     *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1105     duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1106     *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1107 chrisfen 597
1108 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1109     *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1110     duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1111     *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1112     duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1113     *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1114 chrisfen 597 endif
1115 gezelter 411 endif
1116     endif
1117 gezelter 439
1118     if (i_is_Quadrupole) then
1119     if (j_is_Charge) then
1120 gezelter 507
1121 gezelter 439 ri2 = riji * riji
1122     ri3 = ri2 * riji
1123 gezelter 440 ri4 = ri2 * ri2
1124 gezelter 439 cx2 = cx_i * cx_i
1125     cy2 = cy_i * cy_i
1126     cz2 = cz_i * cz_i
1127 gezelter 507
1128 chrisfen 710 !!$ if (summationMethod .eq. UNDAMPED_WOLF) then
1129     !!$ pref = pre14 * q_j / 3.0_dp
1130     !!$ vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1131     !!$ qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1132     !!$ qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1133     !!$ vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1134     !!$ qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1135     !!$ qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1136     !!$ vpair = vpair + ( vterm1 - vterm2 )
1137     !!$ epot = epot + sw*( vterm1 - vterm2 )
1138     !!$
1139     !!$ dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1140     !!$ sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1141     !!$ qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1142     !!$ (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1143     !!$ qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1144     !!$ (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1145     !!$ qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1146     !!$ dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1147     !!$ sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1148     !!$ qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1149     !!$ (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1150     !!$ qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1151     !!$ (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1152     !!$ qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1153     !!$ dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1154     !!$ sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1155     !!$ qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1156     !!$ (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1157     !!$ qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1158     !!$ (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1159     !!$ qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1160     !!$
1161     !!$ dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1162     !!$ rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1163     !!$ dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1164     !!$ rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1165     !!$ dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1166     !!$ rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1167     !!$
1168     !!$ duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1169     !!$ rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1170     !!$ duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1171     !!$ rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1172     !!$ duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1173     !!$ rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1174     !!$
1175     !!$ duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1176     !!$ rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1177     !!$ duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1178     !!$ rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1179     !!$ duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1180     !!$ rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1181     !!$
1182     !!$ else
1183 chrisfen 626 pref = pre14 * q_j / 3.0_dp
1184 chrisfen 597 vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1185     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1186     qzz_i * (3.0_dp*cz2 - 1.0_dp))
1187 chrisfen 626 vpair = vpair + vterm
1188     epot = epot + sw*vterm
1189 chrisfen 597
1190 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1191 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1192     qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1193     qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1194 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1195 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1196     qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1197     qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1198 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1199 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1200     qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1201     qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1202    
1203 chrisfen 626 dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1204     dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1205     dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1206 chrisfen 597
1207 chrisfen 626 duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1208     duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1209     duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1210 chrisfen 597
1211 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1212     duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1213     duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1214 chrisfen 710 !!$ endif
1215 gezelter 439 endif
1216     endif
1217 gezelter 507
1218    
1219 gezelter 411 if (do_pot) then
1220     #ifdef IS_MPI
1221 chuckv 656 pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1222     pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1223 gezelter 411 #else
1224     pot = pot + epot
1225     #endif
1226     endif
1227 gezelter 507
1228 gezelter 411 #ifdef IS_MPI
1229     f_Row(1,atom1) = f_Row(1,atom1) + dudx
1230     f_Row(2,atom1) = f_Row(2,atom1) + dudy
1231     f_Row(3,atom1) = f_Row(3,atom1) + dudz
1232 gezelter 507
1233 gezelter 411 f_Col(1,atom2) = f_Col(1,atom2) - dudx
1234     f_Col(2,atom2) = f_Col(2,atom2) - dudy
1235     f_Col(3,atom2) = f_Col(3,atom2) - dudz
1236 gezelter 507
1237 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1238 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1239     t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1240     t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1241 gezelter 411 endif
1242 gezelter 439 if (i_is_Quadrupole) then
1243     t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1244     t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1245     t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1246 gezelter 411
1247 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1248     t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1249     t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1250     endif
1251    
1252 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1253 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1254     t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1255     t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1256 gezelter 411 endif
1257 gezelter 439 if (j_is_Quadrupole) then
1258     t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1259     t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1260     t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1261 gezelter 411
1262 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1263     t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1264     t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1265     endif
1266    
1267 gezelter 411 #else
1268     f(1,atom1) = f(1,atom1) + dudx
1269     f(2,atom1) = f(2,atom1) + dudy
1270     f(3,atom1) = f(3,atom1) + dudz
1271 gezelter 507
1272 gezelter 411 f(1,atom2) = f(1,atom2) - dudx
1273     f(2,atom2) = f(2,atom2) - dudy
1274     f(3,atom2) = f(3,atom2) - dudz
1275 gezelter 507
1276 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1277 gezelter 439 t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1278     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1279     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1280 gezelter 411 endif
1281 gezelter 439 if (i_is_Quadrupole) then
1282     t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1283     t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1284     t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1285    
1286     t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1287     t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1288     t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1289     endif
1290    
1291 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1292 gezelter 439 t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1293     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1294     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1295 gezelter 411 endif
1296 gezelter 439 if (j_is_Quadrupole) then
1297     t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1298     t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1299     t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1300    
1301     t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1302     t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1303     t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1304     endif
1305    
1306 gezelter 411 #endif
1307 gezelter 507
1308 gezelter 411 #ifdef IS_MPI
1309     id1 = AtomRowToGlobal(atom1)
1310     id2 = AtomColToGlobal(atom2)
1311     #else
1312     id1 = atom1
1313     id2 = atom2
1314     #endif
1315    
1316     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1317 gezelter 507
1318 gezelter 411 fpair(1) = fpair(1) + dudx
1319     fpair(2) = fpair(2) + dudy
1320     fpair(3) = fpair(3) + dudz
1321    
1322     endif
1323    
1324     return
1325     end subroutine doElectrostaticPair
1326 chuckv 492
1327     subroutine destroyElectrostaticTypes()
1328    
1329 gezelter 507 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1330    
1331 chuckv 492 end subroutine destroyElectrostaticTypes
1332    
1333 chrisfen 703 subroutine self_self(atom1, eFrame, mypot, t, do_pot)
1334 chrisfen 695 logical, intent(in) :: do_pot
1335 chrisfen 682 integer, intent(in) :: atom1
1336 chrisfen 695 integer :: atid1
1337 chrisfen 682 real(kind=dp), dimension(9,nLocal) :: eFrame
1338 chrisfen 695 real(kind=dp), dimension(3,nLocal) :: t
1339 chrisfen 703 real(kind=dp) :: mu1, c1
1340     real(kind=dp) :: preVal, epot, mypot
1341 chrisfen 695 real(kind=dp) :: eix, eiy, eiz
1342 chrisfen 682
1343 chrisfen 695 ! this is a local only array, so we use the local atom type id's:
1344     atid1 = atid(atom1)
1345 chrisfen 703
1346     if (.not.summationMethodChecked) then
1347     call checkSummationMethod()
1348     endif
1349 chrisfen 695
1350 chrisfen 703 if (summationMethod .eq. REACTION_FIELD) then
1351     if (ElectrostaticMap(atid1)%is_Dipole) then
1352     mu1 = getDipoleMoment(atid1)
1353    
1354     preVal = pre22 * preRF2 * mu1*mu1
1355     mypot = mypot - 0.5d0*preVal
1356    
1357     ! The self-correction term adds into the reaction field vector
1358    
1359     eix = preVal * eFrame(3,atom1)
1360     eiy = preVal * eFrame(6,atom1)
1361     eiz = preVal * eFrame(9,atom1)
1362    
1363     ! once again, this is self-self, so only the local arrays are needed
1364     ! even for MPI jobs:
1365    
1366     t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1367     eFrame(9,atom1)*eiy
1368     t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1369     eFrame(3,atom1)*eiz
1370     t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1371     eFrame(6,atom1)*eix
1372    
1373     endif
1374    
1375 chrisfen 717 elseif (summationMethod .eq. SHIFTED_FORCE) then
1376     if (ElectrostaticMap(atid1)%is_Charge) then
1377     c1 = getCharge(atid1)
1378    
1379     if (screeningMethod .eq. DAMPED) then
1380     mypot = mypot - (f0c * rcuti * 0.5_dp + &
1381     dampingAlpha*invRootPi) * c1 * c1
1382    
1383     else
1384     mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1385    
1386     endif
1387     endif
1388     endif
1389 chrisfen 695
1390 chrisfen 682 return
1391 chrisfen 703 end subroutine self_self
1392 chrisfen 682
1393 chrisfen 703 subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, myPot, &
1394 chrisfen 700 f, t, do_pot)
1395     logical, intent(in) :: do_pot
1396     integer, intent(in) :: atom1
1397     integer, intent(in) :: atom2
1398     logical :: i_is_Charge, j_is_Charge
1399     logical :: i_is_Dipole, j_is_Dipole
1400     integer :: atid1
1401     integer :: atid2
1402     real(kind=dp), intent(in) :: rij
1403     real(kind=dp), intent(in) :: sw
1404     real(kind=dp), intent(in), dimension(3) :: d
1405     real(kind=dp), intent(inout) :: vpair
1406     real(kind=dp), dimension(9,nLocal) :: eFrame
1407     real(kind=dp), dimension(3,nLocal) :: f
1408     real(kind=dp), dimension(3,nLocal) :: t
1409     real (kind = dp), dimension(3) :: duduz_i
1410     real (kind = dp), dimension(3) :: duduz_j
1411     real (kind = dp), dimension(3) :: uz_i
1412     real (kind = dp), dimension(3) :: uz_j
1413     real(kind=dp) :: q_i, q_j, mu_i, mu_j
1414     real(kind=dp) :: xhat, yhat, zhat
1415     real(kind=dp) :: ct_i, ct_j
1416     real(kind=dp) :: ri2, ri3, riji, vterm
1417 chrisfen 703 real(kind=dp) :: pref, preVal, rfVal, myPot
1418 chrisfen 700 real(kind=dp) :: dudx, dudy, dudz, dudr
1419    
1420 chrisfen 703 if (.not.summationMethodChecked) then
1421     call checkSummationMethod()
1422 chrisfen 700 endif
1423    
1424     dudx = 0.0d0
1425     dudy = 0.0d0
1426     dudz = 0.0d0
1427    
1428     riji = 1.0d0/rij
1429    
1430     xhat = d(1) * riji
1431     yhat = d(2) * riji
1432     zhat = d(3) * riji
1433    
1434     ! this is a local only array, so we use the local atom type id's:
1435     atid1 = atid(atom1)
1436     atid2 = atid(atom2)
1437     i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1438     j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1439     i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1440     j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1441    
1442     if (i_is_Charge.and.j_is_Charge) then
1443     q_i = ElectrostaticMap(atid1)%charge
1444     q_j = ElectrostaticMap(atid2)%charge
1445    
1446     preVal = pre11 * q_i * q_j
1447     rfVal = preRF*rij*rij
1448     vterm = preVal * rfVal
1449    
1450 chrisfen 703 myPot = myPot + sw*vterm
1451    
1452 chrisfen 700 dudr = sw*preVal * 2.0d0*rfVal*riji
1453 chrisfen 703
1454 chrisfen 700 dudx = dudx + dudr * xhat
1455     dudy = dudy + dudr * yhat
1456     dudz = dudz + dudr * zhat
1457 chrisfen 703
1458 chrisfen 700 elseif (i_is_Charge.and.j_is_Dipole) then
1459     q_i = ElectrostaticMap(atid1)%charge
1460     mu_j = ElectrostaticMap(atid2)%dipole_moment
1461     uz_j(1) = eFrame(3,atom2)
1462     uz_j(2) = eFrame(6,atom2)
1463     uz_j(3) = eFrame(9,atom2)
1464     ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1465 chrisfen 703
1466 chrisfen 700 ri2 = riji * riji
1467     ri3 = ri2 * riji
1468    
1469     pref = pre12 * q_i * mu_j
1470     vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1471 chrisfen 703 myPot = myPot + sw*vterm
1472    
1473     dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1474     - preRF2*uz_j(1) )
1475     dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1476     - preRF2*uz_j(2) )
1477     dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1478     - preRF2*uz_j(3) )
1479    
1480 chrisfen 700 duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1481     duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1482     duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1483 chrisfen 703
1484 chrisfen 700 elseif (i_is_Dipole.and.j_is_Charge) then
1485     mu_i = ElectrostaticMap(atid1)%dipole_moment
1486     q_j = ElectrostaticMap(atid2)%charge
1487     uz_i(1) = eFrame(3,atom1)
1488     uz_i(2) = eFrame(6,atom1)
1489     uz_i(3) = eFrame(9,atom1)
1490     ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1491 chrisfen 703
1492 chrisfen 700 ri2 = riji * riji
1493     ri3 = ri2 * riji
1494    
1495     pref = pre12 * q_j * mu_i
1496     vterm = pref * ct_i * ( ri2 - preRF2*rij )
1497 chrisfen 703 myPot = myPot + sw*vterm
1498 chrisfen 700
1499 chrisfen 703 dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1500     - preRF2*uz_i(1) )
1501     dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1502     - preRF2*uz_i(2) )
1503     dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1504     - preRF2*uz_i(3) )
1505 chrisfen 700
1506     duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1507     duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1508     duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1509    
1510     endif
1511 chrisfen 703
1512    
1513     ! accumulate the forces and torques resulting from the self term
1514 chrisfen 700 f(1,atom1) = f(1,atom1) + dudx
1515     f(2,atom1) = f(2,atom1) + dudy
1516     f(3,atom1) = f(3,atom1) + dudz
1517    
1518     f(1,atom2) = f(1,atom2) - dudx
1519     f(2,atom2) = f(2,atom2) - dudy
1520     f(3,atom2) = f(3,atom2) - dudz
1521    
1522     if (i_is_Dipole) then
1523     t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1524     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1525     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1526     elseif (j_is_Dipole) then
1527     t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1528     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1529     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1530     endif
1531    
1532     return
1533     end subroutine rf_self_excludes
1534    
1535 gezelter 411 end module electrostatic_module