ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/electrostatic.F90
Revision: 691
Committed: Wed Oct 19 19:24:40 2005 UTC (19 years, 7 months ago) by chrisfen
File size: 53027 byte(s)
Log Message: 
Still had some globals toUpper problems - these changes should fix those...

File Contents

# User Rev Content
1 gezelter 411 !!
2     !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     !!
4     !! The University of Notre Dame grants you ("Licensee") a
5     !! non-exclusive, royalty free, license to use, modify and
6     !! redistribute this software in source and binary code form, provided
7     !! that the following conditions are met:
8     !!
9     !! 1. Acknowledgement of the program authors must be made in any
10     !! publication of scientific results based in part on use of the
11     !! program. An acceptable form of acknowledgement is citation of
12     !! the article in which the program was described (Matthew
13     !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     !! Parallel Simulation Engine for Molecular Dynamics,"
16     !! J. Comput. Chem. 26, pp. 252-271 (2005))
17     !!
18     !! 2. Redistributions of source code must retain the above copyright
19     !! notice, this list of conditions and the following disclaimer.
20     !!
21     !! 3. Redistributions in binary form must reproduce the above copyright
22     !! notice, this list of conditions and the following disclaimer in the
23     !! documentation and/or other materials provided with the
24     !! distribution.
25     !!
26     !! This software is provided "AS IS," without a warranty of any
27     !! kind. All express or implied conditions, representations and
28     !! warranties, including any implied warranty of merchantability,
29     !! fitness for a particular purpose or non-infringement, are hereby
30     !! excluded. The University of Notre Dame and its licensors shall not
31     !! be liable for any damages suffered by licensee as a result of
32     !! using, modifying or distributing the software or its
33     !! derivatives. In no event will the University of Notre Dame or its
34     !! licensors be liable for any lost revenue, profit or data, or for
35     !! direct, indirect, special, consequential, incidental or punitive
36     !! damages, however caused and regardless of the theory of liability,
37     !! arising out of the use of or inability to use software, even if the
38     !! University of Notre Dame has been advised of the possibility of
39     !! such damages.
40     !!
41    
42     module electrostatic_module
43 gezelter 507
44 gezelter 411 use force_globals
45     use definitions
46     use atype_module
47     use vector_class
48     use simulation
49     use status
50     #ifdef IS_MPI
51     use mpiSimulation
52     #endif
53     implicit none
54    
55     PRIVATE
56    
57 chuckv 656
58 gezelter 602 #define __FORTRAN90
59 chuckv 656 #include "UseTheForce/DarkSide/fInteractionMap.h"
60 gezelter 602 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61    
62 chuckv 656
63 gezelter 434 !! these prefactors convert the multipole interactions into kcal / mol
64     !! all were computed assuming distances are measured in angstroms
65     !! Charge-Charge, assuming charges are measured in electrons
66 gezelter 411 real(kind=dp), parameter :: pre11 = 332.0637778_dp
67 gezelter 434 !! Charge-Dipole, assuming charges are measured in electrons, and
68     !! dipoles are measured in debyes
69     real(kind=dp), parameter :: pre12 = 69.13373_dp
70     !! Dipole-Dipole, assuming dipoles are measured in debyes
71     real(kind=dp), parameter :: pre22 = 14.39325_dp
72     !! Charge-Quadrupole, assuming charges are measured in electrons, and
73     !! quadrupoles are measured in 10^-26 esu cm^2
74     !! This unit is also known affectionately as an esu centi-barn.
75     real(kind=dp), parameter :: pre14 = 69.13373_dp
76 gezelter 411
77 gezelter 602 !! variables to handle different summation methods for long-range electrostatics:
78     integer, save :: summationMethod = NONE
79 chrisfen 603 logical, save :: summationMethodChecked = .false.
80 gezelter 602 real(kind=DP), save :: defaultCutoff = 0.0_DP
81 chrisfen 682 real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82 gezelter 602 logical, save :: haveDefaultCutoff = .false.
83     real(kind=DP), save :: dampingAlpha = 0.0_DP
84     logical, save :: haveDampingAlpha = .false.
85 chrisfen 682 real(kind=DP), save :: dielectric = 1.0_DP
86 gezelter 602 logical, save :: haveDielectric = .false.
87     real(kind=DP), save :: constERFC = 0.0_DP
88     real(kind=DP), save :: constEXP = 0.0_DP
89     logical, save :: haveDWAconstants = .false.
90 chrisfen 682 real(kind=dp), save :: rcuti = 0.0_DP
91     real(kind=dp), save :: rcuti2 = 0.0_DP
92     real(kind=dp), save :: rcuti3 = 0.0_DP
93     real(kind=dp), save :: rcuti4 = 0.0_DP
94     real(kind=dp), save :: alphaPi = 0.0_DP
95     real(kind=dp), save :: invRootPi = 0.0_DP
96     real(kind=dp), save :: rrf = 1.0_DP
97     real(kind=dp), save :: rt = 1.0_DP
98     real(kind=dp), save :: rrfsq = 1.0_DP
99     real(kind=dp), save :: preRF = 0.0_DP
100     logical, save :: preRFCalculated = .false.
101    
102 chuckv 632 #ifdef __IFC
103     ! error function for ifc version > 7.
104 chuckv 631 double precision, external :: derfc
105 chuckv 632 #endif
106    
107 gezelter 602 public :: setElectrostaticSummationMethod
108     public :: setElectrostaticCutoffRadius
109     public :: setDampedWolfAlpha
110     public :: setReactionFieldDielectric
111 chrisfen 682 public :: setReactionFieldPrefactor
112 gezelter 411 public :: newElectrostaticType
113     public :: setCharge
114     public :: setDipoleMoment
115     public :: setSplitDipoleDistance
116     public :: setQuadrupoleMoments
117     public :: doElectrostaticPair
118     public :: getCharge
119     public :: getDipoleMoment
120 chrisfen 445 public :: pre22
121 chuckv 492 public :: destroyElectrostaticTypes
122 chrisfen 682 public :: accumulate_rf
123     public :: accumulate_self_rf
124     public :: reaction_field_final
125     public :: rf_correct_forces
126 gezelter 411
127     type :: Electrostatic
128     integer :: c_ident
129     logical :: is_Charge = .false.
130     logical :: is_Dipole = .false.
131     logical :: is_SplitDipole = .false.
132     logical :: is_Quadrupole = .false.
133 chrisfen 532 logical :: is_Tap = .false.
134 gezelter 411 real(kind=DP) :: charge = 0.0_DP
135     real(kind=DP) :: dipole_moment = 0.0_DP
136     real(kind=DP) :: split_dipole_distance = 0.0_DP
137     real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP
138     end type Electrostatic
139    
140     type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
141    
142     contains
143    
144 gezelter 602 subroutine setElectrostaticSummationMethod(the_ESM)
145     integer, intent(in) :: the_ESM
146    
147     if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
148     call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
149     endif
150    
151 chrisfen 610 summationMethod = the_ESM
152 chrisfen 626
153 gezelter 602 end subroutine setElectrostaticSummationMethod
154    
155 chrisfen 682 subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
156 gezelter 602 real(kind=dp), intent(in) :: thisRcut
157 chrisfen 682 real(kind=dp), intent(in) :: thisRsw
158 gezelter 602 defaultCutoff = thisRcut
159 chrisfen 682 rrf = defaultCutoff
160     rt = thisRsw
161 gezelter 602 haveDefaultCutoff = .true.
162     end subroutine setElectrostaticCutoffRadius
163    
164     subroutine setDampedWolfAlpha(thisAlpha)
165     real(kind=dp), intent(in) :: thisAlpha
166     dampingAlpha = thisAlpha
167     haveDampingAlpha = .true.
168     end subroutine setDampedWolfAlpha
169    
170     subroutine setReactionFieldDielectric(thisDielectric)
171     real(kind=dp), intent(in) :: thisDielectric
172     dielectric = thisDielectric
173     haveDielectric = .true.
174     end subroutine setReactionFieldDielectric
175    
176 chrisfen 682 subroutine setReactionFieldPrefactor
177     if (haveDefaultCutoff .and. haveDielectric) then
178     defaultCutoff2 = defaultCutoff*defaultCutoff
179     preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
180     ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
181     preRFCalculated = .true.
182     else if (.not.haveDefaultCutoff) then
183     call handleError("setReactionFieldPrefactor", "Default cutoff not set")
184     else
185     call handleError("setReactionFieldPrefactor", "Dielectric not set")
186     endif
187     end subroutine setReactionFieldPrefactor
188    
189 gezelter 411 subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
190 chrisfen 532 is_SplitDipole, is_Quadrupole, is_Tap, status)
191 gezelter 507
192 gezelter 411 integer, intent(in) :: c_ident
193     logical, intent(in) :: is_Charge
194     logical, intent(in) :: is_Dipole
195     logical, intent(in) :: is_SplitDipole
196     logical, intent(in) :: is_Quadrupole
197 chrisfen 532 logical, intent(in) :: is_Tap
198 gezelter 411 integer, intent(out) :: status
199     integer :: nAtypes, myATID, i, j
200    
201     status = 0
202     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
203 gezelter 507
204 gezelter 411 !! Be simple-minded and assume that we need an ElectrostaticMap that
205     !! is the same size as the total number of atom types
206    
207     if (.not.allocated(ElectrostaticMap)) then
208 gezelter 507
209 gezelter 411 nAtypes = getSize(atypes)
210 gezelter 507
211 gezelter 411 if (nAtypes == 0) then
212     status = -1
213     return
214     end if
215 gezelter 507
216 gezelter 411 if (.not. allocated(ElectrostaticMap)) then
217     allocate(ElectrostaticMap(nAtypes))
218     endif
219 gezelter 507
220 gezelter 411 end if
221    
222     if (myATID .gt. size(ElectrostaticMap)) then
223     status = -1
224     return
225     endif
226 gezelter 507
227 gezelter 411 ! set the values for ElectrostaticMap for this atom type:
228    
229     ElectrostaticMap(myATID)%c_ident = c_ident
230     ElectrostaticMap(myATID)%is_Charge = is_Charge
231     ElectrostaticMap(myATID)%is_Dipole = is_Dipole
232     ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
233     ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
234 chrisfen 532 ElectrostaticMap(myATID)%is_Tap = is_Tap
235 gezelter 507
236 gezelter 411 end subroutine newElectrostaticType
237    
238     subroutine setCharge(c_ident, charge, status)
239     integer, intent(in) :: c_ident
240     real(kind=dp), intent(in) :: charge
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 setCharge!")
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 setCharge!")
255     status = -1
256     return
257     endif
258    
259     if (.not.ElectrostaticMap(myATID)%is_Charge) then
260     call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
261     status = -1
262     return
263 gezelter 507 endif
264 gezelter 411
265     ElectrostaticMap(myATID)%charge = charge
266     end subroutine setCharge
267    
268     subroutine setDipoleMoment(c_ident, dipole_moment, status)
269     integer, intent(in) :: c_ident
270     real(kind=dp), intent(in) :: dipole_moment
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 setDipoleMoment!")
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 setDipoleMoment!")
285     status = -1
286     return
287     endif
288    
289     if (.not.ElectrostaticMap(myATID)%is_Dipole) then
290     call handleError("electrostatic", "Attempt to setDipoleMoment of an atom type that is not a dipole!")
291     status = -1
292     return
293     endif
294    
295     ElectrostaticMap(myATID)%dipole_moment = dipole_moment
296     end subroutine setDipoleMoment
297    
298     subroutine setSplitDipoleDistance(c_ident, split_dipole_distance, status)
299     integer, intent(in) :: c_ident
300     real(kind=dp), intent(in) :: split_dipole_distance
301     integer, intent(out) :: status
302     integer :: myATID
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 setSplitDipoleDistance!")
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 setSplitDipoleDistance!")
315     status = -1
316     return
317     endif
318    
319     if (.not.ElectrostaticMap(myATID)%is_SplitDipole) then
320     call handleError("electrostatic", "Attempt to setSplitDipoleDistance of an atom type that is not a splitDipole!")
321     status = -1
322     return
323     endif
324    
325     ElectrostaticMap(myATID)%split_dipole_distance = split_dipole_distance
326     end subroutine setSplitDipoleDistance
327    
328     subroutine setQuadrupoleMoments(c_ident, quadrupole_moments, status)
329     integer, intent(in) :: c_ident
330     real(kind=dp), intent(in), dimension(3) :: quadrupole_moments
331     integer, intent(out) :: status
332     integer :: myATID, i, j
333    
334     status = 0
335     myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
336    
337     if (.not.allocated(ElectrostaticMap)) then
338     call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
339     status = -1
340     return
341     end if
342    
343     if (myATID .gt. size(ElectrostaticMap)) then
344     call handleError("electrostatic", "ElectrostaticMap was found to be too small during setQuadrupoleMoments!")
345     status = -1
346     return
347     endif
348    
349     if (.not.ElectrostaticMap(myATID)%is_Quadrupole) then
350     call handleError("electrostatic", "Attempt to setQuadrupoleMoments of an atom type that is not a quadrupole!")
351     status = -1
352     return
353     endif
354 gezelter 507
355 gezelter 411 do i = 1, 3
356 gezelter 507 ElectrostaticMap(myATID)%quadrupole_moments(i) = &
357     quadrupole_moments(i)
358     enddo
359 gezelter 411
360     end subroutine setQuadrupoleMoments
361    
362 gezelter 507
363 gezelter 411 function getCharge(atid) result (c)
364     integer, intent(in) :: atid
365     integer :: localError
366     real(kind=dp) :: c
367 gezelter 507
368 gezelter 411 if (.not.allocated(ElectrostaticMap)) then
369     call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
370     return
371     end if
372 gezelter 507
373 gezelter 411 if (.not.ElectrostaticMap(atid)%is_Charge) then
374     call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
375     return
376     endif
377 gezelter 507
378 gezelter 411 c = ElectrostaticMap(atid)%charge
379     end function getCharge
380    
381     function getDipoleMoment(atid) result (dm)
382     integer, intent(in) :: atid
383     integer :: localError
384     real(kind=dp) :: dm
385 gezelter 507
386 gezelter 411 if (.not.allocated(ElectrostaticMap)) then
387     call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
388     return
389     end if
390 gezelter 507
391 gezelter 411 if (.not.ElectrostaticMap(atid)%is_Dipole) then
392     call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
393     return
394     endif
395 gezelter 507
396 gezelter 411 dm = ElectrostaticMap(atid)%dipole_moment
397     end function getDipoleMoment
398    
399 gezelter 602 subroutine checkSummationMethod()
400    
401 chrisfen 607 if (.not.haveDefaultCutoff) then
402     call handleError("checkSummationMethod", "no Default Cutoff set!")
403     endif
404    
405     rcuti = 1.0d0 / defaultCutoff
406     rcuti2 = rcuti*rcuti
407     rcuti3 = rcuti2*rcuti
408     rcuti4 = rcuti2*rcuti2
409    
410 gezelter 602 if (summationMethod .eq. DAMPED_WOLF) then
411     if (.not.haveDWAconstants) then
412    
413     if (.not.haveDampingAlpha) then
414     call handleError("checkSummationMethod", "no Damping Alpha set!")
415     endif
416    
417 chrisfen 603 if (.not.haveDefaultCutoff) then
418     call handleError("checkSummationMethod", "no Default Cutoff set!")
419     endif
420    
421     constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
422 chuckv 631 constERFC = derfc(dampingAlpha*defaultCutoff)
423 chrisfen 640 invRootPi = 0.56418958354775628695d0
424     alphaPi = 2*dampingAlpha*invRootPi
425 chrisfen 644
426 gezelter 602 haveDWAconstants = .true.
427     endif
428     endif
429    
430 chrisfen 603 if (summationMethod .eq. REACTION_FIELD) then
431     if (.not.haveDielectric) then
432     call handleError("checkSummationMethod", "no reaction field Dielectric set!")
433     endif
434     endif
435    
436     summationMethodChecked = .true.
437 gezelter 602 end subroutine checkSummationMethod
438    
439    
440    
441 gezelter 411 subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
442 chrisfen 607 vpair, fpair, pot, eFrame, f, t, do_pot)
443 gezelter 507
444 gezelter 411 logical, intent(in) :: do_pot
445 gezelter 507
446 gezelter 411 integer, intent(in) :: atom1, atom2
447     integer :: localError
448    
449 chrisfen 607 real(kind=dp), intent(in) :: rij, r2, sw
450 gezelter 411 real(kind=dp), intent(in), dimension(3) :: d
451     real(kind=dp), intent(inout) :: vpair
452     real(kind=dp), intent(inout), dimension(3) :: fpair
453    
454 chrisfen 626 real( kind = dp ) :: pot
455 gezelter 411 real( kind = dp ), dimension(9,nLocal) :: eFrame
456     real( kind = dp ), dimension(3,nLocal) :: f
457     real( kind = dp ), dimension(3,nLocal) :: t
458 gezelter 507
459 gezelter 439 real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
460     real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
461     real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
462     real (kind = dp), dimension(3) :: dudux_j, duduy_j, duduz_j
463 gezelter 411
464     logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
465     logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
466 chrisfen 532 logical :: i_is_Tap, j_is_Tap
467 gezelter 411 integer :: me1, me2, id1, id2
468     real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
469 gezelter 439 real (kind=dp) :: qxx_i, qyy_i, qzz_i
470     real (kind=dp) :: qxx_j, qyy_j, qzz_j
471     real (kind=dp) :: cx_i, cy_i, cz_i
472     real (kind=dp) :: cx_j, cy_j, cz_j
473     real (kind=dp) :: cx2, cy2, cz2
474 gezelter 411 real (kind=dp) :: ct_i, ct_j, ct_ij, a1
475 gezelter 421 real (kind=dp) :: riji, ri, ri2, ri3, ri4
476 chrisfen 597 real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
477 gezelter 421 real (kind=dp) :: xhat, yhat, zhat
478 gezelter 411 real (kind=dp) :: dudx, dudy, dudz
479 chrisfen 626 real (kind=dp) :: scale, sc2, bigR
480 chrisfen 640 real (kind=dp) :: varERFC, varEXP
481 chrisfen 644 real (kind=dp) :: limScale
482 gezelter 411
483     if (.not.allocated(ElectrostaticMap)) then
484     call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
485     return
486     end if
487    
488 gezelter 602 if (.not.summationMethodChecked) then
489     call checkSummationMethod()
490 chrisfen 611
491 gezelter 602 endif
492    
493    
494 gezelter 411 #ifdef IS_MPI
495     me1 = atid_Row(atom1)
496     me2 = atid_Col(atom2)
497     #else
498     me1 = atid(atom1)
499     me2 = atid(atom2)
500     #endif
501    
502     !! some variables we'll need independent of electrostatic type:
503    
504     riji = 1.0d0 / rij
505 chrisfen 644
506 gezelter 421 xhat = d(1) * riji
507     yhat = d(2) * riji
508     zhat = d(3) * riji
509 gezelter 411
510     !! logicals
511     i_is_Charge = ElectrostaticMap(me1)%is_Charge
512     i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
513     i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
514     i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
515 chrisfen 532 i_is_Tap = ElectrostaticMap(me1)%is_Tap
516 gezelter 411
517     j_is_Charge = ElectrostaticMap(me2)%is_Charge
518     j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
519     j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
520     j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
521 chrisfen 532 j_is_Tap = ElectrostaticMap(me2)%is_Tap
522 gezelter 411
523     if (i_is_Charge) then
524     q_i = ElectrostaticMap(me1)%charge
525     endif
526 gezelter 507
527 gezelter 411 if (i_is_Dipole) then
528     mu_i = ElectrostaticMap(me1)%dipole_moment
529     #ifdef IS_MPI
530 gezelter 439 uz_i(1) = eFrame_Row(3,atom1)
531     uz_i(2) = eFrame_Row(6,atom1)
532     uz_i(3) = eFrame_Row(9,atom1)
533 gezelter 411 #else
534 gezelter 439 uz_i(1) = eFrame(3,atom1)
535     uz_i(2) = eFrame(6,atom1)
536     uz_i(3) = eFrame(9,atom1)
537 gezelter 411 #endif
538 gezelter 439 ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
539 gezelter 411
540     if (i_is_SplitDipole) then
541     d_i = ElectrostaticMap(me1)%split_dipole_distance
542     endif
543 gezelter 507
544 gezelter 411 endif
545    
546 gezelter 439 if (i_is_Quadrupole) then
547     qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1)
548     qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2)
549     qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3)
550     #ifdef IS_MPI
551     ux_i(1) = eFrame_Row(1,atom1)
552     ux_i(2) = eFrame_Row(4,atom1)
553     ux_i(3) = eFrame_Row(7,atom1)
554     uy_i(1) = eFrame_Row(2,atom1)
555     uy_i(2) = eFrame_Row(5,atom1)
556     uy_i(3) = eFrame_Row(8,atom1)
557     uz_i(1) = eFrame_Row(3,atom1)
558     uz_i(2) = eFrame_Row(6,atom1)
559     uz_i(3) = eFrame_Row(9,atom1)
560     #else
561     ux_i(1) = eFrame(1,atom1)
562     ux_i(2) = eFrame(4,atom1)
563     ux_i(3) = eFrame(7,atom1)
564     uy_i(1) = eFrame(2,atom1)
565     uy_i(2) = eFrame(5,atom1)
566     uy_i(3) = eFrame(8,atom1)
567     uz_i(1) = eFrame(3,atom1)
568     uz_i(2) = eFrame(6,atom1)
569     uz_i(3) = eFrame(9,atom1)
570     #endif
571     cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
572     cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
573     cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
574     endif
575    
576 gezelter 411 if (j_is_Charge) then
577     q_j = ElectrostaticMap(me2)%charge
578     endif
579 gezelter 507
580 gezelter 411 if (j_is_Dipole) then
581     mu_j = ElectrostaticMap(me2)%dipole_moment
582     #ifdef IS_MPI
583 gezelter 439 uz_j(1) = eFrame_Col(3,atom2)
584     uz_j(2) = eFrame_Col(6,atom2)
585     uz_j(3) = eFrame_Col(9,atom2)
586 gezelter 411 #else
587 gezelter 439 uz_j(1) = eFrame(3,atom2)
588     uz_j(2) = eFrame(6,atom2)
589     uz_j(3) = eFrame(9,atom2)
590 gezelter 411 #endif
591 chrisfen 465 ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
592 gezelter 411
593     if (j_is_SplitDipole) then
594     d_j = ElectrostaticMap(me2)%split_dipole_distance
595     endif
596     endif
597    
598 gezelter 439 if (j_is_Quadrupole) then
599     qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1)
600     qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2)
601     qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3)
602     #ifdef IS_MPI
603     ux_j(1) = eFrame_Col(1,atom2)
604     ux_j(2) = eFrame_Col(4,atom2)
605     ux_j(3) = eFrame_Col(7,atom2)
606     uy_j(1) = eFrame_Col(2,atom2)
607     uy_j(2) = eFrame_Col(5,atom2)
608     uy_j(3) = eFrame_Col(8,atom2)
609     uz_j(1) = eFrame_Col(3,atom2)
610     uz_j(2) = eFrame_Col(6,atom2)
611     uz_j(3) = eFrame_Col(9,atom2)
612     #else
613     ux_j(1) = eFrame(1,atom2)
614     ux_j(2) = eFrame(4,atom2)
615     ux_j(3) = eFrame(7,atom2)
616     uy_j(1) = eFrame(2,atom2)
617     uy_j(2) = eFrame(5,atom2)
618     uy_j(3) = eFrame(8,atom2)
619     uz_j(1) = eFrame(3,atom2)
620     uz_j(2) = eFrame(6,atom2)
621     uz_j(3) = eFrame(9,atom2)
622     #endif
623     cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
624     cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
625     cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
626     endif
627 chrisfen 554
628 gezelter 411 epot = 0.0_dp
629     dudx = 0.0_dp
630     dudy = 0.0_dp
631     dudz = 0.0_dp
632    
633 gezelter 439 dudux_i = 0.0_dp
634     duduy_i = 0.0_dp
635     duduz_i = 0.0_dp
636 gezelter 411
637 gezelter 439 dudux_j = 0.0_dp
638     duduy_j = 0.0_dp
639     duduz_j = 0.0_dp
640 gezelter 411
641     if (i_is_Charge) then
642    
643     if (j_is_Charge) then
644 gezelter 507
645 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
646 chrisfen 626
647 chrisfen 597 vterm = pre11 * q_i * q_j * (riji - rcuti)
648     vpair = vpair + vterm
649 chrisfen 626 epot = epot + sw*vterm
650 chrisfen 597
651 chrisfen 644 dudr = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
652 chrisfen 597
653     dudx = dudx + dudr * d(1)
654     dudy = dudy + dudr * d(2)
655     dudz = dudz + dudr * d(3)
656 gezelter 411
657 chrisfen 640 elseif (summationMethod .eq. DAMPED_WOLF) then
658    
659     varERFC = derfc(dampingAlpha*rij)
660     varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
661     vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
662     vpair = vpair + vterm
663     epot = epot + sw*vterm
664    
665 chrisfen 644 dudr = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
666     + alphaPi*varEXP) &
667     - (constERFC*rcuti2 &
668     + alphaPi*constEXP)) )
669 chrisfen 640
670     dudx = dudx + dudr * d(1)
671     dudy = dudy + dudr * d(2)
672     dudz = dudz + dudr * d(3)
673    
674 chrisfen 597 else
675 chrisfen 626
676 chrisfen 597 vterm = pre11 * q_i * q_j * riji
677     vpair = vpair + vterm
678 chrisfen 626 epot = epot + sw*vterm
679 chrisfen 597
680     dudr = - sw * vterm * riji
681    
682     dudx = dudx + dudr * xhat
683     dudy = dudy + dudr * yhat
684     dudz = dudz + dudr * zhat
685    
686     endif
687    
688 gezelter 411 endif
689    
690     if (j_is_Dipole) then
691    
692 chrisfen 626 pref = pre12 * q_i * mu_j
693 gezelter 411
694 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
695 chrisfen 597 ri2 = riji * riji
696     ri3 = ri2 * riji
697 gezelter 507
698 chrisfen 626 pref = pre12 * q_i * mu_j
699 chrisfen 597 vterm = - pref * ct_j * (ri2 - rcuti2)
700 chrisfen 626 vpair = vpair + vterm
701     epot = epot + sw*vterm
702 chrisfen 597
703     !! this has a + sign in the () because the rij vector is
704     !! r_j - r_i and the charge-dipole potential takes the origin
705     !! as the point dipole, which is atom j in this case.
706    
707 chrisfen 626 dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
708 chrisfen 597 - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
709 chrisfen 626 dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
710 chrisfen 597 - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
711 chrisfen 626 dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
712 chrisfen 597 - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
713    
714 chrisfen 626 duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
715     duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
716     duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
717 gezelter 411
718 chrisfen 597 else
719     if (j_is_SplitDipole) then
720     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
721     ri = 1.0_dp / BigR
722     scale = rij * ri
723     else
724     ri = riji
725     scale = 1.0_dp
726     endif
727    
728     ri2 = ri * ri
729     ri3 = ri2 * ri
730     sc2 = scale * scale
731 chrisfen 626
732     pref = pre12 * q_i * mu_j
733 chrisfen 597 vterm = - pref * ct_j * ri2 * scale
734 chrisfen 626 vpair = vpair + vterm
735     epot = epot + sw*vterm
736 chrisfen 597
737     !! this has a + sign in the () because the rij vector is
738     !! r_j - r_i and the charge-dipole potential takes the origin
739     !! as the point dipole, which is atom j in this case.
740    
741 chrisfen 626 dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
742     dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
743     dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
744 chrisfen 597
745 chrisfen 626 duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
746     duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
747     duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
748 gezelter 411
749 chrisfen 597 endif
750 gezelter 411 endif
751 gezelter 421
752 gezelter 439 if (j_is_Quadrupole) then
753     ri2 = riji * riji
754     ri3 = ri2 * riji
755 gezelter 440 ri4 = ri2 * ri2
756 gezelter 439 cx2 = cx_j * cx_j
757     cy2 = cy_j * cy_j
758     cz2 = cz_j * cz_j
759    
760 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
761 chrisfen 626 pref = pre14 * q_i / 3.0_dp
762 chrisfen 597 vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
763     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
764     qzz_j * (3.0_dp*cz2 - 1.0_dp) )
765     vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
766     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
767     qzz_j * (3.0_dp*cz2 - 1.0_dp) )
768 chrisfen 626 vpair = vpair + ( vterm1 - vterm2 )
769     epot = epot + sw*( vterm1 - vterm2 )
770 chrisfen 597
771     dudx = dudx - (5.0_dp * &
772 chrisfen 626 (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
773 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
774     qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
775     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
776     qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
777     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
778     qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
779     dudy = dudy - (5.0_dp * &
780 chrisfen 626 (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
781 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
782     qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
783     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
784     qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
785     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
786     qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
787     dudz = dudz - (5.0_dp * &
788 chrisfen 626 (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
789 chrisfen 597 (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
790     qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
791     (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
792     qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
793     (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
794     qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
795    
796 chrisfen 626 dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
797 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
798 chrisfen 626 dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
799 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
800 chrisfen 626 dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
801 chrisfen 597 rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
802    
803 chrisfen 626 duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
804 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
805 chrisfen 626 duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
806 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
807 chrisfen 626 duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
808 chrisfen 597 rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
809    
810 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
811 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
812 chrisfen 626 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
813 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
814 chrisfen 626 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
815 chrisfen 597 rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
816    
817     else
818 chrisfen 626 pref = pre14 * q_i / 3.0_dp
819 chrisfen 597 vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
820     qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
821     qzz_j * (3.0_dp*cz2 - 1.0_dp))
822 chrisfen 626 vpair = vpair + vterm
823     epot = epot + sw*vterm
824 chrisfen 597
825 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
826 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
827     qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
828     qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
829 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
830 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
831     qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
832     qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
833 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
834 chrisfen 597 qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
835     qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
836     qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
837    
838 chrisfen 626 dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
839     dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
840     dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
841 chrisfen 597
842 chrisfen 626 duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
843     duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
844     duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
845 chrisfen 597
846 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
847     duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
848     duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
849 chrisfen 597
850     endif
851 gezelter 439 endif
852 gezelter 411 endif
853 gezelter 507
854 gezelter 411 if (i_is_Dipole) then
855 gezelter 507
856 gezelter 411 if (j_is_Charge) then
857 chrisfen 626
858     pref = pre12 * q_j * mu_i
859    
860 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
861 chrisfen 597 ri2 = riji * riji
862     ri3 = ri2 * riji
863 gezelter 507
864 chrisfen 626 pref = pre12 * q_j * mu_i
865 chrisfen 597 vterm = pref * ct_i * (ri2 - rcuti2)
866 chrisfen 626 vpair = vpair + vterm
867     epot = epot + sw*vterm
868 chrisfen 597
869     !! this has a + sign in the () because the rij vector is
870     !! r_j - r_i and the charge-dipole potential takes the origin
871     !! as the point dipole, which is atom j in this case.
872    
873 chrisfen 626 dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
874 chrisfen 597 - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
875 chrisfen 626 dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
876 chrisfen 597 - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
877 chrisfen 626 dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
878 chrisfen 597 - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
879    
880 chrisfen 626 duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
881     duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
882     duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
883 gezelter 411
884 chrisfen 597 else
885     if (i_is_SplitDipole) then
886 gezelter 421 BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
887     ri = 1.0_dp / BigR
888 chrisfen 597 scale = rij * ri
889     else
890 gezelter 421 ri = riji
891     scale = 1.0_dp
892     endif
893 chrisfen 597
894     ri2 = ri * ri
895     ri3 = ri2 * ri
896     sc2 = scale * scale
897 chrisfen 626
898     pref = pre12 * q_j * mu_i
899 chrisfen 597 vterm = pref * ct_i * ri2 * scale
900 chrisfen 626 vpair = vpair + vterm
901     epot = epot + sw*vterm
902 chrisfen 597
903 chrisfen 626 dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
904     dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
905     dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
906 chrisfen 597
907 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
908     duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
909     duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
910 gezelter 421 endif
911 chrisfen 597 endif
912 chrisfen 626
913 chrisfen 597 if (j_is_Dipole) then
914 gezelter 421
915 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
916 chrisfen 597 ri2 = riji * riji
917     ri3 = ri2 * riji
918     ri4 = ri2 * ri2
919 gezelter 507
920 chrisfen 626 pref = pre22 * mu_i * mu_j
921 chrisfen 597 vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
922 chrisfen 626 vpair = vpair + vterm
923     epot = epot + sw*vterm
924 chrisfen 597
925     a1 = 5.0d0 * ct_i * ct_j - ct_ij
926    
927 chrisfen 626 dudx = dudx + sw*pref*3.0d0*ri4 &
928     * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
929     - sw*pref*3.0d0*rcuti4 &
930     * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
931     dudy = dudy + sw*pref*3.0d0*ri4 &
932     * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
933     - sw*pref*3.0d0*rcuti4 &
934     * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
935     dudz = dudz + sw*pref*3.0d0*ri4 &
936     * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
937     - sw*pref*3.0d0*rcuti4 &
938     * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
939 chrisfen 597
940 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
941 chrisfen 597 - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
942 chrisfen 626 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
943 chrisfen 597 - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
944 chrisfen 626 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
945 chrisfen 597 - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
946 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
947 chrisfen 597 - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
948 chrisfen 626 duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
949 chrisfen 597 - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
950 chrisfen 626 duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
951 chrisfen 597 - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
952 chrisfen 626
953 chrisfen 597 else
954     if (i_is_SplitDipole) then
955     if (j_is_SplitDipole) then
956     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
957     else
958     BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
959     endif
960     ri = 1.0_dp / BigR
961     scale = rij * ri
962     else
963     if (j_is_SplitDipole) then
964     BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
965     ri = 1.0_dp / BigR
966     scale = rij * ri
967     else
968     ri = riji
969     scale = 1.0_dp
970     endif
971     endif
972    
973     ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
974    
975     ri2 = ri * ri
976     ri3 = ri2 * ri
977     ri4 = ri2 * ri2
978     sc2 = scale * scale
979    
980 chrisfen 626 pref = pre22 * mu_i * mu_j
981 chrisfen 597 vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
982 chrisfen 626 vpair = vpair + vterm
983     epot = epot + sw*vterm
984 chrisfen 597
985     a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
986    
987 chrisfen 626 dudx = dudx + sw*pref*3.0d0*ri4*scale &
988     *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989     dudy = dudy + sw*pref*3.0d0*ri4*scale &
990     *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991     dudz = dudz + sw*pref*3.0d0*ri4*scale &
992     *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
993 chrisfen 597
994 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
995     *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
996     duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
997     *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
998     duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
999     *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1000 chrisfen 597
1001 chrisfen 626 duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1002     *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1003     duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1004     *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1005     duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1006     *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1007 chrisfen 597 endif
1008 gezelter 411 endif
1009     endif
1010 gezelter 439
1011     if (i_is_Quadrupole) then
1012     if (j_is_Charge) then
1013 gezelter 507
1014 gezelter 439 ri2 = riji * riji
1015     ri3 = ri2 * riji
1016 gezelter 440 ri4 = ri2 * ri2
1017 gezelter 439 cx2 = cx_i * cx_i
1018     cy2 = cy_i * cy_i
1019     cz2 = cz_i * cz_i
1020 gezelter 507
1021 chrisfen 611 if (summationMethod .eq. UNDAMPED_WOLF) then
1022 chrisfen 626 pref = pre14 * q_j / 3.0_dp
1023 chrisfen 597 vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1024     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1025     qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1026     vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1027     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1028     qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1029 chrisfen 626 vpair = vpair + ( vterm1 - vterm2 )
1030     epot = epot + sw*( vterm1 - vterm2 )
1031 chrisfen 597
1032 chrisfen 626 dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1033     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1034 chrisfen 597 qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1035     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1036     qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1037     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1038     qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1039 chrisfen 626 dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1040     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1041 chrisfen 597 qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1042     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1043     qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1044     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1045     qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1046 chrisfen 626 dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1047     sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1048 chrisfen 597 qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1049     (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1050     qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1051     (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1052     qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1053    
1054 chrisfen 626 dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1055 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1056 chrisfen 626 dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1057 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1058 chrisfen 626 dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1059 chrisfen 597 rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1060    
1061 chrisfen 626 duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1062 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1063 chrisfen 626 duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1064 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1065 chrisfen 626 duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1066 chrisfen 597 rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1067    
1068 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1069 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1070 chrisfen 626 duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1071 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1072 chrisfen 626 duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1073 chrisfen 597 rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1074 gezelter 507
1075 chrisfen 597 else
1076 chrisfen 626 pref = pre14 * q_j / 3.0_dp
1077 chrisfen 597 vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1078     qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1079     qzz_i * (3.0_dp*cz2 - 1.0_dp))
1080 chrisfen 626 vpair = vpair + vterm
1081     epot = epot + sw*vterm
1082 chrisfen 597
1083 chrisfen 626 dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1084 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1085     qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1086     qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1087 chrisfen 626 dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1088 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1089     qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1090     qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1091 chrisfen 626 dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1092 chrisfen 597 qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1093     qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1094     qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1095    
1096 chrisfen 626 dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1097     dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1098     dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1099 chrisfen 597
1100 chrisfen 626 duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1101     duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1102     duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1103 chrisfen 597
1104 chrisfen 626 duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1105     duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1106     duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1107 chrisfen 597 endif
1108 gezelter 439 endif
1109     endif
1110 gezelter 507
1111    
1112 gezelter 411 if (do_pot) then
1113     #ifdef IS_MPI
1114 chuckv 656 pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1115     pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1116 gezelter 411 #else
1117     pot = pot + epot
1118     #endif
1119     endif
1120 gezelter 507
1121 gezelter 411 #ifdef IS_MPI
1122     f_Row(1,atom1) = f_Row(1,atom1) + dudx
1123     f_Row(2,atom1) = f_Row(2,atom1) + dudy
1124     f_Row(3,atom1) = f_Row(3,atom1) + dudz
1125 gezelter 507
1126 gezelter 411 f_Col(1,atom2) = f_Col(1,atom2) - dudx
1127     f_Col(2,atom2) = f_Col(2,atom2) - dudy
1128     f_Col(3,atom2) = f_Col(3,atom2) - dudz
1129 gezelter 507
1130 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1131 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1132     t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1133     t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1134 gezelter 411 endif
1135 gezelter 439 if (i_is_Quadrupole) then
1136     t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1137     t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1138     t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1139 gezelter 411
1140 gezelter 439 t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1141     t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1142     t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1143     endif
1144    
1145 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1146 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1147     t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1148     t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1149 gezelter 411 endif
1150 gezelter 439 if (j_is_Quadrupole) then
1151     t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1152     t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1153     t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1154 gezelter 411
1155 gezelter 439 t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1156     t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1157     t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1158     endif
1159    
1160 gezelter 411 #else
1161     f(1,atom1) = f(1,atom1) + dudx
1162     f(2,atom1) = f(2,atom1) + dudy
1163     f(3,atom1) = f(3,atom1) + dudz
1164 gezelter 507
1165 gezelter 411 f(1,atom2) = f(1,atom2) - dudx
1166     f(2,atom2) = f(2,atom2) - dudy
1167     f(3,atom2) = f(3,atom2) - dudz
1168 gezelter 507
1169 gezelter 411 if (i_is_Dipole .or. i_is_Quadrupole) then
1170 gezelter 439 t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1171     t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1172     t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1173 gezelter 411 endif
1174 gezelter 439 if (i_is_Quadrupole) then
1175     t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
1176     t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
1177     t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)
1178    
1179     t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
1180     t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
1181     t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
1182     endif
1183    
1184 gezelter 411 if (j_is_Dipole .or. j_is_Quadrupole) then
1185 gezelter 439 t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1186     t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1187     t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1188 gezelter 411 endif
1189 gezelter 439 if (j_is_Quadrupole) then
1190     t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
1191     t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
1192     t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)
1193    
1194     t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
1195     t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
1196     t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
1197     endif
1198    
1199 gezelter 411 #endif
1200 gezelter 507
1201 gezelter 411 #ifdef IS_MPI
1202     id1 = AtomRowToGlobal(atom1)
1203     id2 = AtomColToGlobal(atom2)
1204     #else
1205     id1 = atom1
1206     id2 = atom2
1207     #endif
1208    
1209     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1210 gezelter 507
1211 gezelter 411 fpair(1) = fpair(1) + dudx
1212     fpair(2) = fpair(2) + dudy
1213     fpair(3) = fpair(3) + dudz
1214    
1215     endif
1216    
1217     return
1218     end subroutine doElectrostaticPair
1219 chuckv 492
1220     subroutine destroyElectrostaticTypes()
1221    
1222 gezelter 507 if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1223    
1224 chuckv 492 end subroutine destroyElectrostaticTypes
1225    
1226 chrisfen 682 subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
1227    
1228     integer, intent(in) :: atom1, atom2
1229     real (kind = dp), intent(in) :: rij
1230     real (kind = dp), dimension(9,nLocal) :: eFrame
1231    
1232     integer :: me1, me2
1233     real (kind = dp), intent(in) :: taper
1234     real (kind = dp):: mu1, mu2
1235     real (kind = dp), dimension(3) :: ul1
1236     real (kind = dp), dimension(3) :: ul2
1237    
1238     integer :: localError
1239    
1240     #ifdef IS_MPI
1241     me1 = atid_Row(atom1)
1242     ul1(1) = eFrame_Row(3,atom1)
1243     ul1(2) = eFrame_Row(6,atom1)
1244     ul1(3) = eFrame_Row(9,atom1)
1245    
1246     me2 = atid_Col(atom2)
1247     ul2(1) = eFrame_Col(3,atom2)
1248     ul2(2) = eFrame_Col(6,atom2)
1249     ul2(3) = eFrame_Col(9,atom2)
1250     #else
1251     me1 = atid(atom1)
1252     ul1(1) = eFrame(3,atom1)
1253     ul1(2) = eFrame(6,atom1)
1254     ul1(3) = eFrame(9,atom1)
1255    
1256     me2 = atid(atom2)
1257     ul2(1) = eFrame(3,atom2)
1258     ul2(2) = eFrame(6,atom2)
1259     ul2(3) = eFrame(9,atom2)
1260     #endif
1261    
1262     mu1 = getDipoleMoment(me1)
1263     mu2 = getDipoleMoment(me2)
1264    
1265     #ifdef IS_MPI
1266     rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
1267     rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
1268     rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
1269    
1270     rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
1271     rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
1272     rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
1273     #else
1274     rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
1275     rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
1276     rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
1277    
1278     rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
1279     rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
1280     rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper
1281     #endif
1282     return
1283     end subroutine accumulate_rf
1284    
1285     subroutine accumulate_self_rf(atom1, mu1, eFrame)
1286    
1287     integer, intent(in) :: atom1
1288     real(kind=dp), intent(in) :: mu1
1289     real(kind=dp), dimension(9,nLocal) :: eFrame
1290    
1291     !! should work for both MPI and non-MPI version since this is not pairwise.
1292     rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
1293     rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
1294     rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
1295    
1296     return
1297     end subroutine accumulate_self_rf
1298    
1299     subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
1300    
1301     integer, intent(in) :: a1
1302     real (kind=dp), intent(in) :: mu1
1303     real (kind=dp), intent(inout) :: rfpot
1304     logical, intent(in) :: do_pot
1305     real (kind = dp), dimension(9,nLocal) :: eFrame
1306     real (kind = dp), dimension(3,nLocal) :: t
1307    
1308     integer :: localError
1309    
1310     if (.not.preRFCalculated) then
1311     call setReactionFieldPrefactor()
1312     endif
1313    
1314     ! compute torques on dipoles:
1315     ! pre converts from mu in units of debye to kcal/mol
1316    
1317     ! The torque contribution is dipole cross reaction_field
1318    
1319     t(1,a1) = t(1,a1) + preRF*mu1*(eFrame(6,a1)*rf(3,a1) - &
1320     eFrame(9,a1)*rf(2,a1))
1321     t(2,a1) = t(2,a1) + preRF*mu1*(eFrame(9,a1)*rf(1,a1) - &
1322     eFrame(3,a1)*rf(3,a1))
1323     t(3,a1) = t(3,a1) + preRF*mu1*(eFrame(3,a1)*rf(2,a1) - &
1324     eFrame(6,a1)*rf(1,a1))
1325    
1326     ! the potential contribution is -1/2 dipole dot reaction_field
1327    
1328     if (do_pot) then
1329     rfpot = rfpot - 0.5d0 * preRF * mu1 * &
1330     (rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + &
1331     rf(3,a1)*eFrame(9,a1))
1332     endif
1333    
1334     return
1335     end subroutine reaction_field_final
1336    
1337     subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
1338    
1339     integer, intent(in) :: atom1, atom2
1340     real(kind=dp), dimension(3), intent(in) :: d
1341     real(kind=dp), intent(in) :: rij, taper
1342     real( kind = dp ), dimension(9,nLocal) :: eFrame
1343     real( kind = dp ), dimension(3,nLocal) :: f
1344     real( kind = dp ), dimension(3), intent(inout) :: fpair
1345    
1346     real (kind = dp), dimension(3) :: ul1
1347     real (kind = dp), dimension(3) :: ul2
1348     real (kind = dp) :: dtdr
1349     real (kind = dp) :: dudx, dudy, dudz, u1dotu2
1350     integer :: me1, me2, id1, id2
1351     real (kind = dp) :: mu1, mu2
1352    
1353     integer :: localError
1354    
1355     if (.not.preRFCalculated) then
1356     call setReactionFieldPrefactor()
1357     endif
1358    
1359     if (rij.le.rrf) then
1360    
1361     if (rij.lt.rt) then
1362     dtdr = 0.0d0
1363     else
1364     ! write(*,*) 'rf correct in taper region'
1365     dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
1366     endif
1367    
1368     #ifdef IS_MPI
1369     me1 = atid_Row(atom1)
1370     ul1(1) = eFrame_Row(3,atom1)
1371     ul1(2) = eFrame_Row(6,atom1)
1372     ul1(3) = eFrame_Row(9,atom1)
1373    
1374     me2 = atid_Col(atom2)
1375     ul2(1) = eFrame_Col(3,atom2)
1376     ul2(2) = eFrame_Col(6,atom2)
1377     ul2(3) = eFrame_Col(9,atom2)
1378     #else
1379     me1 = atid(atom1)
1380     ul1(1) = eFrame(3,atom1)
1381     ul1(2) = eFrame(6,atom1)
1382     ul1(3) = eFrame(9,atom1)
1383    
1384     me2 = atid(atom2)
1385     ul2(1) = eFrame(3,atom2)
1386     ul2(2) = eFrame(6,atom2)
1387     ul2(3) = eFrame(9,atom2)
1388     #endif
1389    
1390     mu1 = getDipoleMoment(me1)
1391     mu2 = getDipoleMoment(me2)
1392    
1393     u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
1394    
1395     dudx = - preRF*mu1*mu2*u1dotu2*dtdr*d(1)/rij
1396     dudy = - preRF*mu1*mu2*u1dotu2*dtdr*d(2)/rij
1397     dudz = - preRF*mu1*mu2*u1dotu2*dtdr*d(3)/rij
1398    
1399     #ifdef IS_MPI
1400     f_Row(1,atom1) = f_Row(1,atom1) + dudx
1401     f_Row(2,atom1) = f_Row(2,atom1) + dudy
1402     f_Row(3,atom1) = f_Row(3,atom1) + dudz
1403    
1404     f_Col(1,atom2) = f_Col(1,atom2) - dudx
1405     f_Col(2,atom2) = f_Col(2,atom2) - dudy
1406     f_Col(3,atom2) = f_Col(3,atom2) - dudz
1407     #else
1408     f(1,atom1) = f(1,atom1) + dudx
1409     f(2,atom1) = f(2,atom1) + dudy
1410     f(3,atom1) = f(3,atom1) + dudz
1411    
1412     f(1,atom2) = f(1,atom2) - dudx
1413     f(2,atom2) = f(2,atom2) - dudy
1414     f(3,atom2) = f(3,atom2) - dudz
1415     #endif
1416    
1417     #ifdef IS_MPI
1418     id1 = AtomRowToGlobal(atom1)
1419     id2 = AtomColToGlobal(atom2)
1420     #else
1421     id1 = atom1
1422     id2 = atom2
1423     #endif
1424    
1425     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1426    
1427     fpair(1) = fpair(1) + dudx
1428     fpair(2) = fpair(2) + dudy
1429     fpair(3) = fpair(3) + dudz
1430    
1431     endif
1432    
1433     end if
1434     return
1435     end subroutine rf_correct_forces
1436    
1437 gezelter 411 end module electrostatic_module