! ! 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/>. ! !> fragment types module mimic_fragments use mimic_cells use mimic_errors use mimic_field_grids use mimic_particles use mimic_precision, only: dp use mimic_properties use mimic_tensors use mimic_utils, only: timer_start, timer_stop implicit none private public :: fragment_type public :: quantum_fragment_type !> generic fragment type :: fragment_type private !> fragment id integer, public :: id !> number of atoms integer, public :: num_atoms !> total charge real(dp), public :: charge !> array of atoms type(atom_type), dimension(:), pointer, public :: atoms contains private procedure, public :: init => init_fragment procedure, public :: centroid => fragment_centroid procedure, public :: center_of_mass => fragment_center_of_mass procedure, public :: center_of_charge => fragment_center_of_charge end type fragment_type !> quantum system representation !> (currently can only be one per simulation) type :: quantum_fragment_type private !> id of the quantum fragment integer, public :: id !> number of nuclei integer, public :: num_nuclei !> total charge real(dp), public :: charge !> nuclear charge real(dp), public :: nuclear_charge !> electron charge real(dp), public :: electronic_charge !> array of nuclei type(nucleus_type), dimension(:), allocatable, public :: nuclei !> simulation cell data type(cell_type), public :: cell !> grid representation of the electron density type(density_type), public :: density !> grid representation of the external potential type(potential_type), public :: potential !> nuclear multipole moments type(multipoles_type), public :: nuclear_multipoles !> electronic multipole moments type(multipoles_type), public :: electronic_multipoles contains private procedure, public :: init => init_quantum_fragment procedure, public :: compute_nuclear_multipoles procedure, public :: compute_electronic_multipoles procedure, public :: centroid => quantum_fragment_centroid procedure, public :: center_of_mass => quantum_fragment_center_of_mass procedure, public :: center_of_charge => quantum_fragment_center_of_charge end type quantum_fragment_type contains !> intialize fragment !> @todo calculate fragment charge (perhaps test for closeness to integer) subroutine init_fragment(this, id, atoms) class(fragment_type), intent(inout) :: this integer, intent(in) :: id type(atom_type), dimension(:), target, intent(in) :: atoms integer :: i this%id = id this%num_atoms = size(atoms) this%atoms => atoms this%charge = 0.0_dp do i = 1, this%num_atoms this%charge = this%charge + this%atoms(i)%charge end do end subroutine init_fragment !> compute the fragment centroid pure function fragment_centroid(this) result(centroid) class(fragment_type), intent(in) :: this real(dp), dimension(3) :: centroid centroid = compute_centroid(this%atoms) end function fragment_centroid !> compute the fragment center-of-mass pure function fragment_center_of_mass(this) result(center_of_mass) class(fragment_type), intent(in) :: this real(dp), dimension(3) :: center_of_mass center_of_mass = compute_center_of_mass(this%atoms) end function fragment_center_of_mass !> compute the fragment center-of-charge pure function fragment_center_of_charge(this) result(center_of_charge) class(fragment_type), intent(in) :: this real(dp), dimension(3) :: center_of_charge center_of_charge = compute_center_of_charge(this%atoms) end function fragment_center_of_charge !> Initialize quantum fragment subroutine init_quantum_fragment(this, id, nuclei, cell, density, potential) class(quantum_fragment_type), intent(inout) :: this integer, intent(in) :: id type(nucleus_type), dimension(:), intent(in) :: nuclei type(cell_type), intent(in) :: cell type(density_type), intent(in) :: density type(potential_type), optional, intent(in) :: potential this%id = id this%num_nuclei = size(nuclei) allocate(this%nuclei(this%num_nuclei)) this%nuclei = nuclei this%cell = cell this%density = density if (present(potential)) then this%potential = potential end if end subroutine init_quantum_fragment !> compute quantum-fragment centroid pure function quantum_fragment_centroid(this) result(centroid) class(quantum_fragment_type), intent(in) :: this real(dp), dimension(3) :: centroid centroid = compute_centroid(this%nuclei) end function quantum_fragment_centroid !> compute quantum-fragment center-of-mass pure function quantum_fragment_center_of_mass(this) result(center_of_mass) class(quantum_fragment_type), intent(in) :: this real(dp), dimension(3) :: center_of_mass center_of_mass = compute_center_of_mass(this%nuclei) end function quantum_fragment_center_of_mass !> compute quantum-fragment center-of-charge pure function quantum_fragment_center_of_charge(this) result(center_of_charge) class(quantum_fragment_type), intent(in) :: this real(dp), dimension(3) :: center_of_charge center_of_charge = compute_center_of_charge(this%nuclei) end function quantum_fragment_center_of_charge !> compute multipole expansion of the electrostatic potential from !> a set of nuclei around an origin subroutine compute_nuclear_multipoles(this, origin, order) class(quantum_fragment_type), intent(inout) :: this !> origin of the multipole expansion real(dp), dimension(3), intent(in) :: origin !> order of the multipole expansion integer, intent(in) :: order integer :: nuc, ord, comp integer :: i, j, k integer :: x, y, z real(dp) :: temp real(dp), dimension(3) :: r real(dp), dimension(:), allocatable :: multipoles call timer_start("compute_nuclear_multipoles") allocate(multipoles((order+3)*(order+2)*(order+1)/6)) multipoles = 0.0_dp do nuc = 1, this%num_nuclei r = this%nuclei(nuc)%coordinate - origin comp = 1 do ord = 0, order do i = ord, 0, -1 do j = ord - i, 0, -1 k = ord - i - j temp = this%nuclei(nuc)%charge 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 multipoles(comp) = multipoles(comp) + temp comp = comp + 1 end do end do end do end do call this%nuclear_multipoles%init(id=0, origin=origin, order=order, values=multipoles) deallocate(multipoles) call timer_stop end subroutine compute_nuclear_multipoles !> compute multipole expansion of the electrostatic potential from !> an electronic density (which is represented on a grid) around an origin subroutine compute_electronic_multipoles(this, expansion_origin, expansion_order, & a_start, a_end) class(quantum_fragment_type), intent(inout) :: this !> origin of the multipole expansion real(dp), dimension(3), intent(in) :: expansion_origin !> order of the multipole expansion integer, intent(in) :: expansion_order !> 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 integer :: i, j, k integer :: x, y, z integer :: a, b, c integer :: a_index integer :: mi integer :: order real(dp) :: dV real(dp) :: temp real(dp), dimension(3) :: ra, rb, rc, r real(dp), dimension(3) :: cell_origin real(dp), dimension(3,3) :: lattice_steps real(dp), dimension(:), allocatable :: multipoles call timer_start("compute_electronic_multipoles") 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 allocate(multipoles(polytensor_size(expansion_order))) multipoles = 0.0_dp cell_origin = this%cell%origin do i = 1, 3 lattice_steps(:,i) = this%cell%lattice(:,i) / real(this%cell%num_points(i), dp) end do dV = this%cell%volume / real(product(this%cell%num_points), dp) associate (density => this%density%field) !$OMP PARALLEL DO PRIVATE(order, i, j, k, c, rc, b, rb, a, a_index, ra, r, temp, & !$OMP& x, y, z, mi) REDUCTION(+:multipoles) do c = 1, this%cell%num_points(3) rc = cell_origin + real(c - 1, dp) * lattice_steps(:,3) do b = 1, this%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 - expansion_origin mi = 1 do order = 0, expansion_order do i = order, 0, -1 do j = order - i, 0, -1 k = order - i - j temp = - density(a_index,b,c) * dV 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 multipoles(mi) = multipoles(mi) + temp mi = mi + 1 end do end do end do end do end do end do !$OMP END PARALLEL DO end associate call this%electronic_multipoles%init(id=0, origin=expansion_origin, order=expansion_order, values=multipoles) deallocate(multipoles) call timer_stop end subroutine compute_electronic_multipoles end module mimic_fragments