ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 700
Committed: Wed Oct 26 23:31:40 2005 UTC (19 years, 7 months ago) by chrisfen
File size: 57031 byte(s)
Log Message: 
reaction field looks to be working now - for charges and dipoles alike

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