ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 644
Committed: Tue Oct 4 19:33:22 2005 UTC (19 years, 7 months ago) by chrisfen
File size: 46599 byte(s)
Log Message: 
maybe some work on wolf

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