species_pot.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 species_pot_oct_m
  use allelectron_oct_m
  use atomic_orbital_oct_m
  use debug_oct_m
  use global_oct_m
  use io_function_oct_m
  use index_oct_m
  use, intrinsic :: iso_fortran_env
  use jellium_oct_m
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use logrid_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use root_solver_oct_m
  use space_oct_m
  use species_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use submesh_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use volume_oct_m
 
  implicit none
 
  private
  public ::                         &
    species_get_long_range_density, &
    species_get_nlcc,               &
    species_get_nlcc_grad,          &
    species_get_local,              &
    species_atom_density,           &
    species_atom_density_np,        &
    species_atom_density_derivative,&
    species_atom_density_derivative_np, &
    species_atom_density_grad
 
  type(mesh_t), pointer :: mesh_p
  real(real64), allocatable :: rho_p(:)
  real(real64), allocatable :: grho_p(:, :)
  real(real64) :: alpha2_p
  real(real64), pointer :: pos_p(:)
 
contains
 
 
  ! ---------------------------------------------------------
  subroutine species_atom_density(species, namespace, space, latt, pos, mesh, spin_channels, rho)
    class(species_t), target, intent(in)   :: species
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    type(mesh_t),            intent(in)    :: mesh
    integer,                 intent(in)    :: spin_channels
    real(real64),            intent(inout) :: rho(:, :)
 
    integer :: isp, ip, in_points, icell
    real(real64) :: rr, x, pos_pc(space%dim), nrm, rmax
    real(real64) :: xx(space%dim), yy(space%dim), rerho, imrho
    real(real64), allocatable :: dorbital(:)
    type(ps_t), pointer :: ps
    type(volume_t) :: volume
    integer :: in_points_red
    type(lattice_iterator_t) :: latt_iter
    integer :: iorb, ii, nn, ll, mm
    real(real64) :: radius, density
    type(submesh_t) :: sphere
 
    push_sub(species_atom_density)
 
    assert(spin_channels == 1 .or. spin_channels == 2)
 
    rho = m_zero
 
    ! build density ...
    select type (species)
    type is(species_from_file_t)
      call generate_uniform_density()
 
    type is(species_user_defined_t)
      call generate_uniform_density()
 
    type is(soft_coulomb_t)
      call generate_uniform_density()
 
    class is(allelectron_t)
 
      do isp = 1, spin_channels
        do iorb = 1, species%get_niwfs()
          call species%get_iwf_ilm(iorb, isp, ii, ll, mm)
          ! For all-electron species, we want to use the principal quantum number
          call species%get_iwf_n(iorb, isp, nn)
 
          radius = species%get_iwf_radius(nn, isp)
          ! make sure that if the spacing is too large, the orbitals fit in a few points at least
          radius = max(radius, m_two*maxval(mesh%spacing))
 
          call submesh_init(sphere, space, mesh, latt, pos, radius)
          safe_allocate(dorbital(1:sphere%np))
 
          call datomic_orbital_get_submesh(species, sphere, nn, ll, mm, isp, dorbital)
          ! The occupations are for one type of orbitals, e.g. 2p gets 6 electrons
          ! So we normalize them by (2*l+1) such that they get distributed evenly
          ! for each value of m
          do ip = 1, sphere%np
            dorbital(ip) = species%conf%occ(ii, isp)/real(2*ll+1, real64) *dorbital(ip)*dorbital(ip)
          end do
          call submesh_add_to_mesh(sphere, dorbital, rho(:, isp))
          safe_deallocate_a(dorbital)
 
          call submesh_end(sphere)
        end do
      end do
 
    type is (jellium_charge_t)
      ! We put, for the electron density, the same as the positive density that
      ! creates the external potential.
      ! This code is repeated in get_density, and should therefore be cleaned!!!!!
 
      call volume_init(volume)
      call volume_read_from_block(volume, namespace, trim(species%rho_string()))
 
      rmax = latt%max_length()
      latt_iter = lattice_iterator_t(latt, rmax)
      rho = m_zero
      do icell = 1, latt_iter%n_cells
        yy = latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = pos, coords = xx)
          xx = xx + yy
          rr = norm2(xx)
 
          rerho = m_zero
          if (volume_in_volume(space, volume, xx)) rerho = m_one
          rho(ip, 1) = rho(ip, 1) + rerho
        end do
      end do
 
      call volume_end(volume)
 
      if (spin_channels > 1) then
        rho(:, 1) = m_half*rho(:, 1)
        rho(:, 2) = rho(:, 1)
      end if
 
      ! rescale to match the valence charge
      do isp = 1, spin_channels
        x = species%get_zval() / dmf_integrate(mesh, rho(:, isp))
        rho(1:mesh%np, isp) = x * rho(1:mesh%np, isp)
      end do
 
    type is (species_charge_density_t)
      ! We put, for the electron density, the same as the positive density that
      ! creates the external potential.
      ! This code is repeated in get_density, and should therefore be cleaned!!!!!
 
      rmax = latt%max_length()
      latt_iter = lattice_iterator_t(latt, rmax)
      rho = m_zero
      do icell = 1, latt_iter%n_cells
        yy = latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = pos, coords = xx)
          xx = xx + yy
          rr = norm2(xx)
 
          rerho = m_zero
          call parse_expression(rerho, imrho, space%dim, xx, rr, m_zero, trim(species%rho_string()))
          rho(ip, 1) = rho(ip, 1) + rerho
        end do
      end do
 
      if (spin_channels > 1) then
        rho(:, 1) = m_half*rho(:, 1)
        rho(:, 2) = rho(:, 1)
      end if
 
      ! rescale to match the valence charge
      do isp = 1, spin_channels
        x = species%get_zval() / dmf_integrate(mesh, rho(:, isp))
        rho(1:mesh%np, isp) = x * rho(1:mesh%np, isp)
      end do
 
 
    type is (jellium_sphere_t) ! ... from jellium
      in_points = 0
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, rr, origin = pos)
        if (rr <= species%radius()) then
          in_points = in_points + 1
        end if
      end do
 
      if (mesh%parallel_in_domains) then
        call mesh%mpi_grp%allreduce(in_points, in_points_red, 1, mpi_integer, mpi_sum)
        in_points = in_points_red
      end if
 
      if (in_points > 0) then
        ! This probably should be done inside the mesh_function_oct_m module.
 
        if (mesh%use_curvilinear) then
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos)
            if (rr <= species%radius()) then
              rho(ip, 1:spin_channels) = species%get_zval() /   &
                (mesh%vol_pp(ip) * real(in_points*spin_channels, real64) )
            end if
          end do
        else
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos)
            if (rr <= species%radius()) then
              rho(ip, 1:spin_channels) = species%get_zval() /   &
                (mesh%vol_pp(1) * real(in_points * spin_channels, real64) )
            end if
          end do
        end if
      end if
 
    type is (jellium_slab_t) ! ... from jellium slab
      density = species%get_density(mesh%box%bounding_box_l) / spin_channels
 
      do ip = 1, mesh%np
        rr = abs(mesh%x(ip, 3) - pos(3))
        if (rr <= species%thickness() / m_two) then
          rho(ip, 1:spin_channels) = density
        end if
      end do
 
    class is (pseudopotential_t)
      ! ...from pseudopotentials
 
      ps => species%ps
 
      if (ps_has_density(ps)) then
 
        assert(allocated(ps%density))
 
        rmax = m_zero
        do isp = 1, spin_channels
          rmax = max(rmax, ps%density(isp)%x_threshold)
        end do
 
        latt_iter = lattice_iterator_t(latt, rmax)
        do icell = 1, latt_iter%n_cells
          pos_pc = pos + latt_iter%get(icell)
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos_pc)
            rr = max(rr, r_small)
 
            do isp = 1, spin_channels
              if (rr >= spline_range_max(ps%density(isp))) cycle
              rho(ip, isp) = rho(ip, isp) + spline_eval(ps%density(isp), rr)
            end do
 
          end do
        end do
 
      else
 
        !we use the square root of the short-range local potential, just to put something that looks like a density
 
        latt_iter = lattice_iterator_t(latt, ps%vl%x_threshold)
        do icell = 1, latt_iter%n_cells
          pos_pc = pos + latt_iter%get(icell)
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos_pc)
            rr = max(rr, r_small)
 
            if (rr >= spline_range_max(ps%vl)) cycle
 
            do isp = 1, spin_channels
              rho(ip, isp) = rho(ip, isp) + sqrt(abs(spline_eval(ps%vl, rr)))
            end do
 
          end do
        end do
 
        ! normalize
        nrm = m_zero
        do isp = 1, spin_channels
          nrm = nrm + dmf_integrate(mesh, rho(:, isp))
        end do
 
        rho(1:mesh%np, 1:spin_channels) = rho(1:mesh%np, 1:spin_channels)*species%get_zval()/nrm
 
      end if
    class default
      assert(.false.)
    end select
 
    pop_sub(species_atom_density)
  contains
    subroutine generate_uniform_density()
      do isp = 1, spin_channels
        rho(1:mesh%np, isp) = m_one
        x = (species%get_zval()/real(spin_channels, real64) ) / dmf_integrate(mesh, rho(:, isp))
        rho(1:mesh%np, isp) = x * rho(1:mesh%np, isp)
      end do
    end subroutine generate_uniform_density
  end subroutine species_atom_density
 
  ! ---------------------------------------------------------
  ! A non periodized version of the routine species_atom_density
  ! This is used for the Hirshfeld routines
  ! TODO: implement it for other approaches than pseudo potentials.
  subroutine species_atom_density_np(species, namespace, pos, mesh, spin_channels, rho)
    class(species_t), target, intent(in)    :: species
    type(namespace_t),        intent(in)    :: namespace
    real(real64),             intent(in)    :: pos(:)
    type(mesh_t),             intent(in)    :: mesh
    integer,                  intent(in)    :: spin_channels
    real(real64),             intent(inout) :: rho(:, :)
 
    integer :: isp, ip
    real(real64) :: rr, nrm
    type(ps_t), pointer :: ps
 
    push_sub(species_atom_density_np)
 
    call profiling_in("SPECIES_ATOM_DEN_NP")
 
    rho = m_zero
    select type(species)
    class is(pseudopotential_t)
      ! ...from pseudopotentials
 
      ps => species%ps
      if (ps_has_density(ps)) then
 
        assert(allocated(ps%density))
 
        !$omp parallel private(ip, rr)
        do isp = 1, spin_channels
          !$omp do
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos)
            if (rr >= spline_range_max(ps%density(isp))) cycle
            rr = max(rr, r_small)
            rho(ip, isp) = rho(ip, isp) + spline_eval(ps%density(isp), rr)
          end do
          !$omp end do nowait
        end do
        !$omp end parallel
 
      else
 
        !we use the square root of the short-range local potential, just to put something that looks like a density
 
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = pos)
          rr = max(rr, r_small)
 
          if (rr >= spline_range_max(ps%vl)) cycle
 
          do isp = 1, spin_channels
            rho(ip, isp) = rho(ip, isp) + sqrt(abs(spline_eval(ps%vl, rr)))
          end do
 
        end do
 
        ! normalize
        nrm = m_zero
        do isp = 1, spin_channels
          nrm = nrm + dmf_integrate(mesh, rho(:, isp))
        end do
 
        rho(1:mesh%np, 1:spin_channels) = rho(1:mesh%np, 1:spin_channels)*species%get_zval()/nrm
 
      end if
    class default
      call messages_not_implemented('species_atom_density_np for non-pseudopotential species', namespace=namespace)
 
    end select
 
    call profiling_out("SPECIES_ATOM_DEN_NP")
 
    pop_sub(species_atom_density_np)
  end subroutine species_atom_density_np
 
 
  ! ---------------------------------------------------------
 
  subroutine species_atom_density_derivative(species, namespace, space, latt, pos, mesh, spin_channels, drho)
    class(species_t), target,intent(in)    :: species
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    type(mesh_t),            intent(in)    :: mesh
    integer,                 intent(in)    :: spin_channels
    real(real64),            intent(inout) :: drho(:, :)
 
    integer :: icell
    real(real64) :: pos_pc(space%dim), range
    type(ps_t), pointer :: ps
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(species_atom_density_derivative)
 
    assert(spin_channels == 1 .or. spin_channels == 2)
 
    drho = m_zero
 
    ! build density ...
    select type(species)
    class is(pseudopotential_t)
      ! ...from pseudopotentials
      ps => species%ps
 
      if (ps_has_density(ps)) then
 
        range = ps%density_der(1)%x_threshold
        if (spin_channels == 2) range = max(range, ps%density_der(2)%x_threshold)
        latt_iter = lattice_iterator_t(latt, range)
 
        do icell = 1, latt_iter%n_cells
          pos_pc = pos + latt_iter%get(icell)
          call species_atom_density_derivative_np(species, namespace, pos_pc, mesh, spin_channels, drho)
        end do
      end if
 
    class default
      call messages_not_implemented('species_atom_density_derivative for non-pseudopotential species', namespace=namespace)
 
    end select
 
    pop_sub(species_atom_density_derivative)
  end subroutine species_atom_density_derivative
 
  ! ---------------------------------------------------------
  !! Non-periodic version of the above routine
  subroutine species_atom_density_derivative_np(species, namespace, pos, mesh, spin_channels, drho)
    class(species_t), target, intent(in)    :: species
    type(namespace_t),        intent(in)    :: namespace
    real(real64),             intent(in)    :: pos(:)
    type(mesh_t),             intent(in)    :: mesh
    integer,                  intent(in)    :: spin_channels
    real(real64),             intent(inout) :: drho(:, :)
 
    integer :: isp, ip
    real(real64) :: rr
    type(ps_t), pointer :: ps
 
    push_sub(species_atom_density_derivative_np)
 
    call profiling_in("SPECIES_ATOM_DEN_DER_NP")
 
    select type(species)
    class is(pseudopotential_t)
      ps => species%ps
 
      if (ps_has_density(ps)) then
        !$omp parallel private(ip, rr, isp)
        do isp = 1, spin_channels
          !$omp do
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos)
            if (rr >= spline_range_max(ps%density_der(isp))) cycle
            rr = max(rr, r_small)
            drho(ip, isp) = drho(ip, isp) + spline_eval(ps%density_der(isp), rr)
          end do
          !$omp end do nowait
        end do
        !$omp end parallel
 
      else
        call messages_write('The pseudopotential for')
        call messages_write(species%get_label())
        call messages_write(' does not contain the density.')
        call messages_fatal(namespace=namespace)
      end if
    class default
      assert(.false.)
    end select
 
    call profiling_out("SPECIES_ATOM_DEN_DER_NP")
 
    pop_sub(species_atom_density_derivative_np)
  end subroutine species_atom_density_derivative_np
 
 
  ! ---------------------------------------------------------
  !! Gradient of the atomic density, if available
  subroutine species_atom_density_grad(species, namespace, space, latt, pos, mesh, spin_channels, drho)
    class(species_t), target, intent(in)   :: species
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    type(mesh_t),            intent(in)    :: mesh
    integer,                 intent(in)    :: spin_channels
    real(real64),            intent(inout) :: drho(:, :, :)
 
    integer :: isp, ip, icell, idir
    real(real64) :: rr, pos_pc(space%dim), range, spline
    type(ps_t), pointer :: ps
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(species_atom_density_grad)
 
    assert(spin_channels == 1 .or. spin_channels == 2)
 
    drho = m_zero
 
    ! build density ...
    select type(species)
    class is(pseudopotential_t)
      ps => species%ps
      ! ...from pseudopotentials
 
      if (ps_has_density(ps)) then
 
        range = ps%density_der(1)%x_threshold
        if (spin_channels == 2) range = max(range, ps%density_der(2)%x_threshold)
        latt_iter = lattice_iterator_t(latt, range)
 
        do icell = 1, latt_iter%n_cells
          pos_pc = pos + latt_iter%get(icell)
 
          do ip = 1, mesh%np
            call mesh_r(mesh, ip, rr, origin = pos_pc)
            rr = max(rr, r_small)
 
            do isp = 1, spin_channels
              if (rr >= spline_range_max(ps%density_der(isp))) cycle
              spline = spline_eval(ps%density_der(isp), rr)
 
              if(abs(spline) < 1e-150_real64) cycle
 
              do idir = 1, space%dim
                drho(ip, isp, idir) = drho(ip, isp, idir) - spline*(mesh%x(ip, idir) - pos_pc(idir))/rr
              end do
            end do
          end do
        end do
 
      else
        call messages_write('The pseudopotential for')
        call messages_write(species%get_label())
        call messages_write(' does not contain the density.')
        call messages_fatal(namespace=namespace)
      end if
 
    class default
      call messages_not_implemented('species_atom_density_grad for non-pseudopotential species', namespace=namespace)
 
    end select
 
    pop_sub(species_atom_density_grad)
  end subroutine species_atom_density_grad
 
  ! ---------------------------------------------------------
 
  subroutine species_get_long_range_density(species, namespace, space, latt, pos, mesh, rho, sphere_inout, nlr_x)
    class(species_t), target,   intent(in)    :: species
    type(namespace_t),          intent(in)    :: namespace
    class(space_t),             intent(in)    :: space
    type(lattice_vectors_t),    intent(in)    :: latt
    real(real64),       target, intent(in)    :: pos(1:space%dim)
    class(mesh_t),      target, intent(in)    :: mesh
    real(real64),               intent(out)   :: rho(:)
    type(submesh_t), optional, target, intent(inout) :: sphere_inout
    real(real64),    optional,  intent(inout) :: nlr_x(:,:)
 
    type(root_solver_t) :: rs
    logical :: conv
    real(real64)   :: startval(space%dim)
    real(real64)   :: delta, alpha, xx(space%dim), yy(space%dim), rr, imrho1, rerho
    real(real64)   :: dist2_min
    integer :: icell, ipos, ip, idir, rankmin
    type(lattice_iterator_t) :: latt_iter
    type(ps_t), pointer :: ps
    type(volume_t) :: volume
    type(submesh_t), target  :: sphere_local
    type(submesh_t), pointer :: sphere
    logical :: have_point
    real(real64),    allocatable :: rho_sphere(:)
    real(real64), parameter      :: threshold = 1e-6_real64
    real(real64)          :: norm_factor, range, radius, radius_nlr, radius_vl
 
    push_sub(species_get_long_range_density)
 
    call profiling_in("SPECIES_LR_DENSITY")
 
    if(present(nlr_x)) then
      assert(species%is_ps())
    end if
 
    select type (species)
    type is(pseudopotential_t)
      ps => species%ps
      radius_nlr = spline_x_threshold(ps%nlr, threshold)
      if (present(sphere_inout)) then
        radius_vl  = ps%vl%x_threshold*1.05_real64
        radius = max(radius_nlr, radius_vl)
        call submesh_init(sphere_inout, space, mesh, latt, pos, radius)
        sphere => sphere_inout
      else
        radius = radius_nlr
        call submesh_init(sphere_local, space, mesh, latt, pos, radius)
        sphere => sphere_local
      endif
 
      safe_allocate(rho_sphere(1:sphere%np))
      if (.not. present(sphere_inout) .and. sphere%np > 0) then
        do ip = 1, sphere%np
          rho_sphere(ip) = sphere%r(ip)
        end do
        call spline_eval_vec(ps%nlr, sphere%np, rho_sphere)
      else
        do ip = 1, sphere%np
          if(sphere%r(ip) <= radius_nlr) then
            rho_sphere(ip) = spline_eval(ps%nlr, sphere%r(ip))
          else
            rho_sphere(ip) = m_zero
          endif
        end do
      end if
 
      rho(1:mesh%np) = m_zero
 
      ! A small amount of charge is missing with the cutoff, we
      ! renormalize so that the long range potential is exact
      norm_factor = abs(species%get_zval()/dsm_integrate(mesh, sphere, rho_sphere))
      do ip = 1, sphere%np
        rho(sphere%map(ip)) = rho(sphere%map(ip)) + norm_factor*rho_sphere(ip)
      end do
 
      if (present(nlr_x)) then
        do idir = 1, space%dim
          do ip = 1, sphere%np
            nlr_x(sphere%map(ip), idir) = nlr_x(sphere%map(ip), idir) + norm_factor*rho_sphere(ip)*sphere%rel_x(idir, ip)
          end do
        end do
      end if
 
      safe_deallocate_a(rho_sphere)
      nullify(ps)
      if ( .not. present(sphere_inout) ) then
        call submesh_end(sphere)
      end if
      nullify(sphere)
 
    type is (full_delta_t)
 
      rho(1:mesh%np) = m_zero
 
      ipos = mesh_nearest_point(mesh, pos, dist2_min, rankmin)
      have_point = .true.
      if (mesh%mpi_grp%rank /= rankmin) have_point = .false.
 
      if (have_point) then
        if (mesh%use_curvilinear) then
          rho(ipos) = -species%get_z()/mesh%vol_pp(ipos)
        else
          rho(ipos) = -species%get_z()/mesh%vol_pp(1)
        end if
      end if
 
      write(message(1), '(3a,f5.2,3a)') &
        "Info: species_full_delta species ", trim(species%get_label()), &
        " atom displaced ", units_from_atomic(units_out%length, sqrt(dist2_min)), &
        " [ ", trim(units_abbrev(units_out%length)), " ]"
      call messages_info(1, namespace=namespace)
 
    type is (full_gaussian_t)
 
      ! periodic copies are not considered in this routine
      if (space%is_periodic()) then
        call messages_not_implemented("species_full_gaussian for periodic systems", namespace=namespace)
      end if
 
      ! We need to work with \xi and \xi_0, not x(\xi) and x(\xi_0) as we do now
      ! for the argument of the Gaussian
      if (mesh%use_curvilinear) then
        call messages_not_implemented("species_full_gaussian with curvilinear coordinates", namespace=namespace)
      end if
 
      ! --------------------------------------------------------------
      ! Constructs density for an all-electron atom with the procedure
      ! sketched in Modine et al. [Phys. Rev. B 55, 10289 (1997)],
      ! section II.B
      ! --------------------------------------------------------------
 
      safe_allocate(rho_p(1:mesh%np))
      safe_allocate(grho_p(1:mesh%np, 1:space%dim))
 
      mesh_p => mesh
      pos_p => pos
 
      ! Initial guess.
      delta   = mesh%spacing(1)
      alpha   = sqrt(m_two)*species%get_sigma()*delta
      alpha2_p = alpha**2  ! global copy of alpha
 
      ! the dim variables are the position of the delta function
      startval(1:space%dim) = pos
 
      ! solve equation
      ! Setting a tolerance such that the distance to the first moment is smaller than 1e-5 Bohr
      call root_solver_init(rs, namespace, space%dim, solver_type=root_newton, maxiter=500, abs_tolerance=1.0e-10_real64)
      call droot_solver_run(rs, func, xx, conv, startval=startval)
 
      if (.not. conv) then
        write(message(1),'(a)') 'Root finding in species_get_density did not converge.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if(debug%info .and. space%dim == 3) then
        write(message(1),'(a,3(f6.3,a))') 'Debug: Gaussian charge position (', xx(1), ', ', xx(2), ', ', xx(3), ')'
        call messages_info(1, namespace=namespace)
      end if
 
      ! we want a charge of -Z
      rho = -species%get_z()*rho_p
 
      nullify(mesh_p)
      nullify(pos_p)
      safe_deallocate_a(grho_p)
      safe_deallocate_a(rho_p)
 
    type is (full_anc_t)
 
      rho = m_zero
 
    type is(jellium_charge_t)
 
      call volume_init(volume)
      call volume_read_from_block(volume, namespace, trim(species%rho_string()))
 
      range = latt%max_length()
      latt_iter = lattice_iterator_t(latt, range)
 
      rho = m_zero
      do icell = 1, latt_iter%n_cells
        yy = latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = pos, coords = xx)
          xx = xx + yy
          rr = norm2(xx)
 
          rerho = m_zero
          if (volume_in_volume(space, volume, xx)) rerho = m_one
          rho(ip) = rho(ip) - rerho
        end do
      end do
 
      call volume_end(volume)
 
    type is(species_charge_density_t)
 
      range = latt%max_length()
      latt_iter = lattice_iterator_t(latt, range)
 
      rho = m_zero
      do icell = 1, latt_iter%n_cells
        yy = latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = pos, coords = xx)
          xx = xx + yy
          rr = norm2(xx)
 
          rerho = m_zero
          call parse_expression(rerho, imrho1, space%dim, xx, rr, m_zero, trim(species%rho_string()))
          rho(ip) = rho(ip) - rerho
        end do
      end do
 
      rr = species%get_zval() / abs(dmf_integrate(mesh, rho(:)))
      rho(1:mesh%np) = rr * rho(1:mesh%np)
 
    class default
      assert(.false.)
    end select
 
    call profiling_out("SPECIES_LR_DENSITY")
    pop_sub(species_get_long_range_density)
  end subroutine species_get_long_range_density
 
 
  ! ---------------------------------------------------------
  subroutine func(xin, ff, jacobian)
    real(real64), intent(in)  :: xin(:)
    real(real64), intent(out) :: ff(:), jacobian(:,:)
 
    real(real64), allocatable :: xrho(:)
    integer :: idir, jdir, dim, ip
 
    push_sub(func)
 
    dim = mesh_p%box%dim
 
    call getrho(dim, xin)
    safe_allocate(xrho(1:mesh_p%np))
 
    ! First, we calculate the function ff.
    do idir = 1, dim
      !$omp parallel do simd
      do ip = 1, mesh_p%np
        xrho(ip) = rho_p(ip) * mesh_p%x(ip, idir)
      end do
      ff(idir) = dmf_integrate(mesh_p, xrho) - pos_p(idir)
    end do
 
    ! Now the jacobian.
    do idir = 1, dim
      do jdir = 1, dim
        !$omp parallel do simd
        do ip = 1, mesh_p%np
          xrho(ip) = grho_p(ip, jdir) * mesh_p%x(ip, idir)
        end do
        jacobian(idir, jdir) = dmf_integrate(mesh_p, xrho)
      end do
    end do
 
    safe_deallocate_a(xrho)
    pop_sub(func)
  end subroutine func
 
  ! ---------------------------------------------------------
  subroutine species_get_nlcc(species, space, latt, pos, mesh, rho_core, accumulate)
    class(species_t), target, intent(in)   :: species
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    type(mesh_t),            intent(in)    :: mesh
    real(real64),            intent(inout) :: rho_core(:)
    logical, optional,       intent(in)    :: accumulate
 
    real(real64) :: center(space%dim), rr
    integer :: icell, ip
    type(lattice_iterator_t) :: latt_iter
    type(ps_t), pointer :: ps
 
    push_sub(species_get_nlcc)
 
    ! only for 3D pseudopotentials, please
    select type(species)
    class is(pseudopotential_t)
      ps => species%ps
      if (.not. optional_default(accumulate, .false.)) rho_core = m_zero
 
      latt_iter = lattice_iterator_t(latt, ps%core%x_threshold)
      do icell = 1, latt_iter%n_cells
        center = pos + latt_iter%get(icell)
        do ip = 1, mesh%np
          rr = norm2(mesh%x(ip, 1:space%dim) - center)
          if (rr < spline_range_max(ps%core)) then
            rho_core(ip) = rho_core(ip) + spline_eval(ps%core, rr)
          end if
        end do
      end do
    class default
      if (.not. optional_default(accumulate, .false.)) rho_core = m_zero
    end select
 
    pop_sub(species_get_nlcc)
  end subroutine species_get_nlcc
 
  ! ---------------------------------------------------------
  subroutine species_get_nlcc_grad(species, space, latt, pos, mesh, rho_core_grad, gnlcc_x)
    class(species_t), target, intent(in)   :: species
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    class(mesh_t),           intent(in)    :: mesh
    real(real64),            intent(out)   :: rho_core_grad(:,:)
    real(real64), optional,         intent(inout) :: gnlcc_x(:,:,:)
 
    real(real64) :: center(space%dim), rr, spline
    integer :: icell, ip, idir, jdir
    type(lattice_iterator_t) :: latt_iter
    type(ps_t), pointer :: ps
 
    push_sub(species_get_nlcc_grad)
 
    rho_core_grad = m_zero
 
    ! only for 3D pseudopotentials, please
    if (.not. species%is_ps()) then
      pop_sub(species_get_nlcc_grad)
      return
    endif
 
    select type(species)
    class is(pseudopotential_t)
      ps => species%ps
      if (.not. ps_has_nlcc(ps)) then
        pop_sub(species_get_nlcc_grad)
        return
      endif
 
      latt_iter = lattice_iterator_t(latt, ps%core_der%x_threshold)
      ! TODO: (#706) These loops should be reformulated as the code here is most likely very slow and inefficient
      do icell = 1, latt_iter%n_cells
        center = pos + latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin = center)
          rr = max(rr, r_small)
          if (rr >= spline_range_max(ps%core_der)) cycle
          spline = spline_eval(ps%core_der, rr)
          if(abs(spline) < 1e-150_real64) cycle
 
          do idir = 1, space%dim
            rho_core_grad(ip, idir) = rho_core_grad(ip, idir) - spline*(mesh%x(ip, idir)-center(idir))/rr
            if (present(gnlcc_x)) then
              do jdir = 1, space%dim
                gnlcc_x(ip, idir, jdir) = gnlcc_x(ip, idir, jdir) &
                  - spline*(mesh%x(ip, idir)-center(idir))/rr*(mesh%x(ip, jdir)-center(jdir))
              end do
            end if
 
          end do
        end do
      end do
    class default
      assert(.false.)
    end select
 
    pop_sub(species_get_nlcc_grad)
  end subroutine species_get_nlcc_grad
 
  ! ---------------------------------------------------------
  ! Return the density of a normalized Gaussian centered on xin
  ! as well as its gradient with respect to the central position
  subroutine getrho(dim, xin)
    integer, intent(in) :: dim
    real(real64),   intent(in) :: xin(1:dim)
 
    integer :: ip, idir
    real(real64)   :: r2, chi(dim), norm, threshold
 
    push_sub(getrho)
 
    ! We set here a threshold of 0.001 for the tail of the Gaussian, similar to what we do for the
    ! pseudopotentials. The value of the threshold corresponds to the default for pseudopotentials
    threshold = -log(0.001_real64)*alpha2_p
 
    do ip = 1, mesh_p%np
      ! This is not correct for curvilinear meshes
      chi(1:dim) = mesh_p%x(ip,1:dim)
      r2 = sum((chi - xin(1:dim))**2)
 
      if (r2 < threshold) then
        rho_p(ip) = exp(-r2/alpha2_p)
      else
        rho_p(ip) = m_zero
      end if
 
      do idir = 1, dim
        grho_p(ip, idir) = (chi(idir) - xin(idir)) * rho_p(ip)
      end do
    end do
 
    norm = dmf_integrate(mesh_p, rho_p)
    call lalg_scal(mesh_p%np, m_one/norm, rho_p)
    call lalg_scal(mesh_p%np, dim, m_two/alpha2_p/norm, grho_p)
 
    pop_sub(getrho)
  end subroutine getrho
 
 
  ! ---------------------------------------------------------
  subroutine species_get_local(species, namespace, space, latt, pos, mesh, vl)
    class(species_t), target, intent(in)   :: species
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(lattice_vectors_t), intent(in)    :: latt
    real(real64),            intent(in)    :: pos(1:space%dim)
    type(mesh_t),            intent(in)    :: mesh
    real(real64),            intent(out)   :: vl(:)
 
    real(real64) :: a1, a2, rb2, range, density ! for jellium
    real(real64) :: xx(space%dim), pos_pc(space%dim), r, r2, threshold
    integer :: ip, err, icell
    type(ps_t), pointer :: ps
    complex(real64) :: zpot
    type(lattice_iterator_t) :: latt_iter
    real(real64) :: aa, bb
 
 
    push_sub(species_get_local)
 
    call profiling_in("SPECIES_GET_LOCAL")
 
    select type(species)
 
    type is (soft_coulomb_t)
 
      call parse_variable(namespace, 'SpeciesProjectorSphereThreshold', 0.001_real64, threshold)
 
      !Assuming that we want to take the contribution from all replica that contributes up to 0.001
      ! to the center of the cell, we arrive to a range of 1000 a.u..
      latt_iter = lattice_iterator_t(latt, species%get_zval() / threshold)
      vl = m_zero
      do icell = 1, latt_iter%n_cells
        pos_pc = pos + latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, r, origin = pos_pc)
          r2 = r*r
          vl(ip) = vl(ip) -species%get_zval()/sqrt(r2+species%get_softening())
        end do
      end do
 
    type is (species_user_defined_t)
      !TODO: we should control the value of 5 by a variable.
      range = 5.0_real64 * latt%max_length()
      latt_iter = lattice_iterator_t(latt, range)
      vl = m_zero
      do icell = 1, latt_iter%n_cells
        pos_pc = pos + latt_iter%get(icell)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, r, origin = pos_pc, coords = xx)
 
          zpot = species%user_pot(space%dim, xx, r)
          vl(ip) = vl(ip) + real(zpot, real64)
        end do
      end do
 
    type is(species_from_file_t)
 
      assert(.not. space%is_periodic())
 
      call dio_function_input(trim(species%get_filename()), namespace, space, mesh, vl, err)
      if (err /= 0) then
        write(message(1), '(a)')    'Error loading file '//trim(species%get_filename())//'.'
        write(message(2), '(a,i4)') 'Error code returned = ', err
        call messages_fatal(2, namespace=namespace)
      end if
 
    type is(jellium_sphere_t)
 
      assert(.not. space%is_periodic())
 
      a1 = species%get_z()/(m_two*species%radius()**3)
      a2 = species%get_z()/species%radius()
      rb2= species%radius()**2
 
      do ip = 1, mesh%np
 
        xx = mesh%x(ip, :) - pos(1:space%dim)
        r = norm2(xx)
 
        if (r <= species%radius()) then
          vl(ip) = (a1*(r*r - rb2) - a2)
        else
          vl(ip) = -species%get_z()/r
        end if
 
      end do
 
    type is (jellium_slab_t)
 
      ! Electrostatic potential from an infinite slab of thickness species%thickness
      ! Potential and electric fields are continuous at +/- L/2
      density = species%get_density(mesh%box%bounding_box_l)
      a1 = m_four * m_pi * density * species%thickness() / m_two
 
      do ip = 1, mesh%np
 
        r = abs(mesh%x(ip, 3) - pos(3))
 
        if (r <= species%thickness()/m_two) then
          vl(ip) = a1 * (r * r / species%thickness() + species%thickness() / m_four)
        else
          vl(ip) = a1 * r
        end if
 
      end do
 
    class is (pseudopotential_t)
 
      assert(.not. space%is_periodic())
 
      ps => species%ps
 
      do ip = 1, mesh%np
        r2 = sum((mesh%x(ip, :) - pos)**2)
        if (r2 < spline_range_max(ps%vlr_sq)) then
          vl(ip) = spline_eval(ps%vlr_sq, r2)
        else
          vl(ip) = p_proton_charge*species%get_zval()/sqrt(r2)
        end if
      end do
 
      nullify(ps)
 
    type is (full_anc_t)
      ! periodic copies are not considered in this routine
      if (space%is_periodic()) then
        call messages_experimental("species_full_anc for periodic systems", namespace=namespace)
      end if
 
      aa = species%a()
      bb = species%b()
      assert(bb < m_zero) ! To be sure it was computed
 
      ! Evaluation of the scaled potential, see Eq. 19
      do ip = 1, mesh%np
        r2 = sum((mesh%x(ip, :) - pos)**2)*(species%get_z()*aa)**2
        if(r2 > r_small**2) then
          r = sqrt(r2)
          vl(ip) = -m_half &
            - (loct_erf(r) + m_two*(aa*bb + m_one/sqrt(m_pi))*r*exp(-r2))/r*aa &
            + (loct_erf(r) + m_two*(aa*bb + m_one/sqrt(m_pi))*r*exp(-r2))**2*m_half &
            + (-m_two*aa**2*bb - m_four*aa/sqrt(m_pi) &
            + m_four*aa*(aa*bb + m_one/sqrt(m_pi))*r2)*exp(-r2)*m_half
        else ! Eq. 10
          vl(ip) = -m_half - m_three * aa**2*bb - 6.0_real64*aa/sqrt(m_pi)
        end if
        vl(ip) = vl(ip) * (species%get_z())**2
      end do
 
    class default
      vl(1:mesh%np) = m_zero
    end select
 
    call profiling_out("SPECIES_GET_LOCAL")
    pop_sub(species_get_local)
  end subroutine species_get_local
 
end module species_pot_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: