external_potential.F90

Go to the documentation of this file.

!! Copyright (C) 2020 N. 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 external_potential_oct_m
  use debug_oct_m
  use global_oct_m
  use iihash_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use interactions_factory_oct_m
  use interaction_partner_oct_m
  use io_function_oct_m
  use iteration_counter_oct_m
  use lorentz_force_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  implicit none
 
  private
  public ::               &
    external_potential_t, &
    external_potential_clone, &
    load_external_potentials
 
  type, extends(interaction_partner_t) :: external_potential_t
    private
 
    type(space_t)   :: space
 
    integer, public :: type
 
    character(len=1024) :: potential_formula
    character(len=200)  :: density_formula
    character(len=MAX_PATH_LEN) :: filename
    real(real64) :: omega
 
    real(real64), allocatable, public :: pot(:)
 
    real(real64), allocatable, public :: b_field(:)
    integer                    :: gauge_2D
    real(real64), allocatable, public :: a_static(:,:)
    real(real64), allocatable, public :: e_field(:)
    !Auxiliary arrays for the electrons only
    !TODO: Suppress once electrons fully use the new framework
    real(real64), allocatable, public :: v_ext(:)
 
  contains
    procedure :: calculate => external_potential_calculate
    procedure :: allocate_memory => external_potential_allocate
    procedure :: deallocate_memory => external_potential_deallocate
    procedure :: init_interaction_as_partner => external_potential_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => external_potential_copy_quantities_to_interaction
    final :: external_potential_finalize
  end type external_potential_t
 
  integer, public, parameter ::  &
    EXTERNAL_POT_USDEF          = 201,           & !< user-defined function for local potential
    external_pot_from_file      = 202,           & 
    external_pot_charge_density = 203,           & 
    external_pot_static_bfield  = 204,           & 
    external_pot_static_efield  = 205              
 
 
  interface external_potential_t
    module procedure external_potential_init
  end interface external_potential_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine external_potential_clone(pot_out, pot_in)
    class(external_potential_t), pointer, intent(out) :: pot_out
    class(external_potential_t),         intent(in)  :: pot_in
 
    type(quantity_iterator_t) :: iter
    class(quantity_t), pointer :: quantity
 
    push_sub(external_potential_clone)
 
    allocate(pot_out)
 
    pot_out%namespace = pot_in%namespace
    pot_out%space = pot_in%space
    pot_out%type = pot_in%type
    pot_out%potential_formula = pot_in%potential_formula
    pot_out%density_formula = pot_in%density_formula
    pot_out%filename = pot_in%filename
    pot_out%omega = pot_in%omega
    pot_out%gauge_2D = pot_in%gauge_2D
 
    if (allocated(pot_in%pot)) then
      allocate(pot_out%pot, source=pot_in%pot)
    end if
    if (allocated(pot_in%b_field)) then
      allocate(pot_out%b_field, source=pot_in%b_field)
    end if
    if (allocated(pot_in%a_static)) then
      allocate(pot_out%a_static, source=pot_in%a_static)
    end if
    if (allocated(pot_in%e_field)) then
      allocate(pot_out%e_field, source=pot_in%e_field)
    end if
    if (allocated(pot_in%v_ext)) then
      allocate(pot_out%v_ext, source=pot_in%v_ext)
    end if
 
    if (allocated(pot_in%supported_interactions_as_partner)) then
      allocate(pot_out%supported_interactions_as_partner, source=pot_in%supported_interactions_as_partner)
    else
      allocate(pot_out%supported_interactions_as_partner(0))
    end if
 
    call iter%start(pot_in%quantities)
    do while (iter%has_next())
      quantity => iter%get_next()
      call pot_out%quantities%add(quantity)
    end do
 
    pop_sub(external_potential_clone)
  end subroutine external_potential_clone
 
  function external_potential_init(namespace) result(this)
    class(external_potential_t), pointer :: this
    type(namespace_t), intent(in) :: namespace
 
    push_sub(external_potential_init)
 
    allocate(this)
 
    this%namespace = namespace_t("ExternalPotential", parent=namespace)
    this%space = space_t(namespace)
 
    this%type = -1
 
    allocate(this%supported_interactions_as_partner(0))
 
    call this%quantities%add(quantity_t("E field", always_available = .true., updated_on_demand = .false., &
      iteration = iteration_counter_t()))
    call this%quantities%add(quantity_t("B field", always_available = .true., updated_on_demand = .false., &
      iteration = iteration_counter_t()))
 
    pop_sub(external_potential_init)
  end function external_potential_init
 
  ! ---------------------------------------------------------
  subroutine external_potential_finalize(this)
    type(external_potential_t), intent(inout) :: this
 
    push_sub(external_potential_finalize)
 
    call this%deallocate_memory()
 
    pop_sub(external_potential_finalize)
  end subroutine external_potential_finalize
 
  ! ---------------------------------------------------------
  subroutine external_potential_allocate(this, mesh)
    class(external_potential_t), intent(inout) :: this
    class(mesh_t),               intent(in)    :: mesh
 
    push_sub(external_potential_allocate)
 
    select case (this%type)
    case (external_pot_usdef, external_pot_from_file, external_pot_charge_density)
      safe_allocate(this%pot(1:mesh%np))
    case (external_pot_static_bfield)
      safe_allocate(this%a_static(1:mesh%np, 1:this%space%dim))
    case (external_pot_static_efield)
      if (this%space%periodic_dim < this%space%dim) then
        safe_allocate(this%pot(1:mesh%np))
        safe_allocate(this%v_ext(1:mesh%np_part))
      end if
    end select
 
    pop_sub(external_potential_allocate)
  end subroutine external_potential_allocate
 
  ! ---------------------------------------------------------
  subroutine external_potential_deallocate(this)
    class(external_potential_t), intent(inout) :: this
 
    push_sub(external_potential_deallocate)
 
    safe_deallocate_a(this%pot)
    safe_deallocate_a(this%b_field)
    safe_deallocate_a(this%a_static)
    safe_deallocate_a(this%e_field)
    safe_deallocate_a(this%v_ext)
 
    pop_sub(external_potential_deallocate)
  end subroutine external_potential_deallocate
 
  ! ---------------------------------------------------------
  subroutine external_potential_init_interaction_as_partner(partner, interaction)
    class(external_potential_t),     intent(in)    :: partner
    class(interaction_surrogate_t),  intent(inout) :: interaction
 
    push_sub(external_potential_init_interaction_as_partner)
 
    select type (interaction)
    type is (lorentz_force_t)
      ! Nothing to be initialized for the Lorentz force.
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(external_potential_init_interaction_as_partner)
  end subroutine external_potential_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine external_potential_copy_quantities_to_interaction(partner, interaction)
    class(external_potential_t),     intent(inout) :: partner
    class(interaction_surrogate_t),  intent(inout) :: interaction
 
    integer :: ip
 
    push_sub(external_potential_copy_quantities_to_interaction)
 
    select type (interaction)
    type is (lorentz_force_t)
      if (partner%type == external_pot_static_efield) then
        do ip = 1, interaction%system_np
          interaction%partner_e_field(:, ip) = partner%e_field
          interaction%partner_b_field(:, ip) = m_zero
        end do
      else if (partner%type == external_pot_static_bfield) then
        do ip = 1, interaction%system_np
          interaction%partner_e_field(:, ip) = m_zero
          interaction%partner_b_field(:, ip) = partner%b_field
        end do
      else
        assert(.false.) !This should never occur.
      end if
 
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(external_potential_copy_quantities_to_interaction)
  end subroutine external_potential_copy_quantities_to_interaction
 
 
  ! ---------------------------------------------------------
  subroutine external_potential_calculate(this, namespace, mesh, poisson)
    class(external_potential_t), intent(inout) :: this
    type(namespace_t),           intent(in)    :: namespace
    class(mesh_t),               intent(in)    :: mesh
    type(poisson_t),             intent(in)    :: poisson
 
    real(real64) :: pot_re, pot_im, r, xx(this%space%dim)
    real(real64), allocatable :: den(:), grx(:)
    integer :: ip, err
 
    push_sub(external_potential_calculate)
 
    select case (this%type)
 
    case (external_pot_usdef)
      assert(allocated(this%pot))
 
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, r, coords = xx)
        call parse_expression(pot_re, pot_im, this%space%dim, xx, r, m_zero, this%potential_formula)
        this%pot(ip) = pot_re
      end do
 
    case (external_pot_from_file)
      assert(allocated(this%pot))
 
      call dio_function_input(trim(this%filename), namespace, this%space, mesh, this%pot, err)
      if (err /= 0) then
        write(message(1), '(a)')    'Error loading file '//trim(this%filename)//'.'
        write(message(2), '(a,i4)') 'Error code returned = ', err
        call messages_fatal(2, namespace=namespace)
      end if
 
    case (external_pot_charge_density)
      assert(allocated(this%pot))
 
      safe_allocate(den(1:mesh%np))
 
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, r, coords = xx)
        call parse_expression(pot_re, pot_im, this%space%dim, xx, r, m_zero, this%potential_formula)
        den(ip) = pot_re
      end do
 
      call dpoisson_solve(poisson, namespace, this%pot, den, all_nodes = .false.)
 
      safe_deallocate_a(den)
 
    case (external_pot_static_bfield)
      assert(allocated(this%b_field))
 
      ! Compute the vector potential from a uniform B field.
      ! The sign is determined by the relation $\vec{B} = \nabla \times \vec{A}$.
      ! This leads to  $\vec{A} = -\frac{1}{2}\vec{r}\times\vec{B}$.
      ! A factor 1/c is already added, to avoid adding it everytime we update the Hamiltonian
      safe_allocate(grx(1:this%space%dim))
 
      select case (this%space%dim)
      case (2)
        select case (this%gauge_2d)
        case (0) ! linear_xy
          if (this%space%periodic_dim == 1) then
            message(1) = "For 2D system, 1D-periodic, StaticMagneticField can only be "
            message(2) = "applied for StaticMagneticField2DGauge = linear_y."
            call messages_fatal(2, namespace=namespace)
          end if
          !$omp parallel do private(grx)
          do ip = 1, mesh%np
            grx(1:this%space%dim) = mesh%x(1:this%space%dim, ip)
            this%a_static(ip, 1) = -m_half/p_c * grx(2) * this%b_field(3)
            this%a_static(ip, 2) =  m_half/p_c * grx(1) * this%b_field(3)
          end do
        case (1) ! linear y
          !$omp parallel do private(grx)
          do ip = 1, mesh%np
            grx(1:this%space%dim) = mesh%x(1:this%space%dim, ip)
            this%a_static(ip, 1) = -m_one/p_c * grx(2) * this%b_field(3)
            this%a_static(ip, 2) = m_zero
          end do
        end select
      case (3)
        !$omp parallel do private(grx)
        do ip = 1, mesh%np
          grx(1:this%space%dim) = mesh%x(1:this%space%dim, ip)
          this%a_static(ip, 1) = -m_half/p_c*(grx(2) * this%b_field(3) - grx(3) * this%b_field(2))
          this%a_static(ip, 2) = -m_half/p_c*(grx(3) * this%b_field(1) - grx(1) * this%b_field(3))
          this%a_static(ip, 3) = -m_half/p_c*(grx(1) * this%b_field(2) - grx(2) * this%b_field(1))
        end do
      case default
        assert(.false.)
      end select
 
      safe_deallocate_a(grx)
 
    case (external_pot_static_efield)
      assert(allocated(this%e_field))
 
      if (this%space%periodic_dim < this%space%dim) then
        ! Compute the scalar potential
        !
        ! Note that the -1 sign is missing. This is because we
        ! consider the electrons with +1 charge. The electric field
        ! however retains the sign because we also consider protons to
        ! have +1 charge when calculating the force.
        !
        ! NTD: This comment is very confusing and prone to error
        ! TODO: Fix this to have physically sound quantities and interactions
        do ip = 1, mesh%np
          this%pot(ip) = sum(mesh%x(this%space%periodic_dim + 1:this%space%dim, ip) &
            * this%e_field(this%space%periodic_dim + 1:this%space%dim))
        end do
        ! The following is needed to make interpolations.
        ! It is used by PCM.
        this%v_ext(1:mesh%np) = this%pot(1:mesh%np)
        do ip = mesh%np+1, mesh%np_part
          this%v_ext(ip) = sum(mesh%x(this%space%periodic_dim + 1:this%space%dim, ip) &
            * this%e_field(this%space%periodic_dim + 1:this%space%dim))
        end do
      end if
 
    end select
 
    pop_sub(external_potential_calculate)
  end subroutine external_potential_calculate
 
  subroutine load_external_potentials(external_potentials, namespace)
    class(partner_list_t), intent(inout)  :: external_potentials
    type(namespace_t),    intent(in)     :: namespace
 
    integer :: n_pot_block, row, read_data
    type(block_t) :: blk
    class(external_potential_t), pointer :: pot
 
    integer :: dim, periodic_dim, idir
 
    push_sub(load_external_potentials)
 
    !%Variable StaticExternalPotentials
    !%Type block
    !%Section System
    !%Description
    !% An static external potential is a model potential added to the local potential of the Hamiltonian
    !%
    !% The format of this block is the following:
    !% The first field defines the type of species (the valid options are detailed
    !% below).
    !%
    !% Then a list of parameters follows. The parameters are specified
    !% by a first field with the parameter name and the field that
    !% follows with the value of the parameter. Some parameters are
    !% specific to a certain species while others are accepted by all
    !% species. These are <tt>mass</tt>, <tt>max_spacing</tt>, and <tt>min_radius</tt>.
    !%
    !% These are examples of possible species:
    !%
    !% <tt>%ExternalPotential
    !% <br>&nbsp;&nbsp; potential_user_defined | potential_formula | "1/2*r^2"
    !% <br>%</tt>
    !%Option potential_from_file  -202
    !% The potential is read from a file. Accepted file formats, detected by extension: obf, ncdf and csv.
    !%Option potential_user_defined -201
    !% Species with user-defined potential. The potential for the
    !% species is defined by the formula given by the <tt>potential_formula</tt>
    !% parameter.
    !%Option potential_charge_density -203
    !% The potential for this species is created from the distribution
    !% of charge given by the <tt>density_formula</tt> parameter.
    !%Option file -10010
    !% The path for the file that describes the species.
    !%Option potential_formula -10012
    !% Mathematical expression that defines the potential for <tt>species_user_defined</tt>. You can use
    !% any of the <i>x</i>, <i>y</i>, <i>z</i> or <i>r</i> variables.
    !%Option density_formula -10013
    !% Mathematical expression that defines the charge density for <tt>species_charge_density</tt>. You can use
    !% any of the <i>x</i>, <i>y</i>, <i>z</i> or <i>r</i> variables.
    !%End
 
    ! First, find out if there is a Species block.
    n_pot_block = 0
    if (parse_block(namespace, 'StaticExternalPotentials', blk) == 0) then
      n_pot_block = parse_block_n(blk)
 
      do row = 0, n_pot_block-1
        !Create a potential
        pot => external_potential_t(namespace)
        !Parse the information from the block
        call read_from_block(pot, namespace, blk, row, read_data)
        assert(read_data > 0)
        !Add this to the list
        call external_potentials%add(pot)
      end do
      call parse_block_end(blk)
    end if
 
 
    !Here I am parsing the variables Dimensions et PeriodicDimensions because we do not have access
    !to this information here.
    !TODO: This needs to be removed and replaced by something better
    call parse_variable(namespace, 'Dimensions', 3, dim)
    call parse_variable(namespace, 'PeriodicDimensions', 0, periodic_dim)
 
 
    !%Variable StaticMagneticField
    !%Type block
    !%Section Hamiltonian
    !%Description
    !% A static constant magnetic field may be added to the usual Hamiltonian,
    !% by setting the block <tt>StaticMagneticField</tt>.
    !% The three possible components of the block (which should only have one
    !% line) are the three components of the magnetic field vector. Note that
    !% if you are running the code in 1D mode, this will not work, and if you
    !% are running the code in 2D mode the magnetic field will have to be in
    !% the <i>z</i>-direction, so that the first two columns should be zero.
    !% Possible in periodic system only in these cases: 2D system, 1D periodic,
    !% with <tt>StaticMagneticField2DGauge = linear_y</tt>;
    !% 3D system, 1D periodic, field is zero in <i>y</i>- and <i>z</i>-directions (given
    !% currently implemented gauges).
    !%
    !% The magnetic field should always be entered in atomic units, regardless
    !% of the <tt>Units</tt> variable. Note that we use the "Gaussian" system
    !% meaning 1 au[B] = <math> 2.350517568\times 10^9</math> Gauss, which corresponds to
    !% <math>2.3505175678\times 10^5</math> Tesla.
    !%End
    if (parse_block(namespace, 'StaticMagneticField', blk) == 0) then
      !Create a potential
      pot => external_potential_t(namespace)
      pot%type = external_pot_static_bfield
      pot%supported_interactions_as_partner = [lorentz_force]
 
      !%Variable StaticMagneticField2DGauge
      !%Type integer
      !%Default linear_xy
      !%Section Hamiltonian
      !%Description
      !% The gauge of the static vector potential <math>A</math> when a magnetic field
      !% <math>B = \left( 0, 0, B_z \right)</math> is applied to a 2D-system.
      !%Option linear_xy 0
      !% Linear gauge with <math>A = \frac{1}{2c} \left( -y, x \right) B_z</math>. (Cannot be used for periodic systems.)
      !%Option linear_y 1
      !% Linear gauge with <math>A = \frac{1}{c} \left( -y, 0 \right) B_z</math>. Can be used for <tt>PeriodicDimensions = 1</tt>
      !% but not <tt>PeriodicDimensions = 2</tt>.
      !%End
      call parse_variable(namespace, 'StaticMagneticField2DGauge', 0, pot%gauge_2d)
      if (.not. varinfo_valid_option('StaticMagneticField2DGauge', pot%gauge_2d)) then
        call messages_input_error(namespace, 'StaticMagneticField2DGauge')
      end if
 
      safe_allocate(pot%b_field(1:3))
      do idir = 1, 3
        call parse_block_float(blk, 0, idir - 1, pot%b_field(idir))
      end do
      select case (dim)
      case (1)
        call messages_input_error(namespace, 'StaticMagneticField')
      case (2)
        if (periodic_dim == 2) then
          message(1) = "StaticMagneticField cannot be applied in a 2D, 2D-periodic system."
          call messages_fatal(1, namespace=namespace)
        end if
        if (pot%b_field(1)**2 + pot%b_field(2)**2 > m_zero) then
          call messages_input_error(namespace, 'StaticMagneticField')
        end if
      case (3)
        ! Consider cross-product below: if grx(1:this%space%periodic_dim) is used, it is not ok.
        ! Therefore, if idir is periodic, b_field for all other directions must be zero.
        ! 1D-periodic: only Bx. 2D-periodic or 3D-periodic: not allowed. Other gauges could allow 2D-periodic case.
        if (periodic_dim >= 2) then
          message(1) = "In 3D, StaticMagneticField cannot be applied when the system is 2D- or 3D-periodic."
          call messages_fatal(1, namespace=namespace)
        else if (periodic_dim == 1 .and. any(abs(pot%b_field(2:3)) > m_zero)) then
          message(1) = "In 3D, 1D-periodic, StaticMagneticField must be zero in the y- and z-directions."
          call messages_fatal(1, namespace=namespace)
        end if
      end select
      call parse_block_end(blk)
 
      if (dim > 3) call messages_not_implemented('Magnetic field for dim > 3', namespace=namespace)
 
      !Add this to the list
      call external_potentials%add(pot)
 
      !The corresponding A field on the mesh is computed in the routine external_potential_calculate
 
    end if
 
    !%Variable StaticElectricField
    !%Type block
    !%Default 0
    !%Section Hamiltonian
    !%Description
    !% A static constant electric field may be added to the usual Hamiltonian,
    !% by setting the block <tt>StaticElectricField</tt>.
    !% The three possible components of the block (which should only have one
    !% line) are the three components of the electric field vector.
    !% It can be applied in a periodic direction of a large supercell via
    !% the single-point Berry phase.
    !%End
    if (parse_block(namespace, 'StaticElectricField', blk) == 0) then
      !Create a potential
      pot => external_potential_t(namespace)
      pot%type = external_pot_static_efield
      pot%supported_interactions_as_partner = [lorentz_force]
 
      safe_allocate(pot%e_field(1:dim))
      do idir = 1, dim
        call parse_block_float(blk, 0, idir - 1, pot%e_field(idir), units_inp%energy / units_inp%length)
 
        !Electron-specific checks (k-points) are done in the hamiltonian_elec.F90 file
        if (idir <= periodic_dim .and. abs(pot%e_field(idir)) > m_epsilon) then
          message(1) = "Applying StaticElectricField in a periodic direction is only accurate for large supercells."
          call messages_warning(1, namespace=namespace)
        end if
      end do
      call parse_block_end(blk)
 
      !Add this to the list
      call external_potentials%add(pot)
 
      !The corresponding A field on the mesh is computed in the routine external_potential_calculate
    end if
 
 
    pop_sub(load_external_potentials)
  end subroutine load_external_potentials
 
  ! ---------------------------------------------------------
  subroutine read_from_block(pot, namespace, blk, row, read_data)
    type(external_potential_t), intent(inout) :: pot
    type(namespace_t),          intent(in)    :: namespace
    type(block_t),              intent(in)    :: blk
    integer,                    intent(in)    :: row
    integer,                    intent(out)   :: read_data
 
    integer :: ncols, icol, flag
    type(iihash_t) :: read_parameters
 
 
    push_sub(read_from_block)
 
    ncols = parse_block_cols(blk, row)
    read_data = 0
 
    call parse_block_integer(blk, row, 0, pot%type)
 
    ! To detect the old species block format, options are represented
    ! as negative values. If we get a non-negative value we know we
    ! are reading a mass.
    if (pot%type >= 0) then
      message(1) = 'Error in reading the ExternalPotentials block'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! now we convert back to positive
    pot%type = -pot%type
 
    read_data = 1
 
    if (pot%type /= external_pot_charge_density .and. pot%type /= external_pot_usdef .and. pot%type /= external_pot_from_file) then
      call messages_input_error(namespace, 'ExternalPotentials', "Unknown type of external potential")
    end if
 
    call iihash_init(read_parameters)
 
    icol = read_data
    do
      if (icol >= ncols) exit
 
      call parse_block_integer(blk, row, icol, flag)
 
      select case (flag)
 
      case (option__staticexternalpotentials__file)
        call check_duplication(option__staticexternalpotentials__file)
        call parse_block_string(blk, row, icol + 1, pot%filename)
 
      case (option__staticexternalpotentials__potential_formula)
        call check_duplication(option__staticexternalpotentials__potential_formula)
        call parse_block_string(blk, row, icol + 1, pot%potential_formula, convert_to_c=.true.)
 
        if (pot%type /= external_pot_usdef) then
          call messages_input_error(namespace, 'ExternalPotentials', 'potential_formula can only be used with user_defined')
        end if
 
      case (option__staticexternalpotentials__density_formula)
        call check_duplication(option__staticexternalpotentials__density_formula)
        call parse_block_string(blk, row, icol + 1, pot%density_formula, convert_to_c=.true.)
 
        if (pot%type /= external_pot_charge_density) then
          call messages_input_error(namespace, 'ExternalPotentials', 'density_formula can only be used with charge_density')
        end if
 
      case default
        call messages_input_error(namespace, 'ExternalPotentials', "Unknown parameter ")
 
      end select
 
      icol = icol + 2
    end do
    ! CHECK THAT WHAT WE PARSED MAKES SENSE
 
 
    if (pot%type == external_pot_usdef .and. &
      .not. parameter_defined(option__staticexternalpotentials__potential_formula)) then
      call messages_input_error(namespace, 'ExternalPotentials', "The 'potential_formula' parameter is missing.")
    end if
 
    if (pot%type == external_pot_charge_density .and. &
      .not. parameter_defined(option__staticexternalpotentials__density_formula)) then
      call messages_input_error(namespace, 'ExternalPotentials', "The 'density_formula' parameter is missing.")
    end if
 
    if (pot%type == external_pot_from_file .and. .not. (parameter_defined(option__staticexternalpotentials__file))) then
      call messages_input_error(namespace, 'ExternalPotentials', "The 'file' parameter is missing.")
    end if
 
    call iihash_end(read_parameters)
 
    pop_sub(read_from_block)
 
  contains
 
    logical function parameter_defined(param) result(defined)
      integer(int64), intent(in) :: param
 
      integer :: tmp
 
      push_sub(read_from_block.parameter_defined)
 
      tmp = iihash_lookup(read_parameters, int(-param), defined)
 
      pop_sub(read_from_block.parameter_defined)
    end function parameter_defined
 
    !------------------------------------------------------
 
    subroutine check_duplication(param)
      integer(int64), intent(in) :: param
 
      push_sub(read_from_block.check_duplication)
 
      if (parameter_defined(param)) then
        call messages_input_error(namespace, 'ExternalPotentials', "Duplicated parameter in external potential.")
      end if
 
      call iihash_insert(read_parameters, int(-param), 1)
 
      pop_sub(read_from_block.check_duplication)
    end subroutine check_duplication
 
  end subroutine read_from_block
  ! ---------------------------------------------------------
 
 
end module external_potential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: