ion_electron_local_potential.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 Nicolas Tancogne-Dejean
!!
!! 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 ion_electron_local_potential_oct_m
  use atom_oct_m
  use comm_oct_m
  use debug_oct_m
  use distributed_oct_m
  use epot_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use poisson_oct_m
  use potential_interaction_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use quantity_oct_m
  use space_oct_m
  use species_oct_m
  use species_pot_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use submesh_oct_m
 
  implicit none
 
  private
  public ::                &
    ion_electron_local_potential_t
 
  type, extends(potential_interaction_t) :: ion_electron_local_potential_t
    private
 
    type(poisson_t), pointer :: psolver
    type(mesh_t),    pointer :: mesh
    type(space_t),   pointer :: space
 
    ! This information should be copied, and not obtained thanks to pointer
    ! This is a temporary change here
    type(distributed_t), pointer, public :: atoms_dist => null()
    type(atom_t), pointer, public :: atom(:) => null()
    real(real64), pointer, public :: pos(:,:) => null()
    type(lattice_vectors_t), pointer :: latt => null()
 
    ! Temporary pointer to namespace
    type(namespace_t), pointer :: namespace
 
    logical :: have_density
 
  contains
    procedure :: init => ion_electron_local_potential_init
    procedure :: bind => ion_electron_local_potential_bind
    procedure :: calculate => ion_electron_local_potential_calculate
    procedure :: calculate_energy => ion_electron_local_potential_calculate_energy
    procedure :: end => ion_electron_local_potential_end
    final :: ion_electron_local_potential_finalize
  end type ion_electron_local_potential_t
 
  interface ion_electron_local_potential_t
    module procedure ion_electron_local_potential_constructor
  end interface ion_electron_local_potential_t
 
