hirshfeld.F90

Go to the documentation of this file.

!! Copyright (C) 2015 X. Andrade
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program 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 General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module hirshfeld_oct_m
  use atom_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_pot_oct_m
  use species_oct_m
  use splines_oct_m
 
  implicit none
 
  private
  public ::                       &
    hirshfeld_t,                  &
    hirshfeld_init,               &
    hirshfeld_end,                &
    hirshfeld_charge,             &
    hirshfeld_volume_ratio,       &
    hirshfeld_density_derivative, &
    hirshfeld_position_derivative
 
  type hirshfeld_t
    private
    class(mesh_t),       pointer     :: mesh
    type(atom_t),        pointer     :: atom(:)
    real(real64),        pointer     :: pos(:,:)
    integer                          :: natoms
    type(lattice_vectors_t), pointer :: latt
    real(real64),        allocatable :: total_density(:)
    real(real64),        allocatable :: free_volume(:)
    real(real64),        allocatable :: free_vol_r3(:,:)
    integer                          :: nspin
    type(namespace_t)                :: namespace
  end type hirshfeld_t
 
  real(real64), parameter, public :: TOL_HIRSHFELD = 1e-9_real64
 
contains
 
  subroutine hirshfeld_init(this, mesh, namespace, space, latt, atom, natoms, pos, nspin)
    type(hirshfeld_t),           intent(out)   :: this
    class(mesh_t),       target, intent(in)    :: mesh
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    type(lattice_vectors_t), target, intent(in):: latt
    type(atom_t), target,        intent(in)    :: atom(:)
    integer,                     intent(in)    :: natoms
    real(real64), target,        intent(in)    :: pos(1:space%dim,1:natoms)
    integer,                     intent(in)    :: nspin
 
    integer :: iatom, ip, isp
    real(real64) :: rr, atom_pos(space%dim), rmax, den
    real(real64), allocatable :: atom_density(:, :)
    type(ps_t), pointer :: ps
    type(lattice_iterator_t) :: latt_iter
    integer :: icell
    push_sub(hirshfeld_init)
 
    call profiling_in("HIRSHFELD_INIT")
 
    this%mesh => mesh
    this%atom => atom
    this%natoms = natoms
    this%nspin = nspin
    this%pos => pos
    this%latt => latt
    this%namespace = namespace
 
    safe_allocate(this%total_density(1:mesh%np))
    safe_allocate(this%free_volume(1:natoms))
    safe_allocate(this%free_vol_r3(1:mesh%np, 1:natoms))
    safe_allocate(atom_density(1:mesh%np, 1:nspin))
 
    !$omp parallel private(iatom)
    !$omp do simd
    do ip = 1, mesh%np
      this%total_density(ip) = m_zero
    end do
    !$omp end do simd nowait
    do iatom = 1, natoms
      !$omp do simd
      do ip = 1, mesh%np
        this%free_vol_r3(ip, iatom) = m_zero
      end do
      !$omp end do simd nowait
    end do
    !$omp end parallel
 
    do iatom = 1, natoms
      select type(spec=>atom(iatom)%species)
      type is(pseudopotential_t)
        ps => spec%ps
      class default
        assert(.false.)
      end select
 
      rmax = m_zero
      do isp = 1, nspin
        rmax = max(rmax, ps%density(isp)%x_threshold)
      end do
 
      latt_iter = lattice_iterator_t(latt, rmax)
      do icell = 1, latt_iter%n_cells
        atom_pos = pos(:, iatom) + latt_iter%get(icell)
        !We get the non periodized density
        !We need to do it to have the r^3 correctly computed for periodic systems
        call species_atom_density_np(atom(iatom)%species, namespace, atom_pos, mesh, nspin, atom_density)
 
        !$omp parallel do private(rr, den)
        do ip = 1, mesh%np
          rr = norm2(mesh%x(:, ip) - atom_pos)
          den = sum(atom_density(ip, 1:nspin))
 
          this%total_density(ip) = this%total_density(ip) + den
          this%free_vol_r3(ip, iatom) = this%free_vol_r3(ip, iatom) + den*rr**3
        end do
      end do
      this%free_volume(iatom) = dmf_integrate(this%mesh, this%free_vol_r3(:, iatom), reduce = .false.)
    end do
 
    call this%mesh%allreduce(this%free_volume)
 
    safe_deallocate_a(atom_density)
 
    call profiling_out("HIRSHFELD_INIT")
 
    pop_sub(hirshfeld_init)
  end subroutine hirshfeld_init
 
  ! ------------------------------------------------
 
  subroutine hirshfeld_end(this)
    type(hirshfeld_t), intent(inout) :: this
 
    push_sub(hirshfeld_end)
 
    safe_deallocate_a(this%total_density)
    safe_deallocate_a(this%free_volume)
    safe_deallocate_a(this%free_vol_r3)
 
    nullify(this%mesh)
    nullify(this%pos)
    nullify(this%atom)
    nullify(this%latt)
 
    pop_sub(hirshfeld_end)
  end subroutine hirshfeld_end
 
  ! -----------------------------------------------
 
  subroutine hirshfeld_charge(this, space, iatom, density, charge)
    type(hirshfeld_t),         intent(in)    :: this
    class(space_t),            intent(in)    :: space
    integer,                   intent(in)    :: iatom
    real(real64),              intent(in)    :: density(:, :)
    real(real64),              intent(out)   :: charge
 
    integer :: ip
    real(real64) :: dens_ip
    real(real64), allocatable :: atom_density(:, :), hirshfeld_density(:)
 
    push_sub(hirshfeld_charge)
 
    call profiling_in("HIRSHFELD_CHARGE")
 
    assert(allocated(this%total_density))
 
    safe_allocate(atom_density(1:this%mesh%np, 1:this%nspin))
    safe_allocate(hirshfeld_density(1:this%mesh%np))
 
    call species_atom_density(this%atom(iatom)%species, this%namespace, space, this%latt, &
      this%pos(:, iatom), this%mesh, this%nspin, atom_density)
 
    !$omp parallel do simd private(dens_ip)
    do ip = 1, this%mesh%np
      dens_ip = sum(atom_density(ip, 1:this%nspin))
      if (abs(dens_ip) > tol_hirshfeld) then
        hirshfeld_density(ip) = sum(density(ip, 1:this%nspin))*dens_ip/this%total_density(ip)
      else
        hirshfeld_density(ip) = 0.0_real64
      end if
    end do
 
    charge = dmf_integrate(this%mesh, hirshfeld_density)
 
    safe_deallocate_a(atom_density)
    safe_deallocate_a(hirshfeld_density)
 
    call profiling_out("HIRSHFELD_CHARGE")
 
    pop_sub(hirshfeld_charge)
  end subroutine hirshfeld_charge
 
  ! -----------------------------------------------
 
  subroutine hirshfeld_volume_ratio(this, iatom, density, volume_ratio)
    type(hirshfeld_t),         intent(in)    :: this
    integer,                   intent(in)    :: iatom
    real(real64),              intent(in)    :: density(:, :)
    real(real64),              intent(out)   :: volume_ratio
 
    integer :: ip
    real(real64), allocatable :: hirshfeld_density(:)
 
 
    push_sub(hirshfeld_volume_ratio)
 
    call profiling_in("HIRSHFELD_VOLUME_RATIO")
 
    assert(allocated(this%total_density))
 
    safe_allocate(hirshfeld_density(1:this%mesh%np))
 
    !$omp parallel do
    do ip = 1, this%mesh%np
      if (this%total_density(ip) > tol_hirshfeld) then
        hirshfeld_density(ip) = this%free_vol_r3(ip, iatom)*sum(density(ip, 1:this%nspin))/this%total_density(ip)
      else
        hirshfeld_density(ip) = 0.0_real64
      end if
    end do
 
    volume_ratio = dmf_integrate(this%mesh, hirshfeld_density)/this%free_volume(iatom)
 
    safe_deallocate_a(hirshfeld_density)
 
    call profiling_out("HIRSHFELD_VOLUME_RATIO")
 
    pop_sub(hirshfeld_volume_ratio)
  end subroutine hirshfeld_volume_ratio
 
  ! -----------------------------------------------
 
  subroutine hirshfeld_density_derivative(this, iatom, ddensity)
    type(hirshfeld_t),         intent(in)    :: this
    integer,                   intent(in)    :: iatom
    real(real64),              intent(out)   :: ddensity(:)
 
    integer :: ip
    real(real64), allocatable :: atom_density(:, :)
 
    push_sub(hirshfeld_density_derivative)
 
    call profiling_in("HIRSHFELD_DENSITY_DER")
 
    safe_allocate(atom_density(1:this%mesh%np, 1:this%nspin))
 
    !$omp parallel do simd
    do ip = 1, this%mesh%np
      if (abs(this%total_density(ip)) > tol_hirshfeld) then
        ddensity(ip) = this%free_vol_r3(ip, iatom)/(this%total_density(ip)*this%free_volume(iatom))
      else
        ddensity(ip) = m_zero
      end if
    end do
 
    safe_deallocate_a(atom_density)
 
    call profiling_out("HIRSHFELD_DENSITY_DER")
 
    pop_sub(hirshfeld_density_derivative)
  end subroutine hirshfeld_density_derivative
 
  ! -----------------------------------------------
  !dvadrr_ij = \frac{\delta V_i}{\delta \vec{x_j}}
  subroutine hirshfeld_position_derivative(this, space, iatom, jatom, density, dposition)
    type(hirshfeld_t),         intent(in)    :: this
    class(space_t),            intent(in)    :: space
    integer,                   intent(in)    :: iatom
    integer,                   intent(in)    :: jatom
    real(real64),              intent(in)    :: density(:, :)
    real(real64), contiguous,         intent(out)   :: dposition(:)
 
    integer :: ip, idir, icell, jcell, isp
    real(real64) :: atom_dens, atom_der,rri, rri2, rrj, tdensity, rij, rmax_isqu, rmax_jsqu
    real(real64) :: pos_i(space%dim), pos_j(space%dim), rmax_i, rmax_j
    real(real64), allocatable :: grad(:, :), atom_density(:, :), atom_derivative(:, :)
    type(lattice_iterator_t) :: latt_iter_i, latt_iter_j
    type(ps_t), pointer :: ps_i, ps_j
    real(real64) :: tmp, xxi(space%dim), xxj(space%dim)
 
    real(real64) :: tol_spacing
 
    push_sub(hirshfeld_position_derivative)
 
    tol_spacing = maxval(this%mesh%spacing(1:space%dim))
 
    call profiling_in("HIRSHFELD_POSITION_DER")
 
    dposition(1:space%dim) = m_zero
 
    select type(spec=> this%atom(iatom)%species)
    type is(pseudopotential_t)
      ps_i => spec%ps
    class default
      assert(.false.)
    end select
    select type(spec=> this%atom(jatom)%species)
    type is(pseudopotential_t)
      ps_j => spec%ps
    class default
      assert(.false.)
    end select
 
    rmax_i = m_zero
    rmax_j = m_zero
    do isp = 1, this%nspin
      rmax_i = max(rmax_i, ps_i%density(isp)%x_threshold)
      rmax_j = max(rmax_j, ps_j%density_der(isp)%x_threshold)
    end do
 
    rmax_isqu = rmax_i**2
    rmax_jsqu = rmax_j**2
 
    safe_allocate(grad(1:this%mesh%np, 1:space%dim))
    grad(1:this%mesh%np, 1:space%dim) = m_zero
    safe_allocate(atom_derivative(1:this%mesh%np, 1:this%nspin))
    safe_allocate(atom_density(1:this%mesh%np, 1:this%nspin))
 
    latt_iter_j = lattice_iterator_t(this%latt, rmax_j)
    do jcell = 1, latt_iter_j%n_cells
 
      pos_j = this%pos(:, jatom) + latt_iter_j%get(jcell)
      atom_derivative(1:this%mesh%np, 1:this%nspin) = m_zero
      call species_atom_density_derivative_np(this%atom(jatom)%species, this%namespace, pos_j, this%mesh, &
        min(this%nspin, 2), atom_derivative(1:this%mesh%np, 1:this%nspin))
 
      ! jcells further away from this distance cannot respect the following 'if' condition with respect to the i atom in this icell
      latt_iter_i = lattice_iterator_t(this%latt, (rmax_j+rmax_i))
      do icell = 1, latt_iter_i%n_cells
 
        pos_i = pos_j + latt_iter_i%get(icell) + (this%pos(:, iatom) - this%pos(:, jatom))
        rij =  norm2(pos_i - pos_j)
 
        if (rij - (rmax_j+rmax_i) < tol_spacing) then
 
          !We get the non periodized density
          !We need to do it to have the r^3 correctly computed for periodic systems
          call species_atom_density_np(this%atom(iatom)%species, this%namespace, pos_i, this%mesh, this%nspin, &
            atom_density(1:this%mesh%np, 1:this%nspin))
 
          !$omp parallel do private(xxi, rri, rri2, xxj, rrj, tdensity, atom_dens, atom_der, tmp, idir)
          do ip = 1, this%mesh%np
            if (this%total_density(ip)< tol_hirshfeld) cycle
 
            xxi = this%mesh%x(:, ip) - pos_i
            rri2 = sum(xxi**2)
            if (rri2 - rmax_isqu > tol_spacing) cycle ! In this case atom_dens = 0
 
            xxj = this%mesh%x(:, ip) - pos_j
            rrj = sum(xxj**2)
            if (rrj - rmax_jsqu > tol_spacing) cycle ! In this case atom_der = 0
 
            rri = sqrt(rri2)
            rrj = sqrt(rrj)
 
            tdensity = sum(density(ip, 1:this%nspin))
            atom_dens = sum(atom_density(ip, 1:this%nspin))
            atom_der = sum(atom_derivative(ip, 1:this%nspin))
 
            if (rrj > tol_hirshfeld) then
              tmp = rri*rri2*atom_dens*tdensity/this%total_density(ip)**2
              do idir = 1, space%dim
                grad(ip, idir) = grad(ip, idir) - tmp*atom_der*xxj(idir)/rrj
              end do
            end if
 
            !Only if we really have the same atoms
            if (iatom == jatom .and. rij < tol_hirshfeld) then
              grad(ip, :) = grad(ip, :) + (m_three*rri*atom_dens + rri2*atom_der)*tdensity/this%total_density(ip)*xxi
            end if
 
          end do
 
        end if
      end do
    end do
 
    dposition = dmf_integrate(this%mesh, space%dim, grad) / this%free_volume(iatom)
 
    safe_deallocate_a(atom_density)
    safe_deallocate_a(atom_derivative)
    safe_deallocate_a(grad)
 
    call profiling_out("HIRSHFELD_POSITION_DER")
 
    pop_sub(hirshfeld_position_derivative)
  end subroutine hirshfeld_position_derivative
 
end module hirshfeld_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: