ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 2399
Committed: Wed Oct 26 23:31:40 2005 UTC (18 years, 10 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 2095 !!
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 2204
44 gezelter 2095 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 2355
58 gezelter 2301 #define __FORTRAN90
59 chuckv 2355 #include "UseTheForce/DarkSide/fInteractionMap.h"
60 gezelter 2301 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61    
62 chuckv 2355
63 gezelter 2118 !! 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 2095 real(kind=dp), parameter :: pre11 = 332.0637778_dp
67 gezelter 2118 !! 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 2095
77 gezelter 2301 !! variables to handle different summation methods for long-range electrostatics:
78     integer, save :: summationMethod = NONE
79 chrisfen 2302 logical, save :: summationMethodChecked = .false.
80 gezelter 2301 real(kind=DP), save :: defaultCutoff = 0.0_DP
81 chrisfen 2381 real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82 gezelter 2301 logical, save :: haveDefaultCutoff = .false.
83     real(kind=DP), save :: dampingAlpha = 0.0_DP
84     logical, save :: haveDampingAlpha = .false.
85 chrisfen 2381 real(kind=DP), save :: dielectric = 1.0_DP
86 gezelter 2301 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 2381 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 2394 real(kind=dp), save :: preRF2 = 0.0_DP
101 chrisfen 2381 logical, save :: preRFCalculated = .false.
102    
103 chuckv 2331 #ifdef __IFC
104     ! error function for ifc version > 7.
105 chuckv 2330 double precision, external :: derfc
106 chuckv 2331 #endif
107    
108 gezelter 2301 public :: setElectrostaticSummationMethod
109     public :: setElectrostaticCutoffRadius
110     public :: setDampedWolfAlpha
111     public :: setReactionFieldDielectric
112 chrisfen 2381 public :: setReactionFieldPrefactor
113 gezelter 2095 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 2189 public :: destroyElectrostaticTypes
122 chrisfen 2394 public :: rf_self_self
123 chrisfen 2399 public :: rf_self_excludes
124 gezelter 2095
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 2229 logical :: is_Tap = .false.
132 gezelter 2095 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 2301 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 2309 summationMethod = the_ESM
150 chrisfen 2325
151 gezelter 2301 end subroutine setElectrostaticSummationMethod
152    
153 chrisfen 2381 subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
154 gezelter 2301 real(kind=dp), intent(in) :: thisRcut
155 chrisfen 2381 real(kind=dp), intent(in) :: thisRsw
156 gezelter 2301 defaultCutoff = thisRcut
157 chrisfen 2381 rrf = defaultCutoff
158     rt = thisRsw
159 gezelter 2301 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 2381 subroutine setReactionFieldPrefactor
175     if (haveDefaultCutoff .and. haveDielectric) then
176     defaultCutoff2 = defaultCutoff*defaultCutoff
177 chrisfen 2394 preRF = (dielectric-1.0d0) / &
178 chrisfen 2381 ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
179 chrisfen 2394 preRF2 = 2.0d0*preRF
180 chrisfen 2381 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 2095 subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
189 chrisfen 2229 is_SplitDipole, is_Quadrupole, is_Tap, status)
190 gezelter 2204
191 gezelter 2095 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 2229 logical, intent(in) :: is_Tap
197 gezelter 2095 integer, intent(out) :: status
198     integer :: nAtypes, myATID, i, j
199    
200     status = 0
201     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
202 gezelter 2204
203 gezelter 2095 !! 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 2204
208 gezelter 2095 nAtypes = getSize(atypes)
209 gezelter 2204
210 gezelter 2095 if (nAtypes == 0) then
211     status = -1
212     return
213     end if
214 gezelter 2204
215 gezelter 2095 if (.not. allocated(ElectrostaticMap)) then
216     allocate(ElectrostaticMap(nAtypes))
217     endif
218 gezelter 2204
219 gezelter 2095 end if
220    
221     if (myATID .gt. size(ElectrostaticMap)) then
222     status = -1
223     return
224     endif
225 gezelter 2204
226 gezelter 2095 ! 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 2229 ElectrostaticMap(myATID)%is_Tap = is_Tap
234 gezelter 2204
235 gezelter 2095 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 2204 endif
263 gezelter 2095
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 2204
354 gezelter 2095 do i = 1, 3
355 gezelter 2204 ElectrostaticMap(myATID)%quadrupole_moments(i) = &
356     quadrupole_moments(i)
357     enddo
358 gezelter 2095
359     end subroutine setQuadrupoleMoments
360    
361 gezelter 2204
362 gezelter 2095 function getCharge(atid) result (c)
363     integer, intent(in) :: atid
364     integer :: localError
365     real(kind=dp) :: c
366 gezelter 2204
367 gezelter 2095 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 2204
372 gezelter 2095 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 2204
377 gezelter 2095 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 2204
385 gezelter 2095 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 2204
390 gezelter 2095 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 2204
395 gezelter 2095 dm = ElectrostaticMap(atid)%dipole_moment
396     end function getDipoleMoment
397    
398 gezelter 2301 subroutine checkSummationMethod()
399    
400 chrisfen 2306 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 2301 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 2302 if (.not.haveDefaultCutoff) then
417     call handleError("checkSummationMethod", "no Default Cutoff set!")
418     endif
419    
420     constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
421 chuckv 2330 constERFC = derfc(dampingAlpha*defaultCutoff)
422 chrisfen 2339 invRootPi = 0.56418958354775628695d0
423     alphaPi = 2*dampingAlpha*invRootPi
424 chrisfen 2343
425 gezelter 2301 haveDWAconstants = .true.
426     endif
427     endif
428    
429 chrisfen 2302 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 2301 end subroutine checkSummationMethod
437    
438    
439    
440 gezelter 2095 subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
441 chrisfen 2399 vpair, fpair, pot, eFrame, f, t, do_pot, indirect_only)
442 gezelter 2204
443 chrisfen 2399 logical, intent(in) :: do_pot, indirect_only
444 gezelter 2204
445 gezelter 2095 integer, intent(in) :: atom1, atom2
446     integer :: localError
447    
448 chrisfen 2306 real(kind=dp), intent(in) :: rij, r2, sw
449 gezelter 2095 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 2325 real( kind = dp ) :: pot
454 gezelter 2095 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 2204
458 gezelter 2123 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 2095
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 2229 logical :: i_is_Tap, j_is_Tap
466 gezelter 2095 integer :: me1, me2, id1, id2
467     real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
468 gezelter 2123 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 2095 real (kind=dp) :: ct_i, ct_j, ct_ij, a1
474 gezelter 2105 real (kind=dp) :: riji, ri, ri2, ri3, ri4
475 chrisfen 2296 real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
476 gezelter 2105 real (kind=dp) :: xhat, yhat, zhat
477 gezelter 2095 real (kind=dp) :: dudx, dudy, dudz
478 chrisfen 2325 real (kind=dp) :: scale, sc2, bigR
479 chrisfen 2339 real (kind=dp) :: varERFC, varEXP
480 chrisfen 2343 real (kind=dp) :: limScale
481 chrisfen 2394 real (kind=dp) :: preVal, rfVal
482 gezelter 2095
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 2301 if (.not.summationMethodChecked) then
489     call checkSummationMethod()
490     endif
491    
492 chrisfen 2394 if (.not.preRFCalculated) then
493     call setReactionFieldPrefactor()
494     endif
495 gezelter 2301
496 gezelter 2095 #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 2343
508 gezelter 2105 xhat = d(1) * riji
509     yhat = d(2) * riji
510     zhat = d(3) * riji
511 gezelter 2095
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 2229 i_is_Tap = ElectrostaticMap(me1)%is_Tap
518 gezelter 2095
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 2229 j_is_Tap = ElectrostaticMap(me2)%is_Tap
524 gezelter 2095
525     if (i_is_Charge) then
526     q_i = ElectrostaticMap(me1)%charge
527     endif
528 gezelter 2204
529 gezelter 2095 if (i_is_Dipole) then
530     mu_i = ElectrostaticMap(me1)%dipole_moment
531     #ifdef IS_MPI
532 gezelter 2123 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 2095 #else
536 gezelter 2123 uz_i(1) = eFrame(3,atom1)
537     uz_i(2) = eFrame(6,atom1)
538     uz_i(3) = eFrame(9,atom1)
539 gezelter 2095 #endif
540 gezelter 2123 ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
541 gezelter 2095
542     if (i_is_SplitDipole) then
543     d_i = ElectrostaticMap(me1)%split_dipole_distance
544     endif
545 gezelter 2204
546 gezelter 2095 endif
547    
548 gezelter 2123 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 2095 if (j_is_Charge) then
579     q_j = ElectrostaticMap(me2)%charge
580     endif
581 gezelter 2204
582 gezelter 2095 if (j_is_Dipole) then
583     mu_j = ElectrostaticMap(me2)%dipole_moment
584     #ifdef IS_MPI
585 gezelter 2123 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 2095 #else
589 gezelter 2123 uz_j(1) = eFrame(3,atom2)
590     uz_j(2) = eFrame(6,atom2)
591     uz_j(3) = eFrame(9,atom2)
592 gezelter 2095 #endif
593 chrisfen 2162 ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
594 gezelter 2095
595     if (j_is_SplitDipole) then
596     d_j = ElectrostaticMap(me2)%split_dipole_distance
597     endif
598     endif
599    
600 gezelter 2123 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 2251
630 gezelter 2095 epot = 0.0_dp
631     dudx = 0.0_dp
632     dudy = 0.0_dp
633     dudz = 0.0_dp
634    
635 gezelter 2123 dudux_i = 0.0_dp
636     duduy_i = 0.0_dp
637     duduz_i = 0.0_dp
638 gezelter 2095
639 gezelter 2123 dudux_j = 0.0_dp
640     duduy_j = 0.0_dp
641     duduz_j = 0.0_dp
642 gezelter 2095
643     if (i_is_Charge) then
644    
645     if (j_is_Charge) then
646 gezelter 2204
647 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
648 chrisfen 2296 vterm = pre11 * q_i * q_j * (riji - rcuti)
649     vpair = vpair + vterm
650 chrisfen 2325 epot = epot + sw*vterm
651 chrisfen 2296
652 chrisfen 2343 dudr = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
653 chrisfen 2296
654     dudx = dudx + dudr * d(1)
655     dudy = dudy + dudr * d(2)
656     dudz = dudz + dudr * d(3)
657 gezelter 2095
658 chrisfen 2339 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 2343 dudr = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
666     + alphaPi*varEXP) &
667     - (constERFC*rcuti2 &
668     + alphaPi*constEXP)) )
669 chrisfen 2339
670     dudx = dudx + dudr * d(1)
671     dudy = dudy + dudr * d(2)
672     dudz = dudz + dudr * d(3)
673    
674 chrisfen 2394 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 2399
679 chrisfen 2394 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 2296 else
689     vterm = pre11 * q_i * q_j * riji
690     vpair = vpair + vterm
691 chrisfen 2325 epot = epot + sw*vterm
692 chrisfen 2296
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 2095 endif
702    
703     if (j_is_Dipole) then
704    
705 chrisfen 2325 pref = pre12 * q_i * mu_j
706 gezelter 2095
707 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
708 chrisfen 2296 ri2 = riji * riji
709     ri3 = ri2 * riji
710 gezelter 2204
711 chrisfen 2325 pref = pre12 * q_i * mu_j
712 chrisfen 2296 vterm = - pref * ct_j * (ri2 - rcuti2)
713 chrisfen 2325 vpair = vpair + vterm
714     epot = epot + sw*vterm
715 chrisfen 2296
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 2325 dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
721 chrisfen 2296 - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
722 chrisfen 2325 dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
723 chrisfen 2296 - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
724 chrisfen 2325 dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
725 chrisfen 2296 - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
726    
727 chrisfen 2325 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 2095
731 chrisfen 2395 elseif (summationMethod .eq. REACTION_FIELD) then
732 chrisfen 2399 ri2 = riji * riji
733     ri3 = ri2 * riji
734 chrisfen 2395
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 2296 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 2325
768     pref = pre12 * q_i * mu_j
769 chrisfen 2296 vterm = - pref * ct_j * ri2 * scale
770 chrisfen 2325 vpair = vpair + vterm
771     epot = epot + sw*vterm
772 chrisfen 2296
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 2325 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 2296
781 chrisfen 2325 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 2095
785 chrisfen 2296 endif
786 gezelter 2095 endif
787 gezelter 2105
788 gezelter 2123 if (j_is_Quadrupole) then
789     ri2 = riji * riji
790     ri3 = ri2 * riji
791 gezelter 2124 ri4 = ri2 * ri2
792 gezelter 2123 cx2 = cx_j * cx_j
793     cy2 = cy_j * cy_j
794     cz2 = cz_j * cz_j
795    
796 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
797 chrisfen 2325 pref = pre14 * q_i / 3.0_dp
798 chrisfen 2296 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 2325 vpair = vpair + ( vterm1 - vterm2 )
805     epot = epot + sw*( vterm1 - vterm2 )
806 chrisfen 2296
807     dudx = dudx - (5.0_dp * &
808 chrisfen 2325 (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
809 chrisfen 2296 (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 2325 (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
817 chrisfen 2296 (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 2325 (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
825 chrisfen 2296 (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 2325 dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
833 chrisfen 2296 rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
834 chrisfen 2325 dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
835 chrisfen 2296 rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
836 chrisfen 2325 dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
837 chrisfen 2296 rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
838    
839 chrisfen 2325 duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
840 chrisfen 2296 rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
841 chrisfen 2325 duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
842 chrisfen 2296 rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
843 chrisfen 2325 duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
844 chrisfen 2296 rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
845    
846 chrisfen 2325 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
847 chrisfen 2296 rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
848 chrisfen 2325 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
849 chrisfen 2296 rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
850 chrisfen 2325 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
851 chrisfen 2296 rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
852    
853     else
854 chrisfen 2325 pref = pre14 * q_i / 3.0_dp
855 chrisfen 2296 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 2325 vpair = vpair + vterm
859     epot = epot + sw*vterm
860 chrisfen 2296
861 chrisfen 2325 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
862 chrisfen 2296 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 2325 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
866 chrisfen 2296 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 2325 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
870 chrisfen 2296 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 2325 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 2296
878 chrisfen 2325 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 2296
882 chrisfen 2325 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 2296
886     endif
887 gezelter 2123 endif
888 gezelter 2095 endif
889 gezelter 2204
890 gezelter 2095 if (i_is_Dipole) then
891 gezelter 2204
892 gezelter 2095 if (j_is_Charge) then
893 chrisfen 2325
894     pref = pre12 * q_j * mu_i
895    
896 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
897 chrisfen 2296 ri2 = riji * riji
898     ri3 = ri2 * riji
899 gezelter 2204
900 chrisfen 2325 pref = pre12 * q_j * mu_i
901 chrisfen 2296 vterm = pref * ct_i * (ri2 - rcuti2)
902 chrisfen 2325 vpair = vpair + vterm
903     epot = epot + sw*vterm
904 chrisfen 2296
905 chrisfen 2325 dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
906 chrisfen 2296 - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
907 chrisfen 2325 dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
908 chrisfen 2296 - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
909 chrisfen 2325 dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
910 chrisfen 2296 - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
911    
912 chrisfen 2399 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 2095
916 chrisfen 2395 elseif (summationMethod .eq. REACTION_FIELD) then
917 chrisfen 2399 ri2 = riji * riji
918     ri3 = ri2 * riji
919 chrisfen 2395
920     pref = pre12 * q_j * mu_i
921 chrisfen 2399 vterm = pref * ct_i * ( ri2 - preRF2*rij )
922 chrisfen 2395 vpair = vpair + vterm
923     epot = epot + sw*vterm
924    
925 chrisfen 2399 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 2395
932 chrisfen 2399 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 2395
936 chrisfen 2296 else
937     if (i_is_SplitDipole) then
938 gezelter 2105 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
939     ri = 1.0_dp / BigR
940 chrisfen 2296 scale = rij * ri
941     else
942 gezelter 2105 ri = riji
943     scale = 1.0_dp
944     endif
945 chrisfen 2296
946     ri2 = ri * ri
947     ri3 = ri2 * ri
948     sc2 = scale * scale
949 chrisfen 2325
950     pref = pre12 * q_j * mu_i
951 chrisfen 2296 vterm = pref * ct_i * ri2 * scale
952 chrisfen 2325 vpair = vpair + vterm
953     epot = epot + sw*vterm
954 chrisfen 2296
955 chrisfen 2325 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 2296
959 chrisfen 2325 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 2105 endif
963 chrisfen 2296 endif
964 chrisfen 2325
965 chrisfen 2296 if (j_is_Dipole) then
966 gezelter 2105
967 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
968 chrisfen 2296 ri2 = riji * riji
969     ri3 = ri2 * riji
970     ri4 = ri2 * ri2
971 gezelter 2204
972 chrisfen 2325 pref = pre22 * mu_i * mu_j
973 chrisfen 2296 vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
974 chrisfen 2325 vpair = vpair + vterm
975     epot = epot + sw*vterm
976 chrisfen 2296
977     a1 = 5.0d0 * ct_i * ct_j - ct_ij
978    
979 chrisfen 2325 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 2296
992 chrisfen 2325 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
993 chrisfen 2296 - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
994 chrisfen 2325 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
995 chrisfen 2296 - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
996 chrisfen 2325 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
997 chrisfen 2296 - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
998 chrisfen 2325 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
999 chrisfen 2296 - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
1000 chrisfen 2325 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1001 chrisfen 2296 - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
1002 chrisfen 2325 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1003 chrisfen 2296 - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
1004 chrisfen 2325
1005 chrisfen 2394 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 2296 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 2325 pref = pre22 * mu_i * mu_j
1069 chrisfen 2296 vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1070 chrisfen 2325 vpair = vpair + vterm
1071     epot = epot + sw*vterm
1072 chrisfen 2296
1073     a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1074    
1075 chrisfen 2325 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 2296
1082 chrisfen 2325 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 2296
1089 chrisfen 2325 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 2296 endif
1096 gezelter 2095 endif
1097     endif
1098 gezelter 2123
1099     if (i_is_Quadrupole) then
1100     if (j_is_Charge) then
1101 gezelter 2204
1102 gezelter 2123 ri2 = riji * riji
1103     ri3 = ri2 * riji
1104 gezelter 2124 ri4 = ri2 * ri2
1105 gezelter 2123 cx2 = cx_i * cx_i
1106     cy2 = cy_i * cy_i
1107     cz2 = cz_i * cz_i
1108 gezelter 2204
1109 chrisfen 2310 if (summationMethod .eq. UNDAMPED_WOLF) then
1110 chrisfen 2325 pref = pre14 * q_j / 3.0_dp
1111 chrisfen 2296 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 2325 vpair = vpair + ( vterm1 - vterm2 )
1118     epot = epot + sw*( vterm1 - vterm2 )
1119 chrisfen 2296
1120 chrisfen 2325 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 2296 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 2325 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 2296 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 2325 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 2296 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 2325 dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1143 chrisfen 2296 rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1144 chrisfen 2325 dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1145 chrisfen 2296 rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1146 chrisfen 2325 dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1147 chrisfen 2296 rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1148    
1149 chrisfen 2325 duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1150 chrisfen 2296 rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1151 chrisfen 2325 duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1152 chrisfen 2296 rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1153 chrisfen 2325 duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1154 chrisfen 2296 rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1155    
1156 chrisfen 2325 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1157 chrisfen 2296 rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1158 chrisfen 2325 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1159 chrisfen 2296 rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1160 chrisfen 2325 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1161 chrisfen 2296 rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1162 gezelter 2204
1163 chrisfen 2296 else
1164 chrisfen 2325 pref = pre14 * q_j / 3.0_dp
1165 chrisfen 2296 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 2325 vpair = vpair + vterm
1169     epot = epot + sw*vterm
1170 chrisfen 2296
1171 chrisfen 2325 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1172 chrisfen 2296 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 2325 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1176 chrisfen 2296 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 2325 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1180 chrisfen 2296 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 2325 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 2296
1188 chrisfen 2325 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 2296
1192 chrisfen 2325 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 2296 endif
1196 gezelter 2123 endif
1197     endif
1198 gezelter 2204
1199    
1200 gezelter 2095 if (do_pot) then
1201     #ifdef IS_MPI
1202 chuckv 2355 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 2095 #else
1205     pot = pot + epot
1206     #endif
1207     endif
1208 gezelter 2204
1209 gezelter 2095 #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 2204
1214 gezelter 2095 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 2204
1218 gezelter 2095 if (i_is_Dipole .or. i_is_Quadrupole) then
1219 gezelter 2123 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 2095 endif
1223 gezelter 2123 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 2095
1228 gezelter 2123 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 2095 if (j_is_Dipole .or. j_is_Quadrupole) then
1234 gezelter 2123 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 2095 endif
1238 gezelter 2123 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 2095
1243 gezelter 2123 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 2095 #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 2204
1253 gezelter 2095 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 2204
1257 gezelter 2095 if (i_is_Dipole .or. i_is_Quadrupole) then
1258 gezelter 2123 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 2095 endif
1262 gezelter 2123 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 2095 if (j_is_Dipole .or. j_is_Quadrupole) then
1273 gezelter 2123 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 2095 endif
1277 gezelter 2123 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 2095 #endif
1288 gezelter 2204
1289 gezelter 2095 #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 2204
1299 gezelter 2095 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 2189
1308     subroutine destroyElectrostaticTypes()
1309    
1310 gezelter 2204 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1311    
1312 chuckv 2189 end subroutine destroyElectrostaticTypes
1313    
1314 chrisfen 2394 subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1315     logical, intent(in) :: do_pot
1316 chrisfen 2381 integer, intent(in) :: atom1
1317 chrisfen 2394 integer :: atid1
1318 chrisfen 2381 real(kind=dp), dimension(9,nLocal) :: eFrame
1319 chrisfen 2394 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 2381
1324 chrisfen 2399 if (.not.preRFCalculated) then
1325     call setReactionFieldPrefactor()
1326     endif
1327    
1328 chrisfen 2394 ! 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 2399
1337 chrisfen 2394 ! 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 2399
1343 chrisfen 2394 ! 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 2399
1353 chrisfen 2381 endif
1354 chrisfen 2394
1355 chrisfen 2381 return
1356 chrisfen 2394 end subroutine rf_self_self
1357 chrisfen 2381
1358 chrisfen 2399 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 2095 end module electrostatic_module