! ! 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: V. Bolnykh !> compute interactions between quantum fragment and short-range classical fragments module mimic_short_range use mimic_cells use mimic_errors use mimic_field_grids use mimic_fragments use mimic_particles use mimic_precision, only: dp use mimic_utils, only: timer_start, timer_stop implicit none public :: compute_sr_energy_ions public :: compute_sr_forces_ions public :: compute_sr_forces_electrons public :: compute_sr_potential_electrons integer, parameter :: ORDER = 4; contains !> compute classical atom - quantum ion interaction energy subroutine compute_sr_energy_ions(atoms, nuclei, & energy, at_st, at_end) !> list of classical atoms type(atom_type), dimension(:), intent(inout) :: atoms !> list of ions to compute interactions with class(nucleus_type), dimension(:), intent(inout) :: nuclei !> resulting energy of interaction real(dp), intent(out) :: energy !> index of the atom to start electrostatic treatment integer, intent(in) :: at_st !> index of the atom to finish electrostatic treatment integer, intent(in) :: at_end integer :: i, j integer :: num_atoms real(dp), dimension(3) :: r_ij real(dp) :: rij real(dp) :: rij_sq real(dp) :: rij_4 real(dp) :: radius_4 real(dp) :: radius_5 real(dp) :: nominator real(dp) :: denominator real(dp) :: inv_denominator call timer_start("compute_sr_energy_ions") if (at_st > at_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "at_st should be smaller than at_end", & __FILE__, __LINE__) endif energy = 0.0_dp num_atoms = size(atoms) !$OMP PARALLEL DO private(i, radius_4, radius_5, j, r_ij, rij_sq, rij, & !$OMP& rij_4, nominator, denominator, inv_denominator) reduction(+:energy) do i = at_st, at_end if (abs(atoms(i)%charge) <= 1.0e-12_dp) cycle radius_4 = atoms(i)%radius**ORDER radius_5 = radius_4 * atoms(i)%radius do j = 1, size(nuclei) r_ij = atoms(i)%coordinate - nuclei(j)%coordinate rij_sq = r_ij(1)*r_ij(1) + r_ij(2)*r_ij(2) + r_ij(3)*r_ij(3) rij = sqrt(rij_sq) rij_4 = rij_sq * rij_sq nominator = radius_4 - rij_4 denominator = radius_5 - rij_4 * rij inv_denominator = 1.0_dp / denominator energy = energy + atoms(i)%charge * nuclei(j)%charge * & nominator * inv_denominator enddo enddo !$OMP END PARALLEL DO call timer_stop end subroutine compute_sr_energy_ions !> compute forces from classical atom - quantum ion interactions !> (stored in FION matrix (through pointers)) subroutine compute_sr_forces_ions(atoms, nuclei, & at_st, at_end) !> list of classical atoms type(atom_type), dimension(:), intent(inout) :: atoms !> list of ions to compute interactions with class(nucleus_type), dimension(:), intent(inout) :: nuclei !> index of the atom to start electrostatic treatment integer, intent(in) :: at_st !> index of the atom to finish electrostatic treatment integer, intent(in) :: at_end integer :: i, j integer :: num_atoms real(dp), dimension(3) :: r_ij real(dp), dimension(3) :: force_temp real(dp) :: rij real(dp) :: rij_sq real(dp) :: rij_4 real(dp) :: radius_4 real(dp) :: radius_5 real(dp) :: nominator real(dp) :: denominator real(dp) :: inv_denominator real(dp), dimension(:), allocatable :: temp_force call timer_start("compute_sr_forces_ions") if (at_st > at_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "at_st should be smaller than at_end", & __FILE__, __LINE__) endif num_atoms = size(atoms) allocate(temp_force(3 * size(nuclei))) temp_force = 0.0_dp !$OMP PARALLEL DO private(i, radius_4, radius_5, j, r_ij, rij_sq, rij, & !$OMP& rij_4, nominator, denominator, inv_denominator, force_temp) & !$OMP& reduction(-:temp_force) do i = at_st, at_end if (abs(atoms(i)%charge) <= 1.0e-12_dp) cycle radius_4 = atoms(i)%radius**ORDER radius_5 = radius_4 * atoms(i)%radius do j = 1, size(nuclei) r_ij = atoms(i)%coordinate - nuclei(j)%coordinate rij_sq = r_ij(1)*r_ij(1) + r_ij(2)*r_ij(2) + r_ij(3)*r_ij(3) rij = sqrt(rij_sq) rij_4 = rij_sq * rij_sq nominator = radius_4 - rij_4 denominator = radius_5 - rij_4 * rij inv_denominator = 1.0_dp / denominator force_temp = r_ij * atoms(i)%charge * nuclei(j)%charge * rij_sq * & (4.0_dp * denominator - 5.0_dp * rij * nominator) * & (inv_denominator * inv_denominator) atoms(i)%force = atoms(i)%force + force_temp temp_force((j-1)*3 + 1 : j*3) = temp_force((j-1)*3 + 1 : j*3) - force_temp enddo enddo !$OMP END PARALLEL DO !$OMP PARALLEL DO private(i) do i = 1, size(nuclei) nuclei(i)%force = nuclei(i)%force + temp_force((i-1)*3 + 1 : i*3) end do !$OMP END PARALLEL DO call timer_stop end subroutine compute_sr_forces_ions !> compute forces from classical atom - quantum electrons interactions !> (stored FION matrix (through pointers)) !> @todo lots of arguments, wrap into a type? @endtodo subroutine compute_sr_forces_electrons(atoms, density, cell, & x_start, x_end, id_start, id_end) USE, INTRINSIC :: ISO_C_BINDING !> list of atoms type(atom_type), dimension(:), intent(in) :: atoms !> electronic density real(dp), dimension(:,:,:), contiguous, intent(in), target :: density !> lattice data type(cell_type), intent(in) :: cell !> index of starting X-plane (CPMD parallelization stuff) integer, intent(in) :: x_start !> index of final X-plane (CPMD parallelization stuff) integer, intent(in) :: x_end !> index of the atom to start electrostatic treatment integer, intent(in) :: id_start !> index of the atom to finish electrostatic treatment integer, intent(in) :: id_end ! cell elementary vectors real(dp), dimension(3) :: da, db, dc ! coordinate of the j-th atom real(dp), dimension(3) :: r_j real(dp) :: force_x, force_y, force_z, temp_force real(dp) :: dV real(dp) :: rij, x_i, y_i, z_i, x_ij, y_ij, z_ij real(dp) :: rij_sq real(dp) :: charge real(dp) :: inv_denominator real(dp) :: point_potential real(dp) :: cov_r_4, cov_r_5 integer :: n_atom integer :: nz_step, n integer :: ny_step integer :: nx_step integer :: n_points integer :: nx, ny, nz integer, dimension(:), allocatable :: x_ind, y_ind real(dp), dimension(:), allocatable :: grid_x, grid_y, grid_z real(dp), dimension(:,:), contiguous, pointer :: density_p call timer_start("compute_sr_forces_electrons") if (id_start > id_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "id_start should be smaller than id_end", & __FILE__, __LINE__) endif if (x_start > x_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "x_start should be smaller than x_end", & __FILE__, __LINE__) endif n_points = cell%num_points(1) * cell%num_points(2) * cell%num_points(3) dV = cell%volume / real(n_points, dp) da = cell%lattice(:,1) / real(cell%num_points(1), dp) db = cell%lattice(:,2) / real(cell%num_points(2), dp) dc = cell%lattice(:,3) / real(cell%num_points(3), dp) nx = cell%num_points_r(1) ny = cell%num_points_r(2) nz = cell%num_points(3) call c_f_pointer(c_loc(density), density_p, [nx * ny, nz]) allocate(x_ind(nx * ny)) allocate(y_ind(nx * ny)) allocate(grid_x(nx * ny)) allocate(grid_y(nx * ny)) allocate(grid_z(nx * ny)) !$OMP PARALLEL DO private(nx_step) do ny_step = 1, ny do nx_step = 1, nx grid_x((ny_step - 1) * nx + nx_step) = db(1) * real(ny_step - 1, dp) + da(1) * (nx_step + x_start - 2) grid_y((ny_step - 1) * nx + nx_step) = db(2) * real(ny_step - 1, dp) + da(2) * (nx_step + x_start - 2) grid_z((ny_step - 1) * nx + nx_step) = db(3) * real(ny_step - 1, dp) + da(3) * (nx_step + x_start - 2) enddo enddo !$OMP END PARALLEL DO ! first calculate action of classical atoms on the electronic grid do n_atom = id_start, id_end if (abs(atoms(n_atom)%charge) <= 1.0e-12_dp) cycle force_x = 0.0_dp force_y = 0.0_dp force_z = 0.0_dp cov_r_4 = atoms(n_atom)%radius**ORDER cov_r_5 = cov_r_4 * atoms(n_atom)%radius r_j = atoms(n_atom)%coordinate charge = atoms(n_atom)%charge ! z-loop over grid !$OMP PARALLEL DO REDUCTION(-:force_x, force_y, force_z) SCHEDULE(DYNAMIC) DEFAULT(SHARED) & !$OMP& private(rij_sq, x_i, y_i, z_i, x_ij, y_ij, z_ij, & !$OMP& rij, inv_denominator, point_potential, temp_force, n) do nz_step = 1,nz x_i = cell%origin(1) + dc(1) * real(nz_step - 1, dp) y_i = cell%origin(2) + dc(2) * real(nz_step - 1, dp) z_i = cell%origin(3) + dc(3) * real(nz_step - 1, dp) ! y-loop over grid !dir$ vector aligned !$OMP SIMD REDUCTION(-:force_x, force_y, force_z) do n = 1, ny * nx ! x-loop over grid ! do nx_step = 1, x_end - x_start + 1 x_ij = r_j(1) - (x_i + grid_x(n)) y_ij = r_j(2) - (y_i + grid_y(n)) z_ij = r_j(3) - (z_i + grid_z(n)) rij_sq = x_ij * x_ij + y_ij * y_ij + z_ij * z_ij rij = sqrt(rij_sq) inv_denominator = 1.0_dp / (cov_r_5 - rij_sq * rij_sq * rij) ! external potential on a grid point point_potential = (cov_r_4 - rij_sq * rij_sq) * inv_denominator temp_force = density_p(n, nz_step) * & rij_sq * (4.0_dp - 5.0_dp * rij * point_potential) & * inv_denominator ! this is not a force yet force_x = force_x - x_ij * temp_force force_y = force_y - y_ij * temp_force force_z = force_z - z_ij * temp_force enddo !end n enddo !end nz_step !$OMP END PARALLEL DO atoms(n_atom)%force(1) = atoms(n_atom)%force(1) + force_x * atoms(n_atom)%charge * dV atoms(n_atom)%force(2) = atoms(n_atom)%force(2) + force_y * atoms(n_atom)%charge * dV atoms(n_atom)%force(3) = atoms(n_atom)%force(3) + force_z * atoms(n_atom)%charge * dV enddo !end n_atom call timer_stop end subroutine compute_sr_forces_electrons !> compute potential classical atom on quantum fragment !> @todo lots of arguments, wrap into a type? @endtodo subroutine compute_sr_potential_electrons(atoms, potential, cell, & x_start, x_end, id_start, id_end) USE, INTRINSIC :: ISO_C_BINDING !> list of SR atoms type(atom_type), dimension(:), intent(in) :: atoms !> potential on the electronic grid real(dp), dimension(:,:,:), contiguous, intent(inout), target :: potential !> lattice data type(cell_type), intent(in) :: cell !> index of starting X-plane (CPMD parallelization stuff) integer, intent(in) :: x_start !> index of final X-plane (CPMD parallelization stuff) integer, intent(in) :: x_end !> Index of the atom to start electrostatic treatment integer, intent(in) :: id_start !> Index of the atom to finish electrostatic treatment integer, intent(in) :: id_end ! cell elementary vectors real(dp), dimension(3) :: da, db, dc ! coordinate of the j-th atom real(dp), dimension(3) :: r_j real(dp) :: dV real(dp) :: rij, x_i, y_i, z_i, x_ij, y_ij, z_ij real(dp) :: rij_sq real(dp) :: charge real(dp) :: inv_denominator real(dp) :: point_potential real(dp) :: cov_r_4, cov_r_5 integer :: n_atom integer :: nz_step, n integer :: ny_step integer :: nx_step integer :: n_points integer :: nx, ny, nz, nx_real integer, dimension(:), allocatable :: x_ind, y_ind real(dp), dimension(:), allocatable :: grid_x, grid_y, grid_z real(dp), dimension(:,:), contiguous, pointer :: potential_p call timer_start("compute_sr_potential_electrons") if (id_start > id_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "id_start should be smaller than id_end", & __FILE__, __LINE__) endif if (x_start > x_end) then call handle_error(SEVERITY_FATAL, & TYPE_INCORRECT_ARG, & "x_start should be smaller than x_end", & __FILE__, __LINE__) endif n_points = cell%num_points(1) * cell%num_points(2) * cell%num_points(3) dV = cell%volume / real(n_points, dp) da = cell%lattice(:,1) / real(cell%num_points(1), dp) db = cell%lattice(:,2) / real(cell%num_points(2), dp) dc = cell%lattice(:,3) / real(cell%num_points(3), dp) nx = cell%num_points_r(1) ny = cell%num_points_r(2) nz = cell%num_points(3) call c_f_pointer(c_loc(potential), potential_p, [nx * ny, nz]) allocate(x_ind(nx * ny)) allocate(y_ind(nx * ny)) allocate(grid_x(nx * ny)) allocate(grid_y(nx * ny)) allocate(grid_z(nx * ny)) !$OMP PARALLEL DO private(nx_step) do ny_step = 1, ny do nx_step = 1, nx grid_x((ny_step - 1) * nx + nx_step) = db(1) * real(ny_step - 1, dp) + da(1) * (nx_step + x_start - 2) grid_y((ny_step - 1) * nx + nx_step) = db(2) * real(ny_step - 1, dp) + da(2) * (nx_step + x_start - 2) grid_z((ny_step - 1) * nx + nx_step) = db(3) * real(ny_step - 1, dp) + da(3) * (nx_step + x_start - 2) enddo enddo !$OMP END PARALLEL DO ! first calculate action of classical atoms on the electronic grid do n_atom = id_start, id_end if (abs(atoms(n_atom)%charge) <= 1.0e-12_dp) cycle cov_r_4 = atoms(n_atom)%radius**ORDER cov_r_5 = cov_r_4 * atoms(n_atom)%radius r_j = atoms(n_atom)%coordinate charge = atoms(n_atom)%charge ! z-loop over grid !$OMP PARALLEL DO SCHEDULE(DYNAMIC) DEFAULT(SHARED) & !$OMP& private(rij_sq, x_i, y_i, z_i, x_ij, y_ij, z_ij, & !$OMP& rij, inv_denominator, point_potential, n) do nz_step = 1,nz x_i = cell%origin(1) + dc(1) * real(nz_step - 1, dp) y_i = cell%origin(2) + dc(2) * real(nz_step - 1, dp) z_i = cell%origin(3) + dc(3) * real(nz_step - 1, dp) ! y-loop over grid !dir$ vector aligned !$OMP SIMD do n = 1, ny * nx ! x-loop over grid ! do nx_step = 1, x_end - x_start + 1 x_ij = r_j(1) - (x_i + grid_x(n)) y_ij = r_j(2) - (y_i + grid_y(n)) z_ij = r_j(3) - (z_i + grid_z(n)) rij_sq = x_ij * x_ij + y_ij * y_ij + z_ij * z_ij rij = sqrt(rij_sq) inv_denominator = 1.0_dp / (cov_r_5 - rij_sq * rij_sq * rij) ! external potential on a grid point point_potential = (cov_r_4 - rij_sq * rij_sq) * inv_denominator potential_p(n, nz_step) = & potential_p(n, nz_step) - point_potential * charge enddo !end ny_step enddo !end nz_step !$OMP END PARALLEL DO enddo !end n_atom nx_real = (x_end - x_start + 1) ! z-loop over grid !$OMP PARALLEL DO SCHEDULE(DYNAMIC) DEFAULT(SHARED) & !$OMP& private(nz_step, ny_step) do nz_step = 1,nz potential(:, ny, nz_step) = 0.0_dp if (nx /= nx_real) then do ny_step = 1, ny potential(nx, ny_step, nz_step) = 0.0_dp end do end if enddo !end nz_step !$OMP END PARALLEL DO call timer_stop end subroutine compute_sr_potential_electrons end module mimic_short_range