UseTheForce/DarkSide/electrostatic.F90

!!
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!

module electrostatic_module

  use force_globals
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none

  PRIVATE


#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"
#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"


  !! these prefactors convert the multipole interactions into kcal / mol
  !! all were computed assuming distances are measured in angstroms
  !! Charge-Charge, assuming charges are measured in electrons
  real(kind=dp), parameter :: pre11 = 332.0637778_dp
  !! Charge-Dipole, assuming charges are measured in electrons, and
  !! dipoles are measured in debyes
  real(kind=dp), parameter :: pre12 = 69.13373_dp
  !! Dipole-Dipole, assuming dipoles are measured in debyes
  real(kind=dp), parameter :: pre22 = 14.39325_dp
  !! Charge-Quadrupole, assuming charges are measured in electrons, and
  !! quadrupoles are measured in 10^-26 esu cm^2
  !! This unit is also known affectionately as an esu centi-barn.
  real(kind=dp), parameter :: pre14 = 69.13373_dp

  !! variables to handle different summation methods for long-range electrostatics:
  integer, save :: summationMethod = NONE
  logical, save :: summationMethodChecked = .false.
  real(kind=DP), save :: defaultCutoff = 0.0_DP
  real(kind=DP), save :: defaultCutoff2 = 0.0_DP
  logical, save :: haveDefaultCutoff = .false.
  real(kind=DP), save :: dampingAlpha = 0.0_DP
  logical, save :: haveDampingAlpha = .false.
  real(kind=DP), save :: dielectric = 1.0_DP
  logical, save :: haveDielectric = .false.
  real(kind=DP), save :: constERFC = 0.0_DP
  real(kind=DP), save :: constEXP = 0.0_DP
  logical, save :: haveDWAconstants = .false.
  real(kind=dp), save :: rcuti = 0.0_DP
  real(kind=dp), save :: rcuti2 = 0.0_DP
  real(kind=dp), save :: rcuti3 = 0.0_DP
  real(kind=dp), save :: rcuti4 = 0.0_DP
  real(kind=dp), save :: alphaPi = 0.0_DP
  real(kind=dp), save :: invRootPi = 0.0_DP
  real(kind=dp), save :: rrf = 1.0_DP
  real(kind=dp), save :: rt = 1.0_DP
  real(kind=dp), save :: rrfsq = 1.0_DP
  real(kind=dp), save :: preRF = 0.0_DP
  logical, save :: preRFCalculated = .false.
 
#ifdef __IFC
! error function for ifc version > 7.
  double precision, external :: derfc
#endif
  
  public :: setElectrostaticSummationMethod
  public :: setElectrostaticCutoffRadius
  public :: setDampedWolfAlpha
  public :: setReactionFieldDielectric
  public :: setReactionFieldPrefactor
  public :: newElectrostaticType
  public :: setCharge
  public :: setDipoleMoment
  public :: setSplitDipoleDistance
  public :: setQuadrupoleMoments
  public :: doElectrostaticPair
  public :: getCharge
  public :: getDipoleMoment
  public :: pre22
  public :: destroyElectrostaticTypes
  public :: accumulate_rf
  public :: accumulate_self_rf
  public :: reaction_field_final
  public :: rf_correct_forces

  type :: Electrostatic
     integer :: c_ident
     logical :: is_Charge = .false.
     logical :: is_Dipole = .false.
     logical :: is_SplitDipole = .false.
     logical :: is_Quadrupole = .false.
     logical :: is_Tap = .false.
     real(kind=DP) :: charge = 0.0_DP
     real(kind=DP) :: dipole_moment = 0.0_DP
     real(kind=DP) :: split_dipole_distance = 0.0_DP
     real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP
  end type Electrostatic

  type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap

contains

  subroutine setElectrostaticSummationMethod(the_ESM)
    integer, intent(in) :: the_ESM    

    if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
       call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
    endif

    summationMethod = the_ESM

  end subroutine setElectrostaticSummationMethod

  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw) 
    real(kind=dp), intent(in) :: thisRcut
    real(kind=dp), intent(in) :: thisRsw
    defaultCutoff = thisRcut
    rrf = defaultCutoff
    rt = thisRsw
    haveDefaultCutoff = .true.
  end subroutine setElectrostaticCutoffRadius

  subroutine setDampedWolfAlpha(thisAlpha)
    real(kind=dp), intent(in) :: thisAlpha
    dampingAlpha = thisAlpha
    haveDampingAlpha = .true.
  end subroutine setDampedWolfAlpha
  
  subroutine setReactionFieldDielectric(thisDielectric)
    real(kind=dp), intent(in) :: thisDielectric
    dielectric = thisDielectric
    haveDielectric = .true.
  end subroutine setReactionFieldDielectric

  subroutine setReactionFieldPrefactor
    if (haveDefaultCutoff .and. haveDielectric) then
       defaultCutoff2 = defaultCutoff*defaultCutoff
       preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
       preRFCalculated = .true.
    else if (.not.haveDefaultCutoff) then
       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
    else 
       call handleError("setReactionFieldPrefactor", "Dielectric not set")
    endif
  end subroutine setReactionFieldPrefactor

  subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
       is_SplitDipole, is_Quadrupole, is_Tap, status)

    integer, intent(in) :: c_ident
    logical, intent(in) :: is_Charge
    logical, intent(in) :: is_Dipole
    logical, intent(in) :: is_SplitDipole
    logical, intent(in) :: is_Quadrupole
    logical, intent(in) :: is_Tap
    integer, intent(out) :: status
    integer :: nAtypes, myATID, i, j

    status = 0
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    !! Be simple-minded and assume that we need an ElectrostaticMap that
    !! is the same size as the total number of atom types

    if (.not.allocated(ElectrostaticMap)) then

       nAtypes = getSize(atypes)

       if (nAtypes == 0) then
          status = -1
          return
       end if

       if (.not. allocated(ElectrostaticMap)) then
          allocate(ElectrostaticMap(nAtypes))
       endif

    end if

    if (myATID .gt. size(ElectrostaticMap)) then
       status = -1
       return
    endif

    ! set the values for ElectrostaticMap for this atom type:

    ElectrostaticMap(myATID)%c_ident = c_ident
    ElectrostaticMap(myATID)%is_Charge = is_Charge
    ElectrostaticMap(myATID)%is_Dipole = is_Dipole
    ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
    ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
    ElectrostaticMap(myATID)%is_Tap = is_Tap

  end subroutine newElectrostaticType

  subroutine setCharge(c_ident, charge, status)
    integer, intent(in) :: c_ident
    real(kind=dp), intent(in) :: charge
    integer, intent(out) :: status
    integer :: myATID

    status = 0
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!")
       status = -1
       return
    end if

    if (myATID .gt. size(ElectrostaticMap)) then
       call handleError("electrostatic", "ElectrostaticMap was found to be too small during setCharge!")
       status = -1
       return
    endif

    if (.not.ElectrostaticMap(myATID)%is_Charge) then
       call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
       status = -1
       return
    endif

    ElectrostaticMap(myATID)%charge = charge
  end subroutine setCharge

  subroutine setDipoleMoment(c_ident, dipole_moment, status)
    integer, intent(in) :: c_ident
    real(kind=dp), intent(in) :: dipole_moment
    integer, intent(out) :: status
    integer :: myATID

    status = 0
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!")
       status = -1
       return
    end if

    if (myATID .gt. size(ElectrostaticMap)) then
       call handleError("electrostatic", "ElectrostaticMap was found to be too small during setDipoleMoment!")
       status = -1
       return
    endif

    if (.not.ElectrostaticMap(myATID)%is_Dipole) then
       call handleError("electrostatic", "Attempt to setDipoleMoment of an atom type that is not a dipole!")
       status = -1
       return
    endif

    ElectrostaticMap(myATID)%dipole_moment = dipole_moment
  end subroutine setDipoleMoment

  subroutine setSplitDipoleDistance(c_ident, split_dipole_distance, status)
    integer, intent(in) :: c_ident
    real(kind=dp), intent(in) :: split_dipole_distance
    integer, intent(out) :: status
    integer :: myATID

    status = 0
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!")
       status = -1
       return
    end if

    if (myATID .gt. size(ElectrostaticMap)) then
       call handleError("electrostatic", "ElectrostaticMap was found to be too small during setSplitDipoleDistance!")
       status = -1
       return
    endif

    if (.not.ElectrostaticMap(myATID)%is_SplitDipole) then
       call handleError("electrostatic", "Attempt to setSplitDipoleDistance of an atom type that is not a splitDipole!")
       status = -1
       return
    endif

    ElectrostaticMap(myATID)%split_dipole_distance = split_dipole_distance
  end subroutine setSplitDipoleDistance

  subroutine setQuadrupoleMoments(c_ident, quadrupole_moments, status)
    integer, intent(in) :: c_ident
    real(kind=dp), intent(in), dimension(3) :: quadrupole_moments
    integer, intent(out) :: status
    integer :: myATID, i, j

    status = 0
    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
       status = -1
       return
    end if

    if (myATID .gt. size(ElectrostaticMap)) then
       call handleError("electrostatic", "ElectrostaticMap was found to be too small during setQuadrupoleMoments!")
       status = -1
       return
    endif

    if (.not.ElectrostaticMap(myATID)%is_Quadrupole) then
       call handleError("electrostatic", "Attempt to setQuadrupoleMoments of an atom type that is not a quadrupole!")
       status = -1
       return
    endif

    do i = 1, 3
       ElectrostaticMap(myATID)%quadrupole_moments(i) = &
            quadrupole_moments(i)
    enddo

  end subroutine setQuadrupoleMoments


  function getCharge(atid) result (c)
    integer, intent(in) :: atid
    integer :: localError
    real(kind=dp) :: c

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
       return
    end if

    if (.not.ElectrostaticMap(atid)%is_Charge) then
       call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
       return
    endif

    c = ElectrostaticMap(atid)%charge
  end function getCharge

  function getDipoleMoment(atid) result (dm)
    integer, intent(in) :: atid
    integer :: localError
    real(kind=dp) :: dm

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
       return
    end if

    if (.not.ElectrostaticMap(atid)%is_Dipole) then
       call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
       return
    endif

    dm = ElectrostaticMap(atid)%dipole_moment
  end function getDipoleMoment

  subroutine checkSummationMethod()

    if (.not.haveDefaultCutoff) then
       call handleError("checkSummationMethod", "no Default Cutoff set!")
    endif

    rcuti = 1.0d0 / defaultCutoff
    rcuti2 = rcuti*rcuti
    rcuti3 = rcuti2*rcuti
    rcuti4 = rcuti2*rcuti2

    if (summationMethod .eq. DAMPED_WOLF) then
       if (.not.haveDWAconstants) then
          
          if (.not.haveDampingAlpha) then
             call handleError("checkSummationMethod", "no Damping Alpha set!")
          endif
          
          if (.not.haveDefaultCutoff) then
             call handleError("checkSummationMethod", "no Default Cutoff set!")
          endif

          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
          constERFC = derfc(dampingAlpha*defaultCutoff)
          invRootPi = 0.56418958354775628695d0
          alphaPi = 2*dampingAlpha*invRootPi
  
          haveDWAconstants = .true.
       endif
    endif

    if (summationMethod .eq. REACTION_FIELD) then
       if (.not.haveDielectric) then
          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
       endif
    endif

    summationMethodChecked = .true.
  end subroutine checkSummationMethod


  subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
       vpair, fpair, pot, eFrame, f, t, do_pot)

    logical, intent(in) :: do_pot

    integer, intent(in) :: atom1, atom2
    integer :: localError

    real(kind=dp), intent(in) :: rij, r2, sw
    real(kind=dp), intent(in), dimension(3) :: d
    real(kind=dp), intent(inout) :: vpair
    real(kind=dp), intent(inout), dimension(3) :: fpair

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(9,nLocal) :: eFrame
    real( kind = dp ), dimension(3,nLocal) :: f
    real( kind = dp ), dimension(3,nLocal) :: t

    real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
    real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
    real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
    real (kind = dp), dimension(3) :: dudux_j, duduy_j, duduz_j

    logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
    logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
    logical :: i_is_Tap, j_is_Tap
    integer :: me1, me2, id1, id2
    real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
    real (kind=dp) :: qxx_i, qyy_i, qzz_i
    real (kind=dp) :: qxx_j, qyy_j, qzz_j
    real (kind=dp) :: cx_i, cy_i, cz_i
    real (kind=dp) :: cx_j, cy_j, cz_j
    real (kind=dp) :: cx2, cy2, cz2
    real (kind=dp) :: ct_i, ct_j, ct_ij, a1
    real (kind=dp) :: riji, ri, ri2, ri3, ri4
    real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2 
    real (kind=dp) :: xhat, yhat, zhat
    real (kind=dp) :: dudx, dudy, dudz
    real (kind=dp) :: scale, sc2, bigR
    real (kind=dp) :: varERFC, varEXP
    real (kind=dp) :: limScale

    if (.not.allocated(ElectrostaticMap)) then
       call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
       return
    end if

    if (.not.summationMethodChecked) then
       call checkSummationMethod()
       
    endif


#ifdef IS_MPI
    me1 = atid_Row(atom1)
    me2 = atid_Col(atom2)
#else
    me1 = atid(atom1)
    me2 = atid(atom2)
#endif

    !! some variables we'll need independent of electrostatic type:

    riji = 1.0d0 / rij
   
    xhat = d(1) * riji
    yhat = d(2) * riji
    zhat = d(3) * riji

    !! logicals
    i_is_Charge = ElectrostaticMap(me1)%is_Charge
    i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
    i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
    i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
    i_is_Tap = ElectrostaticMap(me1)%is_Tap

    j_is_Charge = ElectrostaticMap(me2)%is_Charge
    j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
    j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
    j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
    j_is_Tap = ElectrostaticMap(me2)%is_Tap

    if (i_is_Charge) then
       q_i = ElectrostaticMap(me1)%charge      
    endif

    if (i_is_Dipole) then
       mu_i = ElectrostaticMap(me1)%dipole_moment
#ifdef IS_MPI
       uz_i(1) = eFrame_Row(3,atom1)
       uz_i(2) = eFrame_Row(6,atom1)
       uz_i(3) = eFrame_Row(9,atom1)
#else
       uz_i(1) = eFrame(3,atom1)
       uz_i(2) = eFrame(6,atom1)
       uz_i(3) = eFrame(9,atom1)
#endif
       ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat

       if (i_is_SplitDipole) then
          d_i = ElectrostaticMap(me1)%split_dipole_distance
       endif

    endif

    if (i_is_Quadrupole) then
       qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1)
       qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2)
       qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3)
#ifdef IS_MPI
       ux_i(1) = eFrame_Row(1,atom1)
       ux_i(2) = eFrame_Row(4,atom1)
       ux_i(3) = eFrame_Row(7,atom1)
       uy_i(1) = eFrame_Row(2,atom1)
       uy_i(2) = eFrame_Row(5,atom1)
       uy_i(3) = eFrame_Row(8,atom1)
       uz_i(1) = eFrame_Row(3,atom1)
       uz_i(2) = eFrame_Row(6,atom1)
       uz_i(3) = eFrame_Row(9,atom1)
#else
       ux_i(1) = eFrame(1,atom1)
       ux_i(2) = eFrame(4,atom1)
       ux_i(3) = eFrame(7,atom1)
       uy_i(1) = eFrame(2,atom1)
       uy_i(2) = eFrame(5,atom1)
       uy_i(3) = eFrame(8,atom1)
       uz_i(1) = eFrame(3,atom1)
       uz_i(2) = eFrame(6,atom1)
       uz_i(3) = eFrame(9,atom1)
#endif
       cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
       cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
       cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
    endif

    if (j_is_Charge) then
       q_j = ElectrostaticMap(me2)%charge      
    endif

    if (j_is_Dipole) then
       mu_j = ElectrostaticMap(me2)%dipole_moment
#ifdef IS_MPI
       uz_j(1) = eFrame_Col(3,atom2)
       uz_j(2) = eFrame_Col(6,atom2)
       uz_j(3) = eFrame_Col(9,atom2)
#else
       uz_j(1) = eFrame(3,atom2)
       uz_j(2) = eFrame(6,atom2)
       uz_j(3) = eFrame(9,atom2)
#endif
       ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat

       if (j_is_SplitDipole) then
          d_j = ElectrostaticMap(me2)%split_dipole_distance
       endif
    endif

    if (j_is_Quadrupole) then
       qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1)
       qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2)
       qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3)
#ifdef IS_MPI
       ux_j(1) = eFrame_Col(1,atom2)
       ux_j(2) = eFrame_Col(4,atom2)
       ux_j(3) = eFrame_Col(7,atom2)
       uy_j(1) = eFrame_Col(2,atom2)
       uy_j(2) = eFrame_Col(5,atom2)
       uy_j(3) = eFrame_Col(8,atom2)
       uz_j(1) = eFrame_Col(3,atom2)
       uz_j(2) = eFrame_Col(6,atom2)
       uz_j(3) = eFrame_Col(9,atom2)
#else
       ux_j(1) = eFrame(1,atom2)
       ux_j(2) = eFrame(4,atom2)
       ux_j(3) = eFrame(7,atom2)
       uy_j(1) = eFrame(2,atom2)
       uy_j(2) = eFrame(5,atom2)
       uy_j(3) = eFrame(8,atom2)
       uz_j(1) = eFrame(3,atom2)
       uz_j(2) = eFrame(6,atom2)
       uz_j(3) = eFrame(9,atom2)
#endif
       cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
       cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
       cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
    endif
   
    epot = 0.0_dp
    dudx = 0.0_dp
    dudy = 0.0_dp
    dudz = 0.0_dp

    dudux_i = 0.0_dp
    duduy_i = 0.0_dp
    duduz_i = 0.0_dp

    dudux_j = 0.0_dp
    duduy_j = 0.0_dp
    duduz_j = 0.0_dp

    if (i_is_Charge) then

       if (j_is_Charge) then

          if (summationMethod .eq. UNDAMPED_WOLF) then

             vterm = pre11 * q_i * q_j * (riji - rcuti)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
             
             dudx = dudx + dudr * d(1)
             dudy = dudy + dudr * d(2)
             dudz = dudz + dudr * d(3)

          elseif (summationMethod .eq. DAMPED_WOLF) then

             varERFC = derfc(dampingAlpha*rij)
             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
                                                  + alphaPi*varEXP) &
                                                 - (constERFC*rcuti2 &
                                                    + alphaPi*constEXP)) )
             
             dudx = dudx + dudr * d(1)
             dudy = dudy + dudr * d(2)
             dudz = dudz + dudr * d(3)

          else

             vterm = pre11 * q_i * q_j * riji
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudr  = - sw * vterm * riji
             
             dudx = dudx + dudr * xhat
             dudy = dudy + dudr * yhat
             dudz = dudz + dudr * zhat

          endif

       endif

       if (j_is_Dipole) then

          pref = pre12 * q_i * mu_j

          if (summationMethod .eq. UNDAMPED_WOLF) then
             ri2 = riji * riji
             ri3 = ri2 * riji

             pref = pre12 * q_i * mu_j
             vterm = - pref * ct_j * (ri2 - rcuti2)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             !! this has a + sign in the () because the rij vector is
             !! r_j - r_i and the charge-dipole potential takes the origin
             !! as the point dipole, which is atom j in this case.
             
             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
             
             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )

          else
             if (j_is_SplitDipole) then
                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
                ri = 1.0_dp / BigR
                scale = rij * ri
             else
                ri = riji
                scale = 1.0_dp
             endif
             
             ri2 = ri * ri
             ri3 = ri2 * ri
             sc2 = scale * scale

             pref = pre12 * q_i * mu_j
             vterm = - pref * ct_j * ri2 * scale
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             !! this has a + sign in the () because the rij vector is
             !! r_j - r_i and the charge-dipole potential takes the origin
             !! as the point dipole, which is atom j in this case.
             
             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
             
             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale

          endif
       endif

       if (j_is_Quadrupole) then
          ri2 = riji * riji
          ri3 = ri2 * riji
          ri4 = ri2 * ri2
          cx2 = cx_j * cx_j
          cy2 = cy_j * cy_j
          cz2 = cz_j * cz_j

          if (summationMethod .eq. UNDAMPED_WOLF) then
             pref =  pre14 * q_i / 3.0_dp
             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
             vpair = vpair + ( vterm1 - vterm2 )
             epot = epot + sw*( vterm1 - vterm2 )
             
             dudx = dudx - (5.0_dp * &
                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) ) 
             dudy = dudy - (5.0_dp * &
                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
             dudz = dudz - (5.0_dp * &
                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) ) 
             
             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
             
             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
             
             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
         
          else
             pref =  pre14 * q_i / 3.0_dp
             vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_j * (3.0_dp*cz2 - 1.0_dp))
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
                  qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
                  qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
                  qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) ) 
             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
                  qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
                  qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
                  qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
                  qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
                  qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
                  qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) ) 
             
             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
             
             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
             
             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
          
          endif
       endif
    endif

    if (i_is_Dipole) then 

       if (j_is_Charge) then
          
          pref = pre12 * q_j * mu_i
          
          if (summationMethod .eq. UNDAMPED_WOLF) then
             ri2 = riji * riji
             ri3 = ri2 * riji

             pref = pre12 * q_j * mu_i
             vterm = pref * ct_i * (ri2 - rcuti2)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             !! this has a + sign in the () because the rij vector is
             !! r_j - r_i and the charge-dipole potential takes the origin
             !! as the point dipole, which is atom j in this case.
             
             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
             
             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )

          else
             if (i_is_SplitDipole) then
                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
                ri = 1.0_dp / BigR
                scale = rij * ri
             else
                ri = riji
                scale = 1.0_dp
             endif
             
             ri2 = ri * ri
             ri3 = ri2 * ri
             sc2 = scale * scale

             pref = pre12 * q_j * mu_i
             vterm = pref * ct_i * ri2 * scale
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
             
             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
          endif
       endif
       
       if (j_is_Dipole) then

          if (summationMethod .eq. UNDAMPED_WOLF) then
             ri2 = riji * riji
             ri3 = ri2 * riji
             ri4 = ri2 * ri2

             pref = pre22 * mu_i * mu_j
             vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             a1 = 5.0d0 * ct_i * ct_j - ct_ij
             
             dudx = dudx + sw*pref*3.0d0*ri4 &
                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
                         - sw*pref*3.0d0*rcuti4 &
                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
             dudy = dudy + sw*pref*3.0d0*ri4 &
                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
                         - sw*pref*3.0d0*rcuti4 &
                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
             dudz = dudz + sw*pref*3.0d0*ri4 &
                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
                         - sw*pref*3.0d0*rcuti4 &
                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
             
             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
                  - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
                  - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
                  - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
                  - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))

          else
             if (i_is_SplitDipole) then
                if (j_is_SplitDipole) then
                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
                else
                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
                endif
                ri = 1.0_dp / BigR
                scale = rij * ri                
             else
                if (j_is_SplitDipole) then
                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
                   ri = 1.0_dp / BigR
                   scale = rij * ri                             
                else                
                   ri = riji
                   scale = 1.0_dp
                endif
             endif
             
             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
             
             ri2 = ri * ri
             ri3 = ri2 * ri
             ri4 = ri2 * ri2
             sc2 = scale * scale
             
             pref = pre22 * mu_i * mu_j
             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
             
             dudx = dudx + sw*pref*3.0d0*ri4*scale &
                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
             dudy = dudy + sw*pref*3.0d0*ri4*scale &
                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
             dudz = dudz + sw*pref*3.0d0*ri4*scale &
                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
             
             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
             
             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
          endif
       endif
    endif

    if (i_is_Quadrupole) then
       if (j_is_Charge) then

          ri2 = riji * riji
          ri3 = ri2 * riji
          ri4 = ri2 * ri2
          cx2 = cx_i * cx_i
          cy2 = cy_i * cy_i
          cz2 = cz_i * cz_i

          if (summationMethod .eq. UNDAMPED_WOLF) then
             pref = pre14 * q_j / 3.0_dp
             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
             vpair = vpair + ( vterm1 - vterm2 )
             epot = epot + sw*( vterm1 - vterm2 )
             
             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) ) 
             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) ) 
             
             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
             
             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
             
             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))

          else
             pref = pre14 * q_j / 3.0_dp
             vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
                  qzz_i * (3.0_dp*cz2 - 1.0_dp))
             vpair = vpair + vterm
             epot = epot + sw*vterm
             
             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
                  qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
                  qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
                  qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) ) 
             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
                  qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
                  qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
                  qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
                  qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
                  qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
                  qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) ) 
             
             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
             
             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
             
             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
          endif
       endif
    endif


    if (do_pot) then
#ifdef IS_MPI 
       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
#else
       pot = pot + epot
#endif
    endif

#ifdef IS_MPI
    f_Row(1,atom1) = f_Row(1,atom1) + dudx
    f_Row(2,atom1) = f_Row(2,atom1) + dudy
    f_Row(3,atom1) = f_Row(3,atom1) + dudz

    f_Col(1,atom2) = f_Col(1,atom2) - dudx
    f_Col(2,atom2) = f_Col(2,atom2) - dudy
    f_Col(3,atom2) = f_Col(3,atom2) - dudz

    if (i_is_Dipole .or. i_is_Quadrupole) then
       t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
       t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
       t_Row(3,atom1)=t_Row(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
    endif
    if (i_is_Quadrupole) then
       t_Row(1,atom1)=t_Row(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
       t_Row(2,atom1)=t_Row(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
       t_Row(3,atom1)=t_Row(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)

       t_Row(1,atom1)=t_Row(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
       t_Row(2,atom1)=t_Row(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
       t_Row(3,atom1)=t_Row(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
    endif

    if (j_is_Dipole .or. j_is_Quadrupole) then
       t_Col(1,atom2)=t_Col(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
       t_Col(2,atom2)=t_Col(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
       t_Col(3,atom2)=t_Col(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
    endif
    if (j_is_Quadrupole) then
       t_Col(1,atom2)=t_Col(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
       t_Col(2,atom2)=t_Col(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
       t_Col(3,atom2)=t_Col(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)

       t_Col(1,atom2)=t_Col(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
       t_Col(2,atom2)=t_Col(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
       t_Col(3,atom2)=t_Col(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
    endif

#else
    f(1,atom1) = f(1,atom1) + dudx
    f(2,atom1) = f(2,atom1) + dudy
    f(3,atom1) = f(3,atom1) + dudz

    f(1,atom2) = f(1,atom2) - dudx
    f(2,atom2) = f(2,atom2) - dudy
    f(3,atom2) = f(3,atom2) - dudz

    if (i_is_Dipole .or. i_is_Quadrupole) then
       t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
       t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
       t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
    endif
    if (i_is_Quadrupole) then
       t(1,atom1)=t(1,atom1) - ux_i(2)*dudux_i(3) + ux_i(3)*dudux_i(2)
       t(2,atom1)=t(2,atom1) - ux_i(3)*dudux_i(1) + ux_i(1)*dudux_i(3)
       t(3,atom1)=t(3,atom1) - ux_i(1)*dudux_i(2) + ux_i(2)*dudux_i(1)

       t(1,atom1)=t(1,atom1) - uy_i(2)*duduy_i(3) + uy_i(3)*duduy_i(2)
       t(2,atom1)=t(2,atom1) - uy_i(3)*duduy_i(1) + uy_i(1)*duduy_i(3)
       t(3,atom1)=t(3,atom1) - uy_i(1)*duduy_i(2) + uy_i(2)*duduy_i(1)
    endif

    if (j_is_Dipole .or. j_is_Quadrupole) then
       t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
       t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
       t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
    endif
    if (j_is_Quadrupole) then
       t(1,atom2)=t(1,atom2) - ux_j(2)*dudux_j(3) + ux_j(3)*dudux_j(2)
       t(2,atom2)=t(2,atom2) - ux_j(3)*dudux_j(1) + ux_j(1)*dudux_j(3)
       t(3,atom2)=t(3,atom2) - ux_j(1)*dudux_j(2) + ux_j(2)*dudux_j(1)

       t(1,atom2)=t(1,atom2) - uy_j(2)*duduy_j(3) + uy_j(3)*duduy_j(2)
       t(2,atom2)=t(2,atom2) - uy_j(3)*duduy_j(1) + uy_j(1)*duduy_j(3)
       t(3,atom2)=t(3,atom2) - uy_j(1)*duduy_j(2) + uy_j(2)*duduy_j(1)
    endif

#endif

#ifdef IS_MPI
    id1 = AtomRowToGlobal(atom1)
    id2 = AtomColToGlobal(atom2)
#else
    id1 = atom1
    id2 = atom2
#endif

    if (molMembershipList(id1) .ne. molMembershipList(id2)) then

       fpair(1) = fpair(1) + dudx
       fpair(2) = fpair(2) + dudy
       fpair(3) = fpair(3) + dudz

    endif

    return
  end subroutine doElectrostaticPair

  !! calculates the switching functions and their derivatives for a given
  subroutine calc_switch(r, mu, scale, dscale)

    real (kind=dp), intent(in) :: r, mu
    real (kind=dp), intent(inout) :: scale, dscale
    real (kind=dp) :: rl, ru, mulow, minRatio, temp, scaleVal

    ! distances must be in angstroms
    rl = 2.75d0
    ru = 3.75d0
    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
    minRatio = mulow / (mu*mu)
    scaleVal = 1.0d0 - minRatio
    
    if (r.lt.rl) then
       scale = minRatio
       dscale = 0.0d0
    elseif (r.gt.ru) then
       scale = 1.0d0
       dscale = 0.0d0
    else
       scale = 1.0d0 - scaleVal*((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) &
                        / ((ru - rl)**3)
       dscale = -scaleVal * 6.0d0 * (r-ru)*(r-rl)/((ru - rl)**3)    
    endif
        
    return
  end subroutine calc_switch

  subroutine destroyElectrostaticTypes()

    if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)

  end subroutine destroyElectrostaticTypes

  subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)

    integer, intent(in) :: atom1, atom2
    real (kind = dp), intent(in) :: rij
    real (kind = dp), dimension(9,nLocal) :: eFrame

    integer :: me1, me2
    real (kind = dp), intent(in) :: taper
    real (kind = dp):: mu1, mu2
    real (kind = dp), dimension(3) :: ul1
    real (kind = dp), dimension(3) :: ul2   

    integer :: localError

#ifdef IS_MPI
    me1 = atid_Row(atom1)
    ul1(1) = eFrame_Row(3,atom1)
    ul1(2) = eFrame_Row(6,atom1)
    ul1(3) = eFrame_Row(9,atom1)

    me2 = atid_Col(atom2)
    ul2(1) = eFrame_Col(3,atom2)
    ul2(2) = eFrame_Col(6,atom2)
    ul2(3) = eFrame_Col(9,atom2)
#else
    me1 = atid(atom1)
    ul1(1) = eFrame(3,atom1)
    ul1(2) = eFrame(6,atom1)
    ul1(3) = eFrame(9,atom1)

    me2 = atid(atom2)
    ul2(1) = eFrame(3,atom2)
    ul2(2) = eFrame(6,atom2)
    ul2(3) = eFrame(9,atom2)
#endif

    mu1 = getDipoleMoment(me1)
    mu2 = getDipoleMoment(me2)

#ifdef IS_MPI
    rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
    rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper

    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
#else
    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper

    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper     
#endif
    return  
  end subroutine accumulate_rf

  subroutine accumulate_self_rf(atom1, mu1, eFrame)

    integer, intent(in) :: atom1
    real(kind=dp), intent(in) :: mu1
    real(kind=dp), dimension(9,nLocal) :: eFrame

    !! should work for both MPI and non-MPI version since this is not pairwise.
    rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
    rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
    rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1

    return
  end subroutine accumulate_self_rf

  subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)

    integer, intent(in) :: a1
    real (kind=dp), intent(in) :: mu1
    real (kind=dp), intent(inout) :: rfpot
    logical, intent(in) :: do_pot
    real (kind = dp), dimension(9,nLocal) :: eFrame
    real (kind = dp), dimension(3,nLocal) :: t

    integer :: localError

    if (.not.preRFCalculated) then
       call setReactionFieldPrefactor()
    endif

    ! compute torques on dipoles:
    ! pre converts from mu in units of debye to kcal/mol

    ! The torque contribution is dipole cross reaction_field   

    t(1,a1) = t(1,a1) + preRF*mu1*(eFrame(6,a1)*rf(3,a1) - &
                                   eFrame(9,a1)*rf(2,a1))
    t(2,a1) = t(2,a1) + preRF*mu1*(eFrame(9,a1)*rf(1,a1) - &
                                   eFrame(3,a1)*rf(3,a1))
    t(3,a1) = t(3,a1) + preRF*mu1*(eFrame(3,a1)*rf(2,a1) - &
                                   eFrame(6,a1)*rf(1,a1))

    ! the potential contribution is -1/2 dipole dot reaction_field

    if (do_pot) then
       rfpot = rfpot - 0.5d0 * preRF * mu1 * &
            (rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + &
             rf(3,a1)*eFrame(9,a1))
    endif

    return
  end subroutine reaction_field_final

  subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)

    integer, intent(in) :: atom1, atom2
    real(kind=dp), dimension(3), intent(in) :: d
    real(kind=dp), intent(in) :: rij, taper
    real( kind = dp ), dimension(9,nLocal) :: eFrame
    real( kind = dp ), dimension(3,nLocal) :: f
    real( kind = dp ), dimension(3), intent(inout) :: fpair

    real (kind = dp), dimension(3) :: ul1
    real (kind = dp), dimension(3) :: ul2
    real (kind = dp) :: dtdr
    real (kind = dp) :: dudx, dudy, dudz, u1dotu2
    integer :: me1, me2, id1, id2
    real (kind = dp) :: mu1, mu2

    integer :: localError

    if (.not.preRFCalculated) then
       call setReactionFieldPrefactor()
    endif

    if (rij.le.rrf) then

       if (rij.lt.rt) then
          dtdr = 0.0d0
       else
          !         write(*,*) 'rf correct in taper region'
          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
       endif

#ifdef IS_MPI
       me1 = atid_Row(atom1)
       ul1(1) = eFrame_Row(3,atom1)
       ul1(2) = eFrame_Row(6,atom1)
       ul1(3) = eFrame_Row(9,atom1)

       me2 = atid_Col(atom2)
       ul2(1) = eFrame_Col(3,atom2)
       ul2(2) = eFrame_Col(6,atom2)
       ul2(3) = eFrame_Col(9,atom2)
#else
       me1 = atid(atom1)
       ul1(1) = eFrame(3,atom1)
       ul1(2) = eFrame(6,atom1)
       ul1(3) = eFrame(9,atom1)

       me2 = atid(atom2)
       ul2(1) = eFrame(3,atom2)
       ul2(2) = eFrame(6,atom2)
       ul2(3) = eFrame(9,atom2)
#endif

       mu1 = getDipoleMoment(me1)
       mu2 = getDipoleMoment(me2)

       u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)

       dudx = - preRF*mu1*mu2*u1dotu2*dtdr*d(1)/rij
       dudy = - preRF*mu1*mu2*u1dotu2*dtdr*d(2)/rij
       dudz = - preRF*mu1*mu2*u1dotu2*dtdr*d(3)/rij

#ifdef IS_MPI
       f_Row(1,atom1) = f_Row(1,atom1) + dudx
       f_Row(2,atom1) = f_Row(2,atom1) + dudy
       f_Row(3,atom1) = f_Row(3,atom1) + dudz

       f_Col(1,atom2) = f_Col(1,atom2) - dudx
       f_Col(2,atom2) = f_Col(2,atom2) - dudy
       f_Col(3,atom2) = f_Col(3,atom2) - dudz
#else
       f(1,atom1) = f(1,atom1) + dudx
       f(2,atom1) = f(2,atom1) + dudy
       f(3,atom1) = f(3,atom1) + dudz

       f(1,atom2) = f(1,atom2) - dudx
       f(2,atom2) = f(2,atom2) - dudy
       f(3,atom2) = f(3,atom2) - dudz
#endif

#ifdef IS_MPI
       id1 = AtomRowToGlobal(atom1)
       id2 = AtomColToGlobal(atom2)
#else
       id1 = atom1
       id2 = atom2
#endif

       if (molMembershipList(id1) .ne. molMembershipList(id2)) then

          fpair(1) = fpair(1) + dudx
          fpair(2) = fpair(2) + dudy
          fpair(3) = fpair(3) + dudz

       endif

    end if
    return
  end subroutine rf_correct_forces

end module electrostatic_module
Revision:	2381
Committed:	Tue Oct 18 15:01:42 2005 UTC (18 years, 10 months ago) by chrisfen
File size:	53872 byte(s)
Log Message:	merged reaction field with electrostatics.F90