mimic_long_range.F90 Source File

!
!    MiMiC: A Framework for Multiscale Modeling in Computational Chemistry
!    Copyright (C) 2015-2022  The MiMiC Contributors (see CONTRIBUTORS file for details).
!
!    This file is part of MiMiC.
!
!    MiMiC is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as
!    published by the Free Software Foundation, either version 3 of
!    the License, or (at your option) any later version.
!
!    MiMiC is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!

!> author: Jógvan Magnus Haugaard Olsen
!> computes energy, forces, and potential from the interactions between the
!> quantum fragment and long-range atoms using a multipole expansion
!> approximation for the quantum fragment
module mimic_long_range

    use mimic_cells
    use mimic_errors
    use mimic_fragments
    use mimic_particles
    use mimic_precision, only: dp
    use mimic_tensors
    use mimic_utils, only: timer_start, timer_stop

    implicit none

    private

    public :: compute_lr_energy
    public :: compute_lr_potential
    public :: compute_lr_forces
    public :: compute_folded_tensors

contains

!> Computes folded tensor used in energy and potential computations
subroutine compute_folded_tensors(tensors, tensor_sums, &
                                  multipole_origin, multipole_order, &
                                  atoms)

    !> Array storing tensors
    real(dp), dimension(:), intent(out) :: tensors
    !> Array stroing tensor sums
    real(dp), dimension(:), intent(out) :: tensor_sums
    !> Coordinates of the origin of multipolar expansion
    real(dp), dimension(3), intent(in) :: multipole_origin
    !> Maximal order of multipoles
    integer, intent(in) :: multipole_order
    !> long-range atoms
    type(atom_type), dimension(:), intent(in) :: atoms

    integer :: i, j, k, order
    real(dp) :: taylor
    real(dp), dimension(3) :: r

    call timer_start("compute_folded_tensors")

    tensors = 0.0_dp
    tensor_sums = 0.0_dp

    !$OMP PARALLEL DO private(i, j, r, k, order, taylor, tensors) &
    !$OMP& reduction(+:tensor_sums)
    do i = 1, size(atoms)
        r = multipole_origin - atoms(i)%coordinate
        k = 0
        do order = 0, multipole_order
            j = k + 1
            k = k + tensor_size(order)
            call folded_interaction_tensor(r, tensors(j:k), order, 0)
            taylor = merge(1.0_dp, - 1.0_dp, (mod(order, 2) == 0)) / factorial(order)
            tensors(j:k) = taylor * tensors(j:k)
        end do
        tensor_sums = tensor_sums + atoms(i)%charge * tensors
    end do
    !$OMP END PARALLEL DO

    call timer_stop

end subroutine compute_folded_tensors

!> compute interaction energy between quantum fragment and long-range atoms
subroutine compute_lr_energy(quantum_fragment, &
                             energy, tensor_sums)

    !> quantum fragment
    type(quantum_fragment_type), intent(in) :: quantum_fragment
    !> interaction energy
    real(dp), intent(out) :: energy
    !> Array stroing tensor sums
    real(dp), dimension(:), intent(in) :: tensor_sums

    call timer_start("compute_lr_energy")

    energy = 0.0_dp

    energy = dot_product(quantum_fragment%nuclear_multipoles%values, tensor_sums)

    call timer_stop

end subroutine compute_lr_energy

!> compute potential from long-range atoms on quantum fragment
subroutine compute_lr_potential(quantum_fragment, &
                                a_start, a_end, tensor_sums)

    !> quantum fragment
    type(quantum_fragment_type), intent(in) :: quantum_fragment
    !> start index of real space mesh of lattice vector a
    integer, intent(in) :: a_start
    !> end index of real space mesh of lattice vector a
    integer, intent(in) :: a_end
    !> Array stroing tensor sums
    real(dp), dimension(:), intent(in) :: tensor_sums

    integer :: i, j, k
    integer :: x, y, z
    integer :: a, b, c
    integer :: a_index
    integer :: comp
    integer :: order
    integer :: multipole_order
    real(dp) :: temp
    real(dp), dimension(3) :: ra, rb, rc, r
    real(dp), dimension(3) :: multipole_origin
    real(dp), dimension(3,3) :: lattice_steps

    call timer_start("compute_lr_potential")

    if (a_start > a_end) then
        call handle_error(SEVERITY_FATAL, &
                          TYPE_INCORRECT_ARG, &
                          "a_start should be smaller than a_end", &
                          __FILE__, __LINE__)
    endif

    multipole_origin = quantum_fragment%nuclear_multipoles%origin
    multipole_order = quantum_fragment%nuclear_multipoles%order

    associate (cell => quantum_fragment%cell)

    do i = 1, 3
        lattice_steps(:,i) = cell%lattice(:,i) / real(cell%num_points(i), dp)
    end do

    associate (potential => quantum_fragment%potential%field)

    !$OMP PARALLEL DO private(c, rc, b, rb, a, a_index, ra, r, &
    !$OMP& comp, order, i, j, k, temp)
    do c = 1, cell%num_points(3)
        rc = cell%origin + real(c - 1, dp) * lattice_steps(:,3)
        do b = 1, cell%num_points(2)
            rb = rc + real(b - 1, dp) * lattice_steps(:,2)
            do a = a_start, a_end
                a_index = a - a_start + 1
                ra = rb + real(a - 1, dp) * lattice_steps(:,1)
                r = ra - multipole_origin
                comp = 1
                do order = 0, multipole_order
                    do i = order, 0, - 1
                        do j = order - i, 0, - 1
                            k = order - i - j
                            temp = 1.0_dp
                            do x = 1, i
                                temp = temp * r(1)
                            end do
                            do y = 1, j
                                temp = temp * r(2)
                            end do
                            do z = 1, k
                                temp = temp * r(3)
                            end do
                            potential(a_index,b,c) = potential(a_index,b,c) + temp * tensor_sums(comp)
                            comp = comp + 1
                        end do
                    end do
                end do
            end do
        end do
    end do
    !$OMP END PARALLEL DO

    end associate

    end associate

    call timer_stop

end subroutine compute_lr_potential

!> compute forces from long-range atoms on quantum fragment and vice versa
subroutine compute_lr_forces(atoms, quantum_fragment, &
                             atom_start, atom_end)

    !> list of atoms in the long-range part
    type(atom_type), dimension(:), intent(in) :: atoms
    !> the quantum fragment
    type(quantum_fragment_type), intent(in) :: quantum_fragment
    !> index of the first atom that current process will handle
    integer, intent(in) :: atom_start
    !> index of the last atom that the current process will handle
    integer, intent(in) :: atom_end

    integer :: i, j
    integer :: x, y, z
    integer :: tx, ty, tz
    integer :: mx, my, mz
    integer :: rx, ry, rz
    integer :: fx, fy, fz
    integer :: mi, fi, ti
    integer :: t_start, t_end, m_start, m_end
    integer :: rank_f, rank_t, rank_m
    integer :: order, expansion_order
    real(dp) :: coeff, temp
    real(dp), dimension(3) :: r
    real(dp), dimension(3) :: force
    real(dp), dimension(3) :: expansion_origin
    real(dp), dimension(:), allocatable :: multipoles
    real(dp), dimension(:), allocatable :: tensor, tensor_sum
    real(dp), dimension(:), allocatable :: polytensor, polytensor_sum

    call timer_start("compute_lr_forces")

    if (atom_start > atom_end) then
        call handle_error(SEVERITY_FATAL, &
                          TYPE_INCORRECT_ARG, &
                          "atom_start should be smaller than atom_end", &
                          __FILE__, __LINE__)
    endif

    expansion_origin = quantum_fragment%nuclear_multipoles%origin
    expansion_order = quantum_fragment%nuclear_multipoles%order

    allocate(multipoles(polytensor_size(expansion_order)))
    multipoles = quantum_fragment%nuclear_multipoles%values + &
                 quantum_fragment%electronic_multipoles%values

    allocate(tensor(tensor_size(expansion_order + 1)))
    allocate(tensor_sum(tensor_size(expansion_order + 1)))

    allocate(polytensor(polytensor_size(expansion_order)))
    allocate(polytensor_sum(polytensor_size(expansion_order)))

    !$OMP PARALLEL DO private(i, r, force, order, m_start, m_end, &
    !$OMP& rank_t, rank_m, rank_f, fx, fy, fz, fi, tx, ty, tz, ti, &
    !$OMP& mx, my, mz, mi, coeff, tensor)
    do i = atom_start, atom_end
        r = expansion_origin - atoms(i)%coordinate
        force = 0.0_dp
        m_end = 0
        do order = 0, expansion_order
            m_start = m_end + 1
            m_end = m_start + tensor_size(order) - 1
            call interaction_tensor(r, tensor(1:tensor_size(order+1)), order, 1)
            rank_t = tensor_rank(tensor(1:tensor_size(order+1)))
            rank_m = tensor_rank(multipoles(m_start:m_end))
            rank_f = tensor_rank(force)
            do fx = rank_f, 0, -1
                do fy = rank_f - fx, 0, -1
                    fz = rank_f - fx - fy
                    fi = convert_index(fx, fy, fz)
                    do tx = rank_t, fx, -1
                        do ty = rank_t - tx, fy, -1
                            tz = rank_t - tx - ty
                            if (tz < fz) cycle
                            ti = convert_index(tx, ty, tz)
                            mx = tx - fx
                            my = ty - fy
                            mz = tz - fz
                            mi = convert_index(mx, my, mz)
                            coeff = real((-1)**(tx+ty+tz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
                            force(fi) = force(fi) + coeff * tensor(ti) * multipoles(m_start-1+mi)
                        end do
                    end do
                end do
            end do
        end do
        atoms(i)%force = atoms(i)%force - atoms(i)%charge * force
    end do
    !$OMP END PARALLEL DO

    tensor_sum = 0.0_dp
    !$OMP PARALLEL DO private(i, r, tensor) reduction(+:tensor_sum)
    do i = atom_start, atom_end
        r = expansion_origin - atoms(i)%coordinate
        call interaction_tensor(r, tensor, expansion_order, 1)
        tensor_sum = tensor_sum + atoms(i)%charge * tensor
    end do
    !$OMP END PARALLEL DO

    force = 0.0_dp
    m_start = polytensor_size(expansion_order - 1) + 1
    m_end = m_start + tensor_size(expansion_order) - 1
    rank_t = tensor_rank(tensor_sum)
    rank_m = tensor_rank(multipoles(m_start:m_end))
    rank_f = tensor_rank(force)
    do fx = rank_f, 0, -1
        do fy = rank_f - fx, 0, -1
            fz = rank_f - fx - fy
            fi = convert_index(fx, fy, fz)
            do tx = rank_t, fx, -1
                do ty = rank_t - tx, fy, -1
                    tz = rank_t - tx - ty
                    if (tz < fz) cycle
                    ti = convert_index(tx, ty, tz)
                    mx = tx - fx
                    my = ty - fy
                    mz = tz - fz
                    mi = convert_index(mx, my, mz)
                    coeff = real((-1)**(mx+my+mz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
                    force(fi) = force(fi) + coeff * tensor_sum(ti) * multipoles(m_start-1+mi)
                end do
            end do
        end do
    end do

    force = force / real(quantum_fragment%num_nuclei, dp)
    do j = 1, quantum_fragment%num_nuclei
        quantum_fragment%nuclei(j)%force = quantum_fragment%nuclei(j)%force - force
    end do

    polytensor = 0.0_dp
    polytensor_sum = 0.0_dp
    !$OMP PARALLEL DO private(i, r, order, t_start, t_end, polytensor) &
    !$OMP& reduction(+:polytensor_sum)
    do i = atom_start, atom_end
        r = expansion_origin - atoms(i)%coordinate
        t_end = 1
        do order = 1, expansion_order
            t_start = t_end + 1
            t_end = t_start + tensor_size(order) - 1
            call interaction_tensor(r, polytensor(t_start:t_end), order, 0)
        end do
        polytensor_sum = polytensor_sum + atoms(i)%charge * polytensor
    end do
    !$OMP END PARALLEL DO

    !$OMP PARALLEL DO private(j, r, multipoles, order, rx, ry, rz, temp, &
    !$OMP& force, t_start, t_end, m_start, m_end, rank_t, rank_m, rank_f, &
    !$OMP& fx, fy, fz, fi, tx, ty, tz, ti, mx, my, mz, mi, coeff)
    do j = 1, quantum_fragment%num_nuclei
        r = quantum_fragment%nuclei(j)%coordinate - expansion_origin
        multipoles = 0.0_dp
        mi = 1
        do order = 0, expansion_order - 1
            do rx = order, 0, -1
                do ry = order - rx, 0, -1
                    rz = order - rx - ry
                    temp = quantum_fragment%nuclei(j)%charge
                    do x = 1, rx
                        temp = temp * r(1)
                    end do
                    do y = 1, ry
                        temp = temp * r(2)
                    end do
                    do z = 1, rz
                        temp = temp * r(3)
                    end do
                    multipoles(mi) = multipoles(mi) + temp
                    mi = mi + 1
                end do
            end do
        end do
        force = 0.0_dp
        t_end = 1
        m_end = 0
        do order = 1, expansion_order
            t_start = t_end + 1
            t_end = t_start + tensor_size(order) - 1
            m_start = m_end + 1
            m_end = m_start + tensor_size(order-1) - 1
            rank_t = tensor_rank(polytensor_sum(t_start:t_end))
            rank_m = tensor_rank(multipoles(m_start:m_end))
            rank_f = tensor_rank(force)
            do fx = rank_f, 0, -1
                do fy = rank_f - fx, 0, -1
                    fz = rank_f - fx - fy
                    fi = convert_index(fx, fy, fz)
                    do tx = rank_t, fx, -1
                        do ty = rank_t - tx, fy, -1
                            tz = rank_t - tx - ty
                            if (tz < fz) cycle
                            ti = convert_index(tx, ty, tz)
                            mx = tx - fx
                            my = ty - fy
                            mz = tz - fz
                            mi = convert_index(mx, my, mz)
                            coeff = real((-1)**(tx+ty+tz), dp) / (factorial(mx) * factorial(my) * factorial(mz))
                            force(fi) = force(fi) + coeff * polytensor_sum(t_start-1+ti) * multipoles(m_start-1+mi)
                        end do
                    end do
                end do
            end do
        end do
        quantum_fragment%nuclei(j)%force = quantum_fragment%nuclei(j)%force - force
    end do
    !$OMP END PARALLEL DO

    deallocate(multipoles)
    deallocate(tensor)
    deallocate(tensor_sum)
    deallocate(polytensor)
    deallocate(polytensor_sum)

    call timer_stop

end subroutine compute_lr_forces

end module mimic_long_range