contains
 
  ! ---------------------------------------------------------
  function ion_electron_local_potential_constructor(partner) result(this)
    class(interaction_partner_t), target, intent(inout) :: partner
    class(ion_electron_local_potential_t),               pointer       :: this
 
    push_sub(ion_electron_local_potential_constructor)
 
    allocate(this)
 
    this%label = "ion-electron local"
 
    this%partner => partner
 
    pop_sub(ion_electron_local_potential_constructor)
  end function ion_electron_local_potential_constructor
 
  ! ---------------------------------------------------------
  subroutine ion_electron_local_potential_init(this, mesh, psolver, ions, namespace)
    class(ion_electron_local_potential_t),    intent(inout) :: this
    class(mesh_t),                    target, intent(in)    :: mesh
    type(poisson_t),                  target, intent(in)    :: psolver
    type(ions_t),                     target, intent(in)    :: ions
    type(namespace_t),                target, intent(in)    :: namespace
 
    integer :: ia
 
    push_sub(ion_electron_local_potential_init)
 
    this%mesh => mesh
    this%psolver => psolver
 
    safe_allocate(this%potential(1:mesh%np, 1))
 
    this%have_density = .false.
    do ia = 1, ions%nspecies
      if (local_potential_has_density(ions%space, ions%species(ia)%s)) then
        this%have_density = .true.
        exit
      end if
    end do
 
    this%atoms_dist => ions%atoms_dist
    this%atom => ions%atom
    this%space => ions%space
    this%pos => ions%pos
    this%latt => ions%latt
 
    this%namespace => namespace
 
    pop_sub(ion_electron_local_potential_init)
  end subroutine ion_electron_local_potential_init
 
  ! ---------------------------------------------------------
  subroutine ion_electron_local_potential_bind(this, psolver, ions)
    class(ion_electron_local_potential_t), intent(inout) :: this
    type(poisson_t),               target, intent(in)    :: psolver
    type(ions_t),                  target, intent(in)    :: ions
 
    push_sub(ion_electron_local_potential_bind)
 
    this%psolver => psolver
    this%atoms_dist => ions%atoms_dist
    this%atom => ions%atom
    this%space => ions%space
    this%pos => ions%pos
    this%latt => ions%latt
 
    pop_sub(ion_electron_local_potential_bind)
  end subroutine ion_electron_local_potential_bind
 
  ! ---------------------------------------------------------
  ! Note: this code is adapted from epot_local_potential. There are code duplication at the moment
  subroutine ion_electron_local_potential_calculate(this)
    class(ion_electron_local_potential_t),             intent(inout) :: this
 
    real(real64), allocatable :: density(:), rho(:), vl(:)
    type(submesh_t) :: sphere
    type(ps_t), pointer :: ps
    integer :: ia, ip
    real(real64) :: radius
    logical :: sphere_init
 
    push_sub(ion_electron_local_potential_calculate)
 
    call profiling_in("ION_ELEC_LOC_INTER")
 
    this%potential(:,1) = m_zero
 
    if (this%have_density) then
      safe_allocate(density(1:this%mesh%np))
      density = m_zero
    end if
 
    sphere_init = .false.
 
    do ia = this%atoms_dist%start, this%atoms_dist%end
      ! Local potential, we can get it by solving the Poisson equation
      ! (for all-electron species or pseudopotentials in periodic
      ! systems) or by applying it directly to the grid
 
      if (local_potential_has_density(this%space, this%atom(ia)%species)) then
 
        safe_allocate(rho(1:this%mesh%np))
        call species_get_long_range_density(this%atom(ia)%species, this%namespace, this%space, this%latt, &
          this%pos(:, ia), this%mesh, rho, sphere)
        call lalg_axpy(this%mesh%np, m_one, rho, density)
        safe_deallocate_a(rho)
        sphere_init = .true.
 
      else
 
        safe_allocate(vl(1:this%mesh%np))
        call species_get_local(this%atom(ia)%species, this%namespace, this%space, this%latt, &
          this%pos(:, ia), this%mesh, vl)
        call lalg_axpy(this%mesh%np, m_one, vl, this%potential(:,1))
        safe_deallocate_a(vl)
 
      end if
 
      !then we add the localized part
      select type(spec=>this%atom(ia)%species)
      type is(pseudopotential_t)
 
        ps => spec%ps
 
        radius = min(ps%vl%x_threshold*1.05_real64, spline_range_max(ps%vl))
        if (sphere_init) then
          if ( .not. submesh_compatible(sphere,radius,this%pos(:,ia), minval(this%mesh%spacing(1:this%space%dim))) ) then
            call submesh_init(sphere, this%space, this%mesh, this%latt, this%pos(:, ia), radius)
          end if
        else
          call submesh_init(sphere, this%space, this%mesh, this%latt, this%pos(:, ia), radius)
        end if
        safe_allocate(vl(1:sphere%np))
 
        do ip = 1, sphere%np
          if(sphere%r(ip) <= radius) then
            vl(ip) = spline_eval(ps%vl, sphere%r(ip))
          else
            vl(ip) = 0
          end if
        end do
 
        call submesh_add_to_mesh(sphere, vl, this%potential(:,1))
 
        safe_deallocate_a(vl)
        nullify(ps)
      end select
      call submesh_end(sphere)
    end do
 
    ! reduce over atoms if required
    if (this%atoms_dist%parallel) then
      call comm_allreduce(this%atoms_dist%mpi_grp, this%potential(:,1), dim = this%mesh%np)
      if (this%have_density) then
        call comm_allreduce(this%atoms_dist%mpi_grp, density, dim = this%mesh%np)
      end if
    end if
 
    ! now we solve the poisson equation with the density of all nodes
    if (this%have_density) then
 
      safe_allocate(vl(1:this%mesh%np_part))
      call dpoisson_solve(this%psolver, this%namespace, vl, density)
      call lalg_axpy(this%mesh%np, m_one, vl, this%potential(:,1))
      safe_deallocate_a(vl)
 
    end if
    safe_deallocate_a(density)
 
    call profiling_out("ION_ELEC_LOC_INTER")
 
    pop_sub(ion_electron_local_potential_calculate)
  end subroutine ion_electron_local_potential_calculate
 
  ! ---------------------------------------------------------
  subroutine ion_electron_local_potential_end(this)
    class(ion_electron_local_potential_t), intent(inout) :: this
 
    push_sub(ion_electron_local_potential_end)
 
    safe_deallocate_a(this%potential)
    nullify(this%mesh)
    nullify(this%psolver)
 
    call interaction_end(this)
 
    pop_sub(ion_electron_local_potential_end)
  end subroutine ion_electron_local_potential_end
 
 
  ! ---------------------------------------------------------
  subroutine ion_electron_local_potential_finalize(this)
    type(ion_electron_local_potential_t), intent(inout) :: this
 
    push_sub(ion_electron_local_potential_finalize)
 
    call ion_electron_local_potential_end(this)
 
    pop_sub(ion_electron_local_potential_finalize)
  end subroutine ion_electron_local_potential_finalize
 
  ! ---------------------------------------------------------
  subroutine ion_electron_local_potential_calculate_energy(this)
    class(ion_electron_local_potential_t),             intent(inout) :: this
 
    push_sub(ion_electron_local_potential_calculate_energy)
 
    assert(.false.)
 
    pop_sub(ion_electron_local_potential_calculate_energy)
 
  end  subroutine ion_electron_local_potential_calculate_energy
 
 
end module ion_electron_local_potential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: