! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X Ewald module
!X
!% This module computes the Ewald sum.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
module Ewald_module

use system_module
use units_module

use functions_module

implicit none

private

integer max_n_r, max_n_q
integer n_r, n_q
integer, allocatable:: r_list(:,:), q_list(:,:)

integer:: last_N_A = 0
integer:: n_pairs
integer, allocatable:: pair_list(:,:)

public add_madelung_matrix, add_dmadelung_matrix, add_dmadelung_matrix_dr

contains

!| subroutine add_madelung_matrix
!| calculate Madelung matrix (interactions of pt charges in 3-D periodic supercell)
!|
!|I integer N_A: number of atoms
!|I real(dp) a1(3), a2(3), a3(3): periodic supercell vectors
!|I real(dp) A(3,N_A): atomic positions
!|O real(dp) madelung_matrix(N_A,N_A): Madelung matrix
!|I logical redo_lattice: if true, redo whole calculation, rather than
!|    just summing over periodic images that had significant contributions
!|    last time
!|
!| Noam Bernstein 9/12/2001
!| From Fundamental Formulas of Physics, p. 599
!% Calculate Madelung matrix (interactions of pt charges in 3-D periodic supercell).
!% It is possible to choose between redoing the whole calculation and summing over those periodic images
!% that had significant contributions last time by setting the logical variable redo_lattice.
subroutine add_madelung_matrix(N_A, a1_in, a2_in, a3_in, A_in, madelung_matrix, redo_lattice)
    integer, intent(in) :: N_A                           !% number of atoms (input)
    real(dp), intent(in) :: a1_in(3), a2_in(3), a3_in(3) !% periodic supercell vectors (input)
    real(dp), intent(in) :: A_in(3,N_A)                  !% atomic positions (input)
    real(dp), intent(inout) :: madelung_matrix(N_A,N_A)  !% Madelung matrix  (output)
    logical, intent(in), optional :: redo_lattice        !% if true, redo whole calculation, rather than just summing over periodic images that had significant contributions last time (input)

    integer i, j, range_a(3)
    integer s_a(3), s2, s3
    real(dp) b1(3), b2(3), b3(3), h(3,3)
    real(dp) vol
    real(dp) K(3), ri(3), rj(3), t, tv(3), Kmagsq, Rmag

    real(dp) :: a1(3), a2(3), a3(3)
    real(dp), allocatable :: A(:,:)

    integer rr, qq

    real(dp) epsilon, fourepsilonsq
    real(dp) recip_space_prefac, K_dep_factor
    real(dp) phase
    real(dp) conv_crit, prev_conv_crit

    integer ni, nip1, nip2
    real(dp) b1mag, b2mag, b3mag
    real(dp) a1mag, a2mag, a3mag

    logical use_this_one
    real(dp) convf

    logical u_redo_lattice

    integer i_pair

    integer nv

    allocate(A(3,N_A))
    a1 = a1_in / Bohr
    a2 = a2_in / Bohr
    a3 = a3_in / Bohr
    A = A_in / Bohr

    u_redo_lattice = .true.
    if (present(redo_lattice)) u_redo_lattice = redo_lattice

    !convf = 1.0D-10
    convf = 1.0D-14

    h(1:3,1) = a1(1:3)
    h(1:3,2) = a2(1:3)
    h(1:3,3) = a3(1:3)
    vol = calc_volume(h)
    call cross_3 (a1, a2, b3)
    call cross_3 (a2, a3, b1)
    call cross_3 (a3, a1, b2)
    b1 = 2.0_dp*PI* b1 / vol
    b2 = 2.0_dp*PI* b2 / vol
    b3 = 2.0_dp*PI* b3 / vol

    !epsilon = 0.7_dp/nn_dist

    epsilon = 1.5_dp / ( 1.0_dp/sqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )
    !epsilon = 2.25_dp / ( 1.0_dp/sqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )

    fourepsilonsq = 4.0_dp*epsilon**2

    recip_space_prefac = 4.0_dp*PI/vol

    prev_conv_crit = 0.0_dp
    conv_crit = prev_conv_crit+1.0_dp

    b1mag = sqrt(sum(b1**2))
    b2mag = sqrt(sum(b2**2))
    b3mag = sqrt(sum(b3**2))

    nv = 0
    range_a = 0

    if (u_redo_lattice) then

        if (N_A .ne. last_N_A) then
            if (allocated(pair_list)) deallocate(pair_list)
            call assign_atom_pairs(N_A, n_pairs, pair_list)
            last_N_A = N_A
        end if

        max_n_r = 2000
        max_n_q = 2000
        if (allocated(r_list)) deallocate(r_list)
        if (allocated(q_list)) deallocate(q_list)
        allocate(r_list(3,max_n_r))
        allocate(q_list(3,max_n_q))
        n_r = 0
        n_q = 0


        ! reciprocal space
        do while (abs(conv_crit-prev_conv_crit) .gt. 1.0e-14_dp)

            !!!!!
            do ni=1,3
            range_a(ni) = range_a(ni) + 1
            nip1 = mod(ni,3)+1
            nip2 = mod(ni+1,3)+1

            s_a(ni) = -range_a(ni)
            do s2 = -range_a(nip1), range_a(nip1)
                s_a(nip1) = s2
                do s3 = -range_a(nip2), range_a(nip2)
                    s_a(nip2) = s3

                    use_this_one = .false.

                    K = s_a(1)*b1 + s_a(2)*b2 + s_a(3)*b3
                    Kmagsq = sum(K**2)
                    K_dep_factor = recip_space_prefac*exp(-Kmagsq/fourepsilonsq)/Kmagsq

                    do i_pair = 1, n_pairs

                        i = pair_list(1,i_pair)
                        j = pair_list(2,i_pair)
                        ri = A(1:3,i)
                        rj = A(1:3,j)

                        phase = sum(K*(ri-rj))
                        t = Hartree*K_dep_factor*cos(phase)
                        madelung_matrix(i,j) = madelung_matrix(i,j) + 2*t

                        if (i /= j) madelung_matrix(j,i) = madelung_matrix(j,i) + 2*t

                        if (t /= 0.0_dp .and. abs((2.0_dp*t)/madelung_matrix(i,j)) > convf) then
! workaround for buggy Cray MTA compiler
#ifndef cray_mta
                            use_this_one = .true.
#endif
                        end if

                    end do

                    if (use_this_one) then
                        n_q = n_q + 1
                        q_list(1:3,n_q) = s_a
                    end if

                    nv = nv + 1

                end do ! s3
            end do ! s2
            end do ! ni

            prev_conv_crit = conv_crit
            !NO_PARALLEL t_conv_crit = sum(abs(madelung_matrix))
            !NO_PARALLEL call global_sum(t_conv_crit, conv_crit)
            conv_crit = sum(abs(madelung_matrix))
            ! print *, "G conv ", range_a, conv_crit, prev_conv_crit

        end do ! while not conv
        ! print *, "nv ", nv


        !NO_PARALLEL t_conv_crit = sum(abs(madelung_matrix))
        !NO_PARALLEL call global_sum(t_conv_crit,prev_conv_crit)
        prev_conv_crit = sum(abs(madelung_matrix))
        conv_crit = prev_conv_crit+1.0_dp

        ! real space shift=0
        a1mag = sqrt(sum(a1**2))
        a2mag = sqrt(sum(a2**2))
        a3mag = sqrt(sum(a3**2))

        tv = 0.0_dp

        do i_pair = 1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            if (i .eq. j) cycle

            ri = A(1:3,i)
            rj = A(1:3,j)

            Rmag = sqrt(sum((rj+tv-ri)**2))
            t = Hartree * erfc(epsilon*Rmag)/Rmag
            madelung_matrix(i,j) = madelung_matrix(i,j) + t
            if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + t
        end do


        prev_conv_crit = conv_crit
        !NO_PARALLEL t_conv_crit = sum(dabs(madelung_matrix))
        !NO_PARALLEL call global_sum(t_conv_crit, conv_crit)
        conv_crit = sum(abs(madelung_matrix))
        ! print *, "R conv ", range, conv_crit, prev_conv_crit

        nv = 0
        range_a = 0

        ! real space shift /= 0
        do while (abs(conv_crit-prev_conv_crit) .gt. 1.0e-14_dp)

            ! amag_a(1) = a1mag*(range_a(1)+1)
            ! amag_a(2) = a2mag*(range_a(2)+1)
            ! amag_a(3) = a3mag*(range_a(3)+1)
            ! next_increment = minloc(amag_a)
            ! ni = next_increment(1)

            !!!!!
            do ni=1,3

            range_a(ni) = range_a(ni) + 1
            nip1 = mod(ni,3)+1
            nip2 = mod(ni+1,3)+1

            s_a(ni) = -range_a(ni)
            do s2 = -range_a(nip1), range_a(nip1)
                s_a(nip1) = s2
                do s3 = -range_a(nip2), range_a(nip2)
                    s_a(nip2) = s3

                    use_this_one = .false.

                    tv = s_a(1)*a1 + s_a(2)*a2 + s_a(3)*a3

                    do i_pair = 1, n_pairs
                        i = pair_list(1,i_pair)
                        j = pair_list(2,i_pair)
                        ri = A(1:3,i)
                        rj = A(1:3,j)

                        Rmag = sqrt(sum((rj+tv-ri)**2))
                        t = Hartree * erfc(epsilon*Rmag)/Rmag
                        madelung_matrix(i,j) = madelung_matrix(i,j) + t
                        if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + t
                        if (t .ne. 0 .and. abs(t/madelung_matrix(i,j)) .gt. convf) then
                            use_this_one = .true.
                        end if

                    end do

                    !NO_PARALLEL call global_or(use_this_one, g_use_this_one)
                    !NO_PARALLEL use_this_one = g_use_this_one

                    if (use_this_one) then
                        n_r = n_r + 1
                        r_list(1:3,n_r) = s_a
                    end if

                    nv = nv + 1
                end do
            end do

            s_a(ni) = range_a(ni)
            do s2 = -range_a(nip1), range_a(nip1)
                s_a(nip1) = s2
                do s3 = -range_a(nip2), range_a(nip2)
                    s_a(nip2) = s3

                    use_this_one = .false.

                    tv = s_a(1)*a1 + s_a(2)*a2 + s_a(3)*a3
                    do i_pair = 1, n_pairs
                        i = pair_list(1,i_pair)
                        j = pair_list(2,i_pair)
                        ri = A(1:3,i)
                        rj = A(1:3,j)

                        Rmag = sqrt(sum((rj+tv-ri)**2))
                        t = Hartree * erfc(epsilon*Rmag)/Rmag
                        madelung_matrix(i,j) = madelung_matrix(i,j) + t
                        if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + t

                        if (t .ne. 0 .and. abs(t/madelung_matrix(i,j)) .gt. convf) then
                            use_this_one = .true.
                        end if

                    end do

                    !NO_PARALLEL call global_or(use_this_one, g_use_this_one)
                    !NO_PARALLEL use_this_one = g_use_this_one

                    if (use_this_one) then
                        n_r = n_r + 1
                        r_list(1:3,n_r) = s_a
                    end if

                    nv = nv + 1
                end do
            end do

            !!!!!
            end do ! ni


            prev_conv_crit = conv_crit
            !NO_PARALLEL t_conv_crit = sum(abs(madelung_matrix))
            !NO_PARALLEL call global_sum(t_conv_crit, conv_crit)
            conv_crit = sum(abs(madelung_matrix))
            ! print *, "R conv ", range, conv_crit, prev_conv_crit

        end do ! while conv_crit

        ! print *, "nv ", nv
        ! print *, "n_r, n_q ", n_r, n_q
    else ! don't redo lattice

        do qq=1, n_q

            s_a = q_list(1:3,qq)

            K = s_a(1)*b1 + s_a(2)*b2 + s_a(3)*b3
            Kmagsq = sum(K**2)
            K_dep_factor = recip_space_prefac*exp(-Kmagsq/fourepsilonsq)/Kmagsq
            do i_pair = 1, n_pairs
                i = pair_list(1,i_pair)
                j = pair_list(2,i_pair)
                ri = A(1:3,i)
                rj = A(1:3,j)

                phase = sum(K*(ri-rj))
                t = Hartree * K_dep_factor*cos(phase)
                madelung_matrix(i,j) = madelung_matrix(i,j) + 2*t
                if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + 2*t
            end do

        end do

        tv = 0.0_dp
        do i_pair = 1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            if (i .eq. j) cycle
            ri = A(1:3,i)
            rj = A(1:3,j)

            Rmag = sqrt(sum((rj+tv-ri)**2))
            t = Hartree * erfc(epsilon*Rmag)/Rmag
            madelung_matrix(i,j) = madelung_matrix(i,j) + t
            if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + t
        end do

        do rr=1, n_r
            s_a = r_list(1:3,rr)

            tv = s_a(1)*a1 + s_a(2)*a2 + s_a(3)*a3
            do i_pair = 1, n_pairs
                i = pair_list(1,i_pair)
                j = pair_list(2,i_pair)
                ri = A(1:3,i)
                rj = A(1:3,j)

                Rmag = sqrt(sum((rj+tv-ri)**2))
                t = Hartree * erfc(epsilon*Rmag)/Rmag
                madelung_matrix(i,j) = madelung_matrix(i,j) + t
                if (i .ne. j) madelung_matrix(j,i) = madelung_matrix(j,i) + t
            end do
        end do

    end if

    conv_crit = sum(abs(madelung_matrix))
    ! print *, "overall conv ", conv_crit

    do i=1, N_A
        madelung_matrix(i,i) = madelung_matrix(i,i) - Hartree * 2.0_dp*epsilon/sqrt(PI)
    end do

    madelung_matrix = madelung_matrix - Hartree * PI/(vol*epsilon**2)
madelung_matrix = madelung_matrix + 0.61_dp

    deallocate(A)

end subroutine add_madelung_matrix

!| subroutine add_dmadelung_matrix
!| calculate derivative of Madelung matrix (interactions of pt charges in
!|    3-D periodic supercell)
!|
!|I integer N_A: number of atoms
!|I real(dp) a1(3), a2(3), a3(3): periodic supercell vectors
!|I real(dp) A(3,N_A): atomic positions
!|O real(dp) dmadelung_matrix(N_A,N_A): derivative of Madelung matrix
!|
!| Noam Bernstein 9/12/2001
!% This subroutine calculate derivative of Madelung matrix (interactions of pt charges in
!%  3-D periodic supercell).
subroutine add_dmadelung_matrix(N_A, a1_in, a2_in, a3_in, A_in, dmadelung_matrix)
    integer  :: N_A                             !% number of atoms
    real(dp) :: a1_in(3), a2_in(3), a3_in(3)   !% periodic supercell vectors
    real(dp) :: A_in(3,N_A)                    !% atomic positions
    real(dp) :: dmadelung_matrix(N_A,N_A,3)    !% derivative of Madelung matrix

    !NO PARALLEL real(dp), allocatable:: t_dmadelung_matrix(:,:,:)

    integer i, j
    integer s_a(3)
    real(dp) b1(3), b2(3), b3(3), h(3,3)
    real(dp) vol
    real(dp) K(3), ri(3), rj(3), dt, tv(3), Kmagsq, Rmag, dr(3)

    real(dp) :: a1(3), a2(3), a3(3)
    real(dp), allocatable :: A(:,:)

    integer i_pair
    integer rr, qq

    real(dp) epsilon, fourepsilonsq
    real(dp) recip_space_prefac, K_dep_factor
    real(dp) phase

    allocate(A(3,N_A))
    a1 = a1_in / Bohr
    a2 = a2_in / Bohr
    a3 = a3_in / Bohr
    A = A_in / Bohr

    h(1:3,1) = a1(1:3)
    h(1:3,2) = a2(1:3)
    h(1:3,3) = a3(1:3)
    vol = calc_volume(h)
    call cross_3 (a1, a2, b3)
    call cross_3 (a2, a3, b1)
    call cross_3 (a3, a1, b2)
    b1 = 2.0_dp*PI* b1 / vol
    b2 = 2.0_dp*PI* b2 / vol
    b3 = 2.0_dp*PI* b3 / vol

    !epsilon = 0.7_dp/nn_dist

    epsilon = 1.5_dp / ( 1.0_dp/sqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )
    !epsilon = 2.25_dp / ( 1.0_dp/sqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )

    fourepsilonsq = 4.0_dp*epsilon**2

    recip_space_prefac = 4.0_dp*PI/vol

    do qq=1, n_q

        s_a = q_list(1:3,qq)

        K = s_a(1)*b1 + s_a(2)*b2 + s_a(3)*b3
        Kmagsq = sum(K**2)
        K_dep_factor = recip_space_prefac*exp(-Kmagsq/fourepsilonsq)/Kmagsq
        do i_pair=1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            ri = A(1:3,i)
            rj = A(1:3,j)

            phase = sum(K*(ri-rj))
            dt = Hartree/Bohr * K_dep_factor*(-sin(phase))
            dmadelung_matrix(i,j,1:3) = dmadelung_matrix(i,j,1:3) + 2.0_dp*dt*K(1:3)
            if (i .ne. j) dmadelung_matrix(j,i,1:3) = dmadelung_matrix(j,i,1:3) - 2.0_dp*dt*K(1:3)
        end do

    end do

    tv = 0.0_dp
    do i_pair=1, n_pairs
        i = pair_list(1,i_pair)
        j = pair_list(2,i_pair)
        if (i .eq. j) cycle
        ri = A(1:3,i)
        rj = A(1:3,j)

        dr = ri-tv-rj
        Rmag = sqrt(sum(dr**2))
        dt = Hartree/Bohr * ( (-2.0_dp/sqrt(PI))*Rmag*epsilon*exp(-(epsilon*Rmag)**2) - erfc(epsilon*Rmag) ) / Rmag**2
        dmadelung_matrix(i,j,1:3) = dmadelung_matrix(i,j,1:3) + dt*dr/Rmag
        if (i .ne. j) dmadelung_matrix(j,i,1:3) = dmadelung_matrix(j,i,1:3) - dt*dr/Rmag

    end do

    do rr=1, n_r
        s_a = r_list(1:3,rr)

        tv = s_a(1)*a1 + s_a(2)*a2 + s_a(3)*a3
        do i_pair=1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            ri = A(1:3,i)
            rj = A(1:3,j)

            dr = ri-tv-rj
            Rmag = sqrt(sum(dr**2))
            dt = Hartree/Bohr * ( (-2.0_dp/sqrt(PI))*Rmag*epsilon*exp(-(epsilon*Rmag)**2) - erfc(epsilon*Rmag) ) / Rmag**2
            dmadelung_matrix(i,j,1:3) = dmadelung_matrix(i,j,1:3) + dt*dr/Rmag
            if (i .ne. j) dmadelung_matrix(j,i,1:3) = dmadelung_matrix(j,i,1:3) - dt*dr/Rmag

        end do
    end do

    !NO_PARALLEL allocate(t_dmadelung_matrix(N_A,N_A,3))
    !NO_PARALLEL t_dmadelung_matrix = dmadelung_matrix
    !NO_PARALLEL dmadelung_matrix = 0.0_dp
    !NO_PARALLEL call global_sum(t_dmadelung_matrix, dmadelung_matrix)
    !NO_PARALLEL deallocate(t_dmadelung_matrix)

    deallocate(A)

end subroutine add_dmadelung_matrix

!| subroutine add_dmadelung_matrix_dr
!% Calculate derivative of Madelung matrix (interactions of pt charges in
!%    3-D periodic supercell) dotted with atomic positions.
!% From Nielsen and Martin, Phys. Rev. B {\bf 32}, (1985).
!|
!|I integer N_A: number of atoms
!|I real(dp) a1(3), a2(3), a3(3): periodic supercell vectors
!|I real(dp) A(3,N_A): atomic positions
!|O real(dp) dmadelung_matrix_dr(N_A,N_A): derivative of Madelung matrix
!|    dotted with atomic positions
!|
!| Noam Bernstein 9/12/2001
!| From Nielsen and Martin, PRB _32_, 1985.

subroutine add_dmadelung_matrix_dr(N_A, a1_in, a2_in, a3_in, A_in, component, dmadelung_matrix_dr)
    integer  :: N_A                             !% number of atoms (input)
    real(dp) :: a1_in(3), a2_in(3), a3_in(3)    !% periodic supercell vectors (input)
    real(dp) :: A_in(3,N_A)                     !% atomic positions (input)
    integer  :: component
    real(dp) :: dmadelung_matrix_dr(N_A,N_A,3)  !% derivative of Madelung matrix dotted with atomic positions (output)

    integer i, j
    integer s_a(3)
    real(dp) b1(3), b2(3), b3(3), h(3,3)
    real(dp) vol
    real(dp) K(3), ri(3), rj(3), dt, tv(3), Kmagsq, Rmag, dr(3)
    real(dp) t_ij(3)

    real(dp) :: a1(3), a2(3), a3(3)
    real(dp), allocatable :: A(:,:)

    integer i_pair
    integer rr, qq

    real(dp) epsilon, fourepsilonsq
    real(dp) recip_space_prefac, K_dep_factor
    real(dp) phase

    allocate(A(3,N_A))
    a1 = a1_in / Bohr
    a2 = a2_in / Bohr
    a3 = a3_in / Bohr
    A = A_in / Bohr

    h(1:3,1) = a1(1:3)
    h(1:3,2) = a2(1:3)
    h(1:3,3) = a3(1:3)
    vol = calc_volume(h)
    call cross_3 (a1, a2, b3)
    call cross_3 (a2, a3, b1)
    call cross_3 (a3, a1, b2)
    b1 = 2.0_dp*PI* b1 / vol
    b2 = 2.0_dp*PI* b2 / vol
    b3 = 2.0_dp*PI* b3 / vol

    !epsilon = 0.7_dp/nn_dist

    epsilon = 1.5_dp / ( 1.0_dp/dsqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )
    !epsilon = 2.25_dp / ( 1.0_dp/dsqrt(PI) * 2.0_dp * ((3.0_dp*vol)/(4.0_dp*PI))**(1.0_dp/3.0_dp) )

    fourepsilonsq = 4.0_dp*epsilon**2

    recip_space_prefac = PI/(vol*epsilon**2)

    do qq=1, n_q

        s_a = q_list(1:3,qq)

        K = s_a(1)*b1 + s_a(2)*b2 + s_a(3)*b3
        Kmagsq = sum(K**2)

        K_dep_factor = recip_space_prefac*dexp(-Kmagsq/fourepsilonsq)/(Kmagsq/fourepsilonsq)

        do i_pair=1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            ri = A(1:3,i)
            rj = A(1:3,j)

            dr = ri-rj
            phase = sum(K*(ri-rj))

            dt = Hartree * K_dep_factor*2.0_dp*dcos(phase)

            if (component .eq. 0) then
                t_ij(1:3) = dt * (2.0_dp*K(1:3)*K(1:3)/Kmagsq * (Kmagsq/fourepsilonsq + 1.0_dp) - 1.0_dp)
                dmadelung_matrix_dr(i,j,1) = dmadelung_matrix_dr(i,j,1) - sum(t_ij(1:3))
                if (i .ne. j) dmadelung_matrix_dr(j,i,1) = dmadelung_matrix_dr(j,i,1) - sum(t_ij(1:3))
            else
                t_ij(1:3) = dt * (2.0_dp*K(component)*K(1:3)/Kmagsq * (Kmagsq/fourepsilonsq + 1.0_dp))
                t_ij(component) = t_ij(component) - dt
                dmadelung_matrix_dr(i,j,1:3) = dmadelung_matrix_dr(i,j,1:3) - t_ij(1:3)
                if (i .ne. j) dmadelung_matrix_dr(j,i,1:3) = dmadelung_matrix_dr(j,i,1:3) - t_ij(1:3)
            end if
        end do

    end do

    tv = 0.0_dp
    do i_pair=1, n_pairs
        i = pair_list(1,i_pair)
        j = pair_list(2,i_pair)
        if (i .eq. j) cycle
        ri = A(1:3,i)
        rj = A(1:3,j)

        dr = ri-tv-rj
        Rmag = sqrt(sum(dr**2))
        if (component .eq. 0) then
            t_ij(1:3) = Hartree * 0.5_dp*epsilon*2.0_dp*Hprime(epsilon*Rmag)*dr(1:3)*dr(1:3)/Rmag**2
            dmadelung_matrix_dr(i,j,1) = dmadelung_matrix_dr(i,j,1) - sum(t_ij(1:3))
            if (i .ne. j) dmadelung_matrix_dr(j,i,1) = dmadelung_matrix_dr(j,i,1) - sum(t_ij(1:3))
        else
            t_ij(1:3) = Hartree * 0.5_dp*epsilon*2.0_dp*Hprime(epsilon*Rmag)*dr(component)*dr(1:3)/Rmag**2
            dmadelung_matrix_dr(i,j,1:3) = dmadelung_matrix_dr(i,j,1:3) - t_ij(1:3)
            if (i .ne. j) dmadelung_matrix_dr(j,i,1:3) = dmadelung_matrix_dr(j,i,1:3) - t_ij(1:3)
        end if
    end do

    do rr=1, n_r
        s_a = r_list(1:3,rr)

        tv = s_a(1)*a1 + s_a(2)*a2 + s_a(3)*a3
        do i_pair=1, n_pairs
            i = pair_list(1,i_pair)
            j = pair_list(2,i_pair)
            ri = A(1:3,i)
            rj = A(1:3,j)

            dr = ri-tv-rj
            Rmag = sqrt(sum(dr**2))
            if (component .eq. 0) then
                t_ij(1:3) = Hartree * 0.5_dp*epsilon*2.0_dp*Hprime(epsilon*Rmag)*dr(1:3)*dr(1:3)/Rmag**2
                dmadelung_matrix_dr(i,j,1) = dmadelung_matrix_dr(i,j,1) - sum(t_ij(1:3))
                if (i .ne. j) dmadelung_matrix_dr(j,i,1) = dmadelung_matrix_dr(j,i,1) - sum(t_ij(1:3))
            else
                t_ij(1:3) = Hartree * 0.5_dp*epsilon*2.0_dp*Hprime(epsilon*Rmag)*dr(component)*dr(1:3)/Rmag**2
                dmadelung_matrix_dr(i,j,1:3) = dmadelung_matrix_dr(i,j,1:3) - t_ij(1:3)
                if (i .ne. j) dmadelung_matrix_dr(j,i,1:3) = dmadelung_matrix_dr(j,i,1:3) - t_ij(1:3)
            end if
        end do
    end do

    if (component .eq. 0) then
        dmadelung_matrix_dr(1:N_A,1:N_A,1) = dmadelung_matrix_dr(1:N_A,1:N_A,1) - &
            Hartree * 3.0_dp*PI/(vol*epsilon**2)
    else
        dmadelung_matrix_dr(1:N_A,1:N_A,1:3) = dmadelung_matrix_dr(1:N_A,1:N_A,1:3) - &
            Hartree * PI/(vol*epsilon**2)
    end if

!NO_PARALLEL    if (component .eq. 0) then
!NO_PARALLEL        allocate(t_dmadelung_matrix_dr(N_A,N_A,1))
!NO_PARALLEL        t_dmadelung_matrix_dr(1:N_A,1:N_A,1) = dmadelung_matrix_dr(1:N_A,1:N_A,1)
!NO_PARALLEL        dmadelung_matrix_dr(1:N_A,1:N_A,1) = 0.0_dp
!NO_PARALLEL        call global_sum(t_dmadelung_matrix_dr, dmadelung_matrix_dr)
!NO_PARALLEL        deallocate(t_dmadelung_matrix_dr)
!NO_PARALLEL    else
!NO_PARALLEL        allocate(t_dmadelung_matrix_dr(N_A,N_A,3))
!NO_PARALLEL        t_dmadelung_matrix_dr = dmadelung_matrix_dr
!NO_PARALLEL        dmadelung_matrix_dr = 0.0_dp
!NO_PARALLEL        call global_sum(t_dmadelung_matrix_dr, dmadelung_matrix_dr)
!NO_PARALLEL        deallocate(t_dmadelung_matrix_dr)
!NO_PARALLEL    endif

    deallocate(A)

end subroutine add_dmadelung_matrix_dr

!| subroutine cross_3
!% Calculate cross product of 2 3-vectors.  a(3), b(3): vectors to be crossed, c(3): cross product.
!|
!|I real(dp) a(3), b(3): vectors to be crossed
!|O real(dp) c(3): cross product
!|
!| Noam Bernstein 9/12/2001

subroutine cross_3(a, b, c)
    real(dp) a(3), b(3), c(3)

    c(1) = a(2)*b(3)-a(3)*b(2)
    c(2) = a(3)*b(1)-a(1)*b(3)
    c(3) = a(1)*b(2)-a(2)*b(1)

end subroutine cross_3

!| real(dp) function calc_volume
!| calculate volume of periodic supercell
!| returns volume of supercell
!|
!|I real(dp) a(3,3): matrix of peridic supercell vectors
!|
!| Noam Bernstein 9/12/2001

real(dp) function calc_volume(mat)
    real(dp) mat(3,3)

    calc_volume = abs(det_3_by_3(mat))
end function calc_volume

!| real(dp) function det_3_by_3
!| calculate determinant of 3x3 matrix
!| returns determinant of matrix
!|
!|I real(dp) a(3,3): matrix
!|
!| Noam Bernstein 9/12/2001

real(dp) function det_3_by_3(a)
    real(dp) a(3,3)

    det_3_by_3 = -a(1,3)*a(2,2)*a(3,1) + a(1,2)*a(2,3)*a(3,1) + a(1,3)*a(2,1)*a(3,2) - &
            a(1,1)*a(2,3)*a(3,2) - a(1,2)*a(2,1)*a(3,3) + a(1,1)*a(2,2)*a(3,3)
end function det_3_by_3

!| real(dp) function Hprime
!| calculate Hprime function for dMadelung_matrix/dr . r
!| returns value of Hprime(x)
!|
!|I real(dp) x: argument
!|
!| Noam Bernstein 9/12/2001

real(dp) function Hprime(x)
    real(dp) x

    Hprime =  (-2.0_dp/sqrt(PI))*exp(-x**2) - erfc(x)/x

end function Hprime

subroutine assign_atom_pairs(N, n_pairs, pair_list)
  integer, intent(in) :: N
  integer, intent(out) :: n_pairs
  integer, allocatable, intent(inout) :: pair_list(:,:)

  integer cur_pair_i, i, j

  n_pairs = N*(N+1)/2

  if (allocated(pair_list)) deallocate(pair_list)
  allocate(pair_list(2,n_pairs))

  cur_pair_i = 1
  do i=1, N
    do j=i, N
      pair_list(1,cur_pair_i) = i
      pair_list(2,cur_pair_i) = j
      cur_pair_i = cur_pair_i + 1
    end do
  end do

end subroutine assign_atom_pairs

end module Ewald_module
