species_factory.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023-2024 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 species_factory_oct_m
  use allelectron_oct_m
  use debug_oct_m
  use element_oct_m
  use global_oct_m
  use iihash_oct_m
  use io_oct_m
  use jellium_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use pseudo_set_oct_m
  use pseudopotential_oct_m
  use species_oct_m
  use share_directory_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                        &
    species_factory_t,             &
    species_factory_init,          &
    species_factory_end
 
  type :: species_factory_t
    private
    logical :: initialized = .false.
    integer :: default_pseudopotential_set_id
    type(pseudo_set_t) :: default_pseudopotential_set
    integer :: default_allelectron_type
    real(real64) :: default_sigma
    real(real64) :: default_anc_a
  contains
    procedure :: init => species_factory_init
    procedure :: end => species_factory_end
    procedure :: create_from_input => species_factory_create_from_input
    procedure :: create_from_block => read_from_block
  end type species_factory_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine species_factory_init(factory, namespace)
    class(species_factory_t),   intent(inout) :: factory
    type(namespace_t),          intent(in)    :: namespace
 
    integer :: ierr, default_val
 
    if (factory%initialized) return
 
    push_sub_with_profile(species_factory_init)
 
    factory%initialized = .true.
 
    call share_directory_set(conf%share)
 
    !%Variable AllElectronType
    !%Type integer
    !%Default no
    !%Section System::Species
    !%Description
    !% Selects the type of all-electron species that applies by default to all
    !% atoms. This is not compatible with <tt>PseudopotentialSet</tt>, but it is
    !% compatible with the <tt>Species</tt> block.
    !%
    !% The three types differ in what they support:
    !% <tt>full_delta</tt> and <tt>full_anc</tt> work with curvilinear coordinates
    !% (see <tt>CurvMethod</tt>), <tt>full_gaussian</tt> works with <tt>curv_gygi</tt> and
    !% <tt>curv_modine</tt> but not with <tt>curv_briggs</tt>. For periodic systems,
    !% <tt>full_gaussian</tt> is not implemented and <tt>full_anc</tt> is experimental.
    !%
    !% The relevant parameter is <tt>AllElectronSigma</tt> for <tt>full_gaussian</tt> and
    !% <tt>AllElectronANCParam</tt> for <tt>full_anc</tt>; <tt>full_delta</tt> takes neither.
    !% A parameter set for a type that does not use it is ignored, and a warning is printed.
    !%
    !%Option no 0
    !% Do not specify any default all-electron type of species. All species must be
    !% specified in the Species block.
    !%Option full_delta 1
    !% All atoms are supposed to be by default of type <tt>species_full_delta</tt>.
    !%Option full_gaussian 2
    !% All atoms are supposed to be by default of type <tt>species_full_gaussian</tt>.
    !%Option full_anc 3
    !% All atoms are supposed to be by default of type <tt>species_full_anc</tt>.
    !%End
    call parse_variable(namespace, 'AllElectronType', option__allelectrontype__no, factory%default_allelectron_type)
    call messages_print_var_option('AllElectronType', factory%default_allelectron_type, namespace=namespace)
 
    !%Variable AllElectronSigma
    !%Type float
    !%Default 0.6
    !%Section System::Species
    !%Description
    !% Default value for the parameter <tt>gaussian_width</tt>. This is useful
    !% for specifying multiple atoms without specifying the species block. The
    !% default value is taken from the recommendation in
    !% <i>Phys. Rev. B</i> <b>55</b>, 10289 (1997).
    !%
    !% Note that this is a width in units of the grid spacing, not a length: the
    !% Gaussian exponent is <math>\zeta = 1/(2\sigma^2\delta^2)</math> with
    !% <math>\delta</math> the spacing. Reducing <tt>Spacing</tt> at fixed
    !% <math>\sigma</math> therefore shrinks the physical size of the nucleus, so a
    !% spacing convergence study is not a convergence study of a single Hamiltonian.
    !%
    !% It is only used by <tt>species_full_gaussian</tt> and is ignored for the other
    !% all-electron types.
    !%End
    call parse_variable(namespace, 'AllElectronSigma', 0.6_real64, factory%default_sigma)
 
    !%Variable AllElectronANCParam
    !%Type float
    !%Default 4.0
    !%Section System::Species
    !%Description
    !% Default value for the parameter <tt>anc_a</tt> of the <tt>species_full_anc</tt>
    !% species, the parameter <i>a</i> of the analytic norm-conserving regularized Coulomb
    !% potential of [Gygi, <i>J. Chem. Theory Comput.</i> <b>19</b>, 1300 (2023)]. This is
    !% useful for specifying multiple atoms without specifying the species block.
    !%
    !% It has the dimension of an inverse length, and the code uses the scaled radius
    !% <math>a Z r</math>, so the regularization length is <math>1/(aZ)</math> bohr and
    !% depends on the atomic number. The default value of 4 has only been validated for
    !% hydrogen; for a heavy atom it puts the regularization far below any usable grid
    !% spacing, and a smaller value should be used.
    !%
    !% It is only used by <tt>species_full_anc</tt> and is ignored for the other
    !% all-electron types.
    !%End
    call parse_variable(namespace, 'AllElectronANCParam', 4.0_real64, factory%default_anc_a)
 
    ! These are valid variable names, so the parser cannot report them as unused.
    if (parse_is_defined(namespace, 'AllElectronSigma') .and. &
      factory%default_allelectron_type /= option__allelectrontype__full_gaussian) then
      message(1) = "AllElectronSigma is only used by species_full_gaussian and will be ignored."
      call messages_warning(1, namespace=namespace)
    end if
    if (parse_is_defined(namespace, 'AllElectronANCParam') .and. &
      factory%default_allelectron_type /= option__allelectrontype__full_anc) then
      message(1) = "AllElectronANCParam is only used by species_full_anc and will be ignored."
      call messages_warning(1, namespace=namespace)
    end if
 
    !%Variable PseudopotentialSet
    !%Type integer
    !%Default standard
    !%Section System::Species
    !%Description
    !% Selects the set of pseudopotentials used by default for species
    !% not defined in the <tt>Species</tt> block.
    !%
    !% These sets of pseudopotentials come from different
    !% sources. Octopus developers have not validated them. We include
    !% them with the code for convenience of the users, but you are
    !% expected to check the quality and suitability of the
    !% pseudopotential for your application.
    !%
    !% Note that not all these pseudopotentials are compatible with spin-orbit coupling
    !% Only hgh_lda, hgh_lda_sc and pseudodojo_pbe_fr provide SOC information.
    !%
    !%Option none 0
    !% Do not load any pseudopotential by default. All species must be
    !% specified in the Species block.
    !%Option standard 1
    !% The standard set of Octopus that provides LDA pseudopotentials
    !% in the PSF format for some elements: H, Li, C, N, O, Na, Si, S, Ti, Se, Cd.
    !%Option sg15 2
    !% The set of Optimized Norm-Conserving Vanderbilt
    !% PBE pseudopotentials. Ref: M. Schlipf and F. Gygi, <i>Comp. Phys. Commun.</i> <b>196</b>, 36 (2015).
    !% This set provides pseudopotentials for elements up to Z = 83
    !% (Bi), excluding Lanthanides.
    !% Current version of the set is 1.2.
    !%Option hgh_lda 3
    !% The set of Hartwigsen-Goedecker-Hutter LDA pseudopotentials for elements from H to Rn.
    !% Ref: C. Hartwigsen, S. Goedecker, and J. Hutter, <i>Phys. Rev. B</i> <b>58</b>, 3641 (1998).
    !%Option hgh_lda_sc 31
    !% The semicore set of Hartwigsen-Goedecker-Hutter LDA pseudopotentials.
    !% Ref: C. Hartwigsen, S. Goedecker, and J. Hutter, <i>Phys. Rev. B</i> <b>58</b>, 3641 (1998).
    !%Option hscv_lda 4
    !% The set of Hamann-Schlueter-Chiang-Vanderbilt (HSCV) potentials
    !% for LDA exchange and correlation downloaded from http://fpmd.ucdavis.edu/potentials/index.htm.
    !% These pseudopotentials were originally intended for the QBox
    !% code. They were generated using the method of Hamann, Schluter
    !% and Chiang. Ref: D. Vanderbilt, <i>Phys. Rev. B</i> <b>32</b>, 8412 (1985).
    !% Warning from the original site: The potentials provided in this
    !% site are distributed without warranty. In most cases,
    !% potentials were not tested. Potentials should be thoroughly
    !% tested before being used in simulations.
    !%Option hscv_pbe 5
    !% PBE version of the HSCV pseudopotentials. Check the
    !% documentation of the option <tt>hscv_lda</tt> for details and warnings.
    !%Option pseudodojo_pbe 100
    !% PBE version of the pseudopotentials of http://pseudo-dojo.org. Version 0.5
    !%Option pseudodojo_lda 103
    !% LDA pseudopotentials of http://pseudo-dojo.org. Version 0.41.
    !%Option pseudodojo_pbesol 105
    !% PBEsol version of the pseudopotentials of http://pseudo-dojo.org. Version 0.41.
    !%Option pseudodojo_pbe_fr 106
    !% Fully-relativistic PBE version of the pseudopotentials of http://pseudo-dojo.org. Version 0.4.
    !%End
    default_val = option__pseudopotentialset__standard
    if(factory%default_allelectron_type /= option__allelectrontype__no) default_val = option__pseudopotentialset__none
    call parse_variable(namespace, 'PseudopotentialSet', default_val, factory%default_pseudopotential_set_id)
    call messages_print_var_option('PseudopotentialSet', factory%default_pseudopotential_set_id, namespace=namespace)
    if (factory%default_pseudopotential_set_id /= option__pseudopotentialset__none) then
      call pseudo_set_init(factory%default_pseudopotential_set, get_set_directory(factory%default_pseudopotential_set_id), ierr)
    else
      call pseudo_set_nullify(factory%default_pseudopotential_set)
    end if
 
    if (factory%default_pseudopotential_set_id /= option__pseudopotentialset__none &
      .and. factory%default_allelectron_type /= option__allelectrontype__no) then
      message(1) = "PseudopotentialSet /= none cannot be used with AllElectronType /= no."
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub_with_profile(species_factory_init)
  end subroutine species_factory_init
 
  ! ---------------------------------------------------------
 
  subroutine species_factory_end(factory)
    class(species_factory_t),   intent(inout) :: factory
 
    push_sub(species_factory_end)
 
    if (factory%initialized) then
      call pseudo_set_end(factory%default_pseudopotential_set)
      factory%initialized = .false.
    end if
 
    pop_sub(species_factory_end)
  end subroutine species_factory_end
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  function species_factory_create_from_input(factory, namespace, label, index) result(spec)
    class(species_factory_t), intent(in) :: factory
    type(namespace_t),        intent(in) :: namespace
    character(len=*),         intent(in) :: label
    integer,                  intent(in) :: index
    class(species_t),         pointer    :: spec
 
    character(len=LABEL_LEN)  :: lab
    integer :: ib, row, n_spec_block, read_data
    type(block_t) :: blk
    type(element_t) :: element
 
    push_sub(species_factory_create_from_input)
 
    read_data   = 0
 
    !%Variable Species
    !%Type block
    !%Section System::Species
    !%Description
    !% A species is by definition either an "ion" (nucleus + core electrons) described
    !% through a pseudopotential, or a model potential.
    !%
    !% Note that some sets of pseudopotentials are distributed with
    !% the code. To use these pseudopotentials, you do not need to define them
    !% explicitly in the <tt>Species</tt> block, as default parameters
    !% are provided.
    !% You can select the set for default pseudopotentials using the
    !% <tt>PseudopotentialSet</tt> variable.
    !%
    !% Supported norm-conserving pseudopotential formats are
    !% detected by the file extension: UPF (<tt>.upf</tt>), PSF (SIESTA, <tt>.psf</tt>), FHI (ABINIT 6, <tt>.fhi</tt>),
    !% CPI (Fritz-Haber, <tt>.cpi</tt>), QSO (quantum-simulation.org, for Qbox, <tt>.xml</tt>),
    !% HGH (Hartwigsen-Goedecker-Hutter, <tt>.hgh</tt>).
    !% PSPIO format can also be used via <tt>species_pspio</tt> if that library is linked.
    !% Note: pseudopotentials may only be used in 3D.
    !%
    !% The format of this block is the following: The first field is a
    !% string that defines the name of the species. The second 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>%Species
    !% <br>&nbsp;&nbsp;'O'       | species_pseudo         | file | 'O.psf' | lmax |  1 | lloc | 1
    !% <br>&nbsp;&nbsp;'H'       | species_pseudo         | file | '../H.hgh'
    !% <br>&nbsp;&nbsp;'Xe'      | species_pseudo         | set | pseudojo_pbe
    !% <br>&nbsp;&nbsp;'C'       | species_pseudo         | file | "carbon.xml"
    !% <br>&nbsp;&nbsp;'jlm'     | species_jellium        | jellium_radius | 5.0
    !% <br>&nbsp;&nbsp;'rho'     | species_charge_density | density_formula | "exp(-r/a)" | mass | 17.0 | valence | 6
    !% <br>&nbsp;&nbsp;'udf'     | species_user_defined   | potential_formula | "1/2*r^2" | valence | 8
    !% <br>&nbsp;&nbsp;'He_all'  | species_full_delta
    !% <br>&nbsp;&nbsp;'H_all'   | species_full_gaussian  |  gaussian_width |  0.2
    !% <br>&nbsp;&nbsp;'Li1D'    | species_soft_coulomb   |  softening | 1.5 | valence | 3
    !% <br>&nbsp;&nbsp;'H_all'   | species_full_anc       |  anc_a | 4
    !% <br>%</tt>
    !%Option species_pseudo  -7
    !% The species is a pseudopotential. How to get the
    !% pseudopotential can be specified by the <tt>file</tt> or
    !% the <tt>set</tt> parameters. If both are missing, the
    !% pseudopotential will be taken from the <tt>PseudopotentialSet</tt>
    !% specified for the run, this is useful if you want to change
    !% some parameters of the pseudo, like the <tt>mass</tt>.
    !%
    !% The optional parameters for this type of species are
    !% <tt>lmax</tt>, that defines the maximum angular momentum
    !% component to be used, and <tt>lloc</tt>, that defines the
    !% angular momentum to be considered as local. When these
    !% parameters are not set, the value for lmax is the maximum
    !% angular component from the pseudopotential file. The default
    !% value for <tt>lloc</tt> is taken from the pseudopotential if
    !% available, if not, it is set to 0. Note that, depending on the
    !% type of pseudopotential, it might not be possible to select
    !% <tt>lmax</tt> and <tt>lloc</tt>, if that is the case the
    !% parameters will be ignored.
    !%
    !%Option species_pspio  -110
    !% (experimental) Alternative method to read pseudopotentials
    !% using the PSPIO library. This species uses the same parameters
    !% as <tt>species_pseudo</tt>.
    !%Option species_user_defined -123
    !% Species with user-defined potential. The potential for the
    !% species is defined by the formula given by the <tt>potential_formula</tt>
    !% parameter.
    !% The
    !% <tt>valence</tt> parameter determines the number of electrons
    !% associated with the species. By default, a valence of 0 is assumed.
    !%Option species_charge_density -125
    !% The potential for this species is created from the distribution
    !% of charge given by the <tt>density_formula</tt> parameter.
    !% The
    !% <tt>valence</tt> parameter determines the number of electrons
    !% associated with the species. By default, a valence of 0 is assumed.
    !%Option species_point  -3
    !%Option species_jellium  -3
    !% Jellium sphere.
    !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
    !%Option species_jellium_slab  -4
    !% A slab of jellium that extends across the simulation box in the
    !% <i>xy</i>-plane. The dimension along the <i>z</i> direction is
    !% determined by the required parameter <tt>thickness</tt>.
    !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
    !%Option species_full_delta   -127
    !% Full atomic potential represented by a delta charge
    !% distribution. The atom will be displaced to the nearest grid
    !% point. The atomic number is determined from the name of the species.
    !%Option species_full_anc     -130
    !% Analytical norm-conserving regularized Coulomb potential from
    !% [Gygi J. Chem. Theory Comput. 2023, 19, 1300−1309].
    !% Unlike <tt>species_full_delta</tt> and <tt>species_full_gaussian</tt>, this species has
    !% no nuclear charge density: the potential is applied directly, without solving a Poisson
    !% equation for the nucleus. The <tt>External</tt> and <tt>Hartree</tt> components reported
    !% in <tt>static/info</tt> are therefore split differently than for the other two types.
    !%Option species_full_gaussian   -124
    !% A full-potential atom is defined by a Gaussian accumulation of
    !% positive charge, in the form:
    !%
    !% <math>\rho(r) = Z (\zeta/\pi)^{3/2} e^{-\zeta (\vec{r}-\vec{r_0})^2}</math>
    !%
    !% with <math>\zeta = 1/(2\sigma^2\delta^2)</math> in bohr<sup>-2</sup>, where
    !% <math>\sigma</math> is set by the parameter <tt>gaussian_width</tt> and
    !% <math>\delta</math> is the grid spacing (the grid spacing in the first dimension,
    !% to be precise). The prefactor maintains proper normalization: the integral of
    !% <math>\rho</math> sums up to <math>Z</math>.
    !% <math>\vec{r_0}</math> is calculated in such a way that the
    !% first moment of <math>\rho(r)/Z</math> is equal to the atomic
    !% position. For a precise description, see N. A. Modine,
    !% <i>Phys. Rev. B</i> <b>55</b>, 10289 (1997).
    !%
    !% The Gaussian is built in the primitive coordinates of the mesh, whose spacing is
    !% uniform, so with curvilinear coordinates it is narrower in real space where the grid
    !% is denser. This works with <tt>curv_gygi</tt> and <tt>curv_modine</tt>, but not with
    !% <tt>curv_briggs</tt>, and not for periodic systems.
    !% The atomic number is determined from the name of the species.
    !%Option species_from_file  -126
    !% The potential is read from a file. Accepted file formats, detected by extension: obf, ncdf and csv.
    !% The
    !% <tt>valence</tt> parameter determines the number of electrons
    !% associated with the species. By default, a valence of 0 is assumed.
    !%Option species_soft_coulomb -128
    !% The potential is a soft-Coulomb function, <i>i.e.</i> a function in the form:
    !%
    !% <math>v(r) = - z_{val} / \sqrt{a^2 + r^2}</math>
    !%
    !% The value of <i>a</i> should be given by the mandatory <tt>softening</tt> parameter.
    !% The charge associated with this species must be given by the <tt>valence</tt> parameter.
    !%Option species_jellium_charge_density -129
    !% The parameter is the name of a volume block specifying the shape of the jellium.
    !%Option lmax -10003
    !% The maximum angular-momentum channel that will be used for the pseudopotential.
    !%Option lloc -10004
    !% The angular-momentum channel of the pseudopotential to be considered local.
    !%Option mass -10005
    !% The mass of the species in atomic mass units, <i>i.e.</i> the mass of a proton is
    !% roughly one. It is set automatically for pseudopotentials from the
    !% <a href=http://www.nist.gov/pml/data/comp.cfm>NIST values</a>.
    !% For other species, the default is 1.0.
    !%Option valence -10006
    !% The number of electrons of the species. It is set automatically from the name of the species.
    !% if it correspond to the name in the periodic table. If not specified and if the name
    !% does not match an atom name, a value of 0 is assumed.
    !%Option jellium_radius -10007
    !% The radius of the sphere for <tt>species_jellium</tt>. If this value is not specified,
    !% the default of 0.5 bohr is used.
    !%Option set -10017
    !% For a <tt>species_pseudo</tt>, get the pseudopotential from a
    !% particular set. This flag must be followed with one of the
    !% valid values for the variable <tt>PseudopotentialSet</tt>.
    !%Option gaussian_width -10008
    !% The width <math>\sigma</math> of the Gaussian (in units of the grid spacing) used to
    !% represent the nuclear charge for <tt>species_full_gaussian</tt>. If not present, the
    !% value of <tt>AllElectronSigma</tt> is used, which itself defaults to 0.6. See
    !% <tt>AllElectronSigma</tt> for the consequences of tying the nuclear size to the grid.
    !%Option softening -10009
    !% The softening parameter <i>a</i> for <tt>species_soft_coulomb</tt> in units of length.
    !%Option file -10010
    !% The path for the file that describes the species.
    !%Option db_file -10011
    !% Obsolete. Use the <tt>set</tt> option of the <tt>PseudopotentialSet</tt> variable instead.
    !%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.
    !%Option thickness -10014
    !% The thickness of the slab for species_jellium_slab. Must be positive.
    !%Option vdw_radius -10015
    !% The van der Waals radius that will be used for this species.
    !%Option volume -10016
    !% Name of a volume block
    !%Option hubbard_l -10018
    !% The angular-momentum for which the effective U will be applied.
    !%Option hubbard_u -10019
    !% The effective U that will be used for the DFT+U calculations.
    !%Option hubbard_j -10020
    !% The value of j (hubbard_l-1/2 or hubbard_l+1/2) on which the effective U is applied.
    !%Option hubbard_alpha -10021
    !% The strength of the potential constraining the occupations of the localized subspace
    !% as defined in PRB 71, 035105 (2005)
    !%Option anc_a -10022
    !% The value of the parameter a of the ANC potential, as defined in [Gygi, JCTC 2023, 19, 1300−1309].
    !% This parameter has the unit of inverse length and determines the range of regularization.
    !% The code scales the radius as <math>a Z r</math>, so the regularization length is
    !% <math>1/(aZ)</math> bohr and shrinks with the atomic number. If not present, the value of
    !% <tt>AllElectronANCParam</tt> is used, which defaults to 4, a value validated for hydrogen
    !% only.
    !%End
 
    call messages_obsolete_variable(namespace, 'SpecieAllElectronSigma', 'Species')
    call messages_obsolete_variable(namespace, 'SpeciesAllElectronSigma', 'Species')
 
    ! First, find out if there is a Species block.
    n_spec_block = 0
    if (parse_block(namespace, 'Species', blk) == 0) then
      n_spec_block = parse_block_n(blk)
    end if
 
    ! Find out if the sought species is in the block
    row = -1
    block: do ib = 1, n_spec_block
      call parse_block_string(blk, ib-1, 0, lab)
      if (trim(lab) == trim(label)) then
        row = ib - 1
        exit block
      end if
    end do block
 
    ! Read whatever may be read from the block
    if (row >= 0) then
      spec => factory%create_from_block(namespace, blk, row, label, index, read_data)
      call parse_block_end(blk)
 
      assert(read_data > 0)
 
      pop_sub(species_factory_create_from_input)
      return
    end if
 
    ! We get here if there is a Species block but it does not contain
    ! the species we are looking for.
    if (n_spec_block > 0) call parse_block_end(blk)
 
    ! Initialize all electron species (except soft-Coulomb) from specified allelectron type
    if(factory%default_allelectron_type /= option__allelectrontype__no) then
      select case(factory%default_allelectron_type)
      case(option__allelectrontype__full_delta)
        spec => full_delta_t(label, index, factory%default_sigma)
      case(option__allelectrontype__full_gaussian)
        spec => full_gaussian_t(label, index, factory%default_sigma)
      case(option__allelectrontype__full_anc)
        spec => full_anc_t(label, index, factory%default_anc_a)
      case default
        assert(.false.)
      end select
 
      ! get the mass, vdw radius and atomic number for this element
      call element_init(element, get_symbol(spec%get_label()))
 
      assert(element_valid(element))
 
      call spec%set_z(real(element_atomic_number(element), real64))
      call spec%set_zval(spec%get_z())
      call spec%set_mass(element_mass(element))
      call spec%set_vdw_radius(element_vdw_radius(element))
 
      call element_end(element)
 
    else ! Pseudopotential from specified set or the default one
      spec => pseudopotential_t(label, index)
      select type(spec)
      type is(pseudopotential_t)
        call read_from_set(spec, factory%default_pseudopotential_set_id, factory%default_pseudopotential_set, read_data)
 
        if (read_data == 0) then
          call messages_write( 'Species '//trim(spec%get_label())//' not found in default pseudopotential set.', new_line=.true. )
          call messages_write('( '//trim(get_set_directory(factory%default_pseudopotential_set_id))//' )')
          call messages_fatal(namespace=namespace)
        end if
      end select
    end if
 
    pop_sub(species_factory_create_from_input)
  end function species_factory_create_from_input
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Private procedures
 
  ! ---------------------------------------------------------
  function read_from_block(factory, namespace, blk, row, label, index, read_data) result(spec)
    class(species_factory_t), intent(in) :: factory
    type(namespace_t),        intent(in) :: namespace
    type(block_t),            intent(in) :: blk
    integer,                  intent(in) :: row
    character(len=*),         intent(in) :: label
    integer,                  intent(in) :: index
    integer,                  intent(out):: read_data
    class(species_t),         pointer    :: spec
 
    integer :: ncols, icol, flag, set_read_data, ierr, type
    type(element_t) :: element
    type(iihash_t) :: read_parameters
    integer :: user_lmax, user_llocal, hubbard_l, pseudopotential_set_id
    real(real64) :: hubbard_u, hubbard_j, hubbard_alpha, mass, z_val, jradius, jthick, vdw_radius, aa
    real(real64) :: sigma, softening
    character(len=MAX_PATH_LEN) :: filename
 
    integer, parameter ::  &
      species_jellium        = 3,             & !< jellium sphere.
      species_jellium_slab   = 4,             & 
      species_jellium_charge_density = 129,   & 
      species_pseudo         = 7,             & 
      species_pspio          = 110,           & 
      species_usdef          = 123,           & 
      species_full_gaussian  = 124,           & 
      species_charge_density = 125,           & 
      species_from_file      = 126,           &
      species_full_delta     = 127,           & 
      species_soft_coulomb   = 128,           & 
      species_full_anc       = 130              
 
    push_sub(read_from_block)
 
    ncols = parse_block_cols(blk, row)
    read_data = 0
 
    call parse_block_integer(blk, row, 1, 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 (type >= 0) then
      call messages_write('Found  a species  with the old format.  Please update', new_line = .true.)
      call messages_write('the Species block to the new format, where the second', new_line = .true.)
      call messages_write('column indicates the type of the species.')
      call messages_fatal(namespace=namespace)
    end if
 
    ! now we convert back to positive
    type = -type
 
    read_data = 2
 
 
    select case (type)
 
    case (species_soft_coulomb)
      spec => soft_coulomb_t(label, index)
 
    case (species_usdef) ! user-defined
      spec => species_user_defined_t(label, index)
 
    case (species_from_file)
      spec => species_from_file_t(label, index)
 
    case (species_jellium)
      spec => jellium_sphere_t(label, index)
 
    case (species_jellium_slab)
      spec => jellium_slab_t(label, index)
 
    case (species_full_delta)
      spec => full_delta_t(label, index, factory%default_sigma)
 
    case (species_full_gaussian)
      spec => full_gaussian_t(label, index, factory%default_sigma)
 
    case (species_full_anc)
      spec => full_anc_t(label, index, factory%default_anc_a)
 
    case (species_charge_density)
      spec => species_charge_density_t(label, index)
 
    case (species_jellium_charge_density)
      spec => jellium_charge_t(label, index)
 
    case (species_pseudo)
      spec => pseudopotential_t(label, index)
 
    case (species_pspio) ! a pseudopotential file to be handled by the pspio library
      spec => pseudopotential_t(label, index)
 
    case default
      call messages_input_error(namespace, 'Species', "Unknown type for species '"//trim(spec%get_label())//"'", row=row, column=1)
    end select
 
    call spec%set_mass(-m_one)
    call spec%set_vdw_radius(-m_one)
    call spec%set_zval(-m_one)
 
    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__species__lmax)
        call check_duplication(option__species__lmax)
        call parse_block_integer(blk, row, icol + 1, user_lmax)
 
        select type(spec)
        class is(pseudopotential_t)
          call spec%set_user_lmax(user_lmax)
        class default
          call messages_input_error(namespace, 'Species', &
            "The 'lmax' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
            row=row, column=icol+1)
        end select
 
        if (user_lmax < 0) then
          call messages_input_error(namespace, 'Species', &
            "The 'lmax' parameter in species "//trim(spec%get_label())//" cannot be negative", &
            row=row, column=icol+1)
        end if
 
      case (option__species__lloc)
        call check_duplication(option__species__lloc)
        call parse_block_integer(blk, row, icol + 1, user_llocal)
 
        select type(spec)
        class is(pseudopotential_t)
          call spec%set_user_lloc(user_llocal)
        class default
          call messages_input_error(namespace, 'Species', &
            "The 'lloc' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
            row=row, column=icol+1)
        end select
 
        if (user_llocal < 0) then
          call messages_input_error(namespace, 'Species', &
            "The 'lloc' parameter in species "//trim(spec%get_label())//" cannot be negative", row=row, column=icol+1)
        end if
 
      case (option__species__hubbard_l)
        call check_duplication(option__species__hubbard_l)
        call parse_block_integer(blk, row, icol + 1, hubbard_l)
 
        select type(spec)
        class is(pseudopotential_t)
          call spec%set_hubbard_l(hubbard_l)
        class default
          call messages_input_error(namespace, 'Species', &
            "The 'hubbard_l' parameter in species "//trim(spec%get_label())//" can only be used with pseudopotential species", &
            row=row, column=icol+1)
        end select
 
        if (hubbard_l < 0) then
          call messages_input_error(namespace, 'Species', &
            "The 'hubbard_l' parameter in species "//trim(spec%get_label())//" cannot be negative", row=row, column=icol+1)
        end if
 
      case (option__species__hubbard_u)
        call check_duplication(option__species__hubbard_u)
        call parse_block_float(blk, row, icol + 1, hubbard_u, unit = units_inp%energy)
        call spec%set_hubbard_u(hubbard_u)
 
      case (option__species__hubbard_alpha)
        call check_duplication(option__species__hubbard_alpha)
        call parse_block_float(blk, row, icol + 1, hubbard_alpha, unit = units_inp%energy)
        call spec%set_hubbard_alpha(hubbard_alpha)
 
      case (option__species__hubbard_j)
        call check_duplication(option__species__hubbard_j)
        call parse_block_float(blk, row, icol + 1, hubbard_j)
        call spec%set_hubbard_j(hubbard_j)
 
        if (abs(abs(spec%get_hubbard_j()-spec%get_hubbard_l())-m_half) <= m_epsilon) then
          call messages_input_error(namespace, 'Species', "The 'hubbard_j' parameter in species "// &
            trim(spec%get_label())//" can only be hubbard_l +/- 1/2", row=row, column=icol+1)
        end if
 
      case (option__species__mass)
        call check_duplication(option__species__mass)
        call parse_block_float(blk, row, icol + 1, mass, unit = units_inp%mass)
        call spec%set_mass(mass)
 
      case (option__species__valence)
        call check_duplication(option__species__valence)
        call parse_block_float(blk, row, icol + 1, z_val)
        call spec%set_zval(z_val)
        call spec%set_z(z_val)
 
      case (option__species__jellium_radius)
        call check_duplication(option__species__jellium_radius)
        call parse_block_float(blk, row, icol + 1, jradius)
        if (jradius <= m_zero) call messages_input_error(namespace, 'Species', 'jellium_radius must be positive', &
          row=row, column=icol+1)
 
        select type(spec)
        type is(jellium_sphere_t)
          call spec%set_radius(jradius)
        class default
          call messages_input_error(namespace, 'Species', 'jellium_radius can only be used with species_jellium', &
            row=row, column=icol+1)
        end select
 
      case (option__species__gaussian_width)
        call check_duplication(option__species__gaussian_width)
        call parse_block_float(blk, row, icol + 1, sigma)
        if (sigma <= m_zero) call messages_input_error(namespace, 'Species', 'gaussian_width must be positive', &
          row=row, column=icol+1)
        select type(spec)
        type is(full_gaussian_t)
          call spec%set_sigma(sigma)
        class default
          call messages_input_error(namespace, 'Species', 'gaussian_width can only be used with species_full_gaussian', &
            row=row, column=icol+1)
        end select
 
      case (option__species__anc_a)
        call check_duplication(option__species__anc_a)
        call parse_block_float(blk, row, icol + 1, aa)
        if (aa <= m_zero) call messages_input_error(namespace, 'Species', 'anc_a must be positive', &
          row=row, column=icol+1)
 
        select type(spec)
        type is(full_anc_t)
          call spec%set_a(aa)
        class default
          call messages_input_error(namespace, 'Species', 'anc_a can only be used with species_full_anc', &
            row=row, column=icol+1)
        end select
 
      case (option__species__softening)
        call check_duplication(option__species__softening)
        call parse_block_float(blk, row, icol + 1, softening)
        softening = softening**2
 
        select type(spec)
        type is(soft_coulomb_t)
          call spec%set_softening2(softening)
        class default
          call messages_input_error(namespace, 'Species', 'softening can only be used with species_soft_coulomb', &
            row=row, column=icol+1)
        end select
 
      case (option__species__file)
        call check_duplication(option__species__file)
        call parse_block_string(blk, row, icol + 1, filename)
        call spec%set_filename(filename)
 
      case (option__species__db_file)
        call messages_write("The 'db_file' option for 'Species' block is obsolete. Please use", new_line = .true.)
        call messages_write("the option 'set' or the variable 'PseudopotentialSet' instead.")
        call messages_fatal(namespace=namespace)
 
      case (option__species__set)
        call check_duplication(option__species__set)
        call parse_block_integer(blk, row, icol + 1, pseudopotential_set_id)
 
        select type(spec)
        type is(pseudopotential_t)
          spec%pseudopotential_set_initialized = .true.
          spec%pseudopotential_set_id = pseudopotential_set_id
          call pseudo_set_init(spec%pseudopotential_set, get_set_directory(spec%pseudopotential_set_id), ierr)
        class default
          call messages_input_error(namespace, 'Species', 'set can only be used with species_pseudo', &
            row=row, column=icol+1)
        end select
 
      case (option__species__potential_formula)
        call check_duplication(option__species__potential_formula)
        select type(spec)
        type is(species_user_defined_t)
          call parse_block_string(blk, row, icol + 1, spec%potential_formula, convert_to_c=.true.)
        class default
          call messages_input_error(namespace, 'Species', 'potential_formula can only be used with species_user_defined', &
            row=row, column=icol+1)
        end select
 
      case (option__species__volume)
        call check_duplication(option__species__volume)
 
        select type(spec)
        type is(jellium_charge_t)
          call parse_block_string(blk, row, icol + 1, spec%density_formula, convert_to_c=.true.)
 
        class default
          call messages_input_error(namespace, 'Species', 'volume can only be used with species_jellium_charge_density', &
            row=row, column=icol+1)
        end select
 
      case (option__species__density_formula)
        call check_duplication(option__species__density_formula)
 
        select type(spec)
        type is(species_charge_density_t)
          call parse_block_string(blk, row, icol + 1, spec%density_formula, convert_to_c=.true.)
 
        class default
          call messages_input_error(namespace, 'Species', 'density_formula can only be used with species_charge_density', &
            row=row, column=icol+1)
        end select
 
      case (option__species__thickness)
        call check_duplication(option__species__thickness)
        call parse_block_float(blk, row, icol + 1, jthick) ! thickness of the jellium slab
 
        if (jthick <= m_zero) then
          call messages_input_error(namespace, 'Species', 'the value of the thickness parameter in species '&
            //trim(spec%get_label())//' must be positive.', row=row, column=icol+1)
        end if
 
        select type(spec)
        type is(jellium_slab_t)
          call spec%set_thickness(jthick)
        class default
          call messages_input_error(namespace, 'Species', 'thickness can only be used with species_jellium_slab', &
            row=row, column=icol+1)
        end select
 
      case (option__species__vdw_radius)
        call check_duplication(option__species__vdw_radius)
        call parse_block_float(blk, row, icol + 1, vdw_radius, unit = units_inp%length)
        call spec%set_vdw_radius(vdw_radius)
 
      case default
        call messages_input_error(namespace, 'Species', "Unknown parameter in species '"//trim(spec%get_label())//"'", &
          row=row, column=icol)
 
      end select
 
      icol = icol + 2
    end do
    ! CHECK THAT WHAT WE PARSED MAKES SENSE
 
    select type(spec)
    type is(soft_coulomb_t)
      if (.not. parameter_defined(option__species__softening)) then
        call messages_input_error(namespace, 'Species', &
          "The 'softening' parameter is missing for species "//trim(spec%get_label()))
      end if
 
    type is(species_user_defined_t)
      if (.not. parameter_defined(option__species__potential_formula)) then
        call messages_input_error(namespace, 'Species', &
          "The 'potential_formula' parameter is missing for species '"//trim(spec%get_label())//"'")
      end if
 
    type is(species_charge_density_t)
      if (.not. parameter_defined(option__species__density_formula)) then
        call messages_input_error(namespace, 'Species', &
          "The 'density_formula' parameter is missing for species '"//trim(spec%get_label())//"'")
      end if
 
    type is(species_from_file_t)
      if( .not. (parameter_defined(option__species__file) .or. parameter_defined(option__species__db_file))) then
        call messages_input_error(namespace, 'Species', &
          "The 'file' or 'db_file' parameter is missing for species '"//trim(spec%get_label())//"'")
      end if
 
    type is(jellium_slab_t)
      if (.not. parameter_defined(option__species__thickness)) then
        call messages_input_error(namespace, 'Species', &
          "The 'thickness' parameter is missing for species '"//trim(spec%get_label())//"'")
      end if
 
    type is(jellium_charge_t)
      if (.not. parameter_defined(option__species__volume)) then
        call messages_input_error(namespace, 'Species', &
          "The 'volume' parameter is missing for species '"//trim(spec%get_label())//"'")
      end if
 
    type is(pseudopotential_t)
      if (parameter_defined(option__species__lmax) .and. parameter_defined(option__species__lloc)) then
        if (spec%get_user_lloc() > spec%get_user_lmax()) then
          call messages_input_error(namespace, 'Species', &
            "the 'lloc' parameter cannot be larger than the 'lmax' parameter in species "//trim(spec%get_label()))
        end if
      end if
 
      if(.not. (parameter_defined(option__species__file) .or. parameter_defined(option__species__db_file))) then
        ! we need to read the species from the pseudopotential set
 
        !if the set was not defined, use the default set
        if (.not. parameter_defined(option__species__set)) then
          spec%pseudopotential_set_id = factory%default_pseudopotential_set_id
          spec%pseudopotential_set = factory%default_pseudopotential_set
        end if
 
        call read_from_set(spec, spec%pseudopotential_set_id, spec%pseudopotential_set, set_read_data)
 
        if (set_read_data == 0) then
          call messages_write('Species '//trim(spec%get_label())//' is not defined in the requested pseudopotential set.', &
            new_line=.true.)
          call messages_write('( '//trim(get_set_directory(spec%pseudopotential_set_id))//' )')
          call messages_fatal(namespace=namespace)
        end if
      end if
 
    end select
 
    select type (spec)
    class is(pseudopotential_t)
      call check_real_atom_species()
 
    class is(full_anc_t)
      call check_real_atom_species()
 
      ! If z_val was not specified, we set it to be z
      if (spec%get_zval() < m_zero) then
        call spec%set_zval(spec%get_z())
      end if
 
    class is(full_gaussian_t)
      call check_real_atom_species()
 
      ! If z_val was not specified, we set it to be z
      if (spec%get_zval() < m_zero) then
        call spec%set_zval(spec%get_z())
      end if
 
    class is(full_delta_t)
      call check_real_atom_species()
 
      ! If z_val was not specified, we set it to be z
      if (spec%get_zval() < m_zero) then
        call spec%set_zval(spec%get_z())
      end if
 
    class default
      if (.not. parameter_defined(option__species__mass)) then
        call spec%set_mass(m_one)
        call messages_write('Info: default mass for species '//trim(spec%get_label())//':')
        call messages_write(spec%get_mass())
        call messages_write(' amu.')
        call messages_info(namespace=namespace)
      end if
 
      if (.not. parameter_defined(option__species__vdw_radius)) then
        call spec%set_vdw_radius(m_zero)
        call messages_write('Info: default vdW radius for species '//trim(spec%get_label())//':')
        call messages_write(spec%get_vdw_radius())
        call messages_write(' [b]')
        call messages_info(namespace=namespace)
      end if
 
      if (.not. parameter_defined(option__species__valence)) then
        if (spec%is_user_defined()) then
          call spec%set_zval(m_zero)
        else
          call messages_input_error(namespace, 'Species', &
            "The 'valence' parameter is missing for species '"//trim(spec%get_label())//"'")
        end if
      end if
 
    end select
 
    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, 'Species', "Duplicated parameter in species '"//trim(spec%get_label())//"'")
      end if
 
      call iihash_insert(read_parameters, int(-param), 1)
 
      pop_sub(read_from_block.check_duplication)
    end subroutine check_duplication
 
 
    !------------------------------------------------------
    subroutine check_real_atom_species()
 
      call element_init(element, get_symbol(spec%get_label()))
 
      if (.not. element_valid(element)) then
        call messages_write('Cannot determine the element for species '//trim(spec%get_label())//'.')
        call messages_fatal(namespace=namespace)
      end if
 
      call spec%set_z(real(element_atomic_number(element), real64))
 
      if (spec%get_mass() < m_zero) then
        call spec%set_mass(element_mass(element))
        call messages_write('Info: default mass for species '//trim(spec%get_label())//':')
        call messages_write(spec%get_mass())
        call messages_write(' amu.')
        call messages_info(namespace=namespace)
      end if
 
      if (spec%get_vdw_radius() < m_zero) then
        call spec%set_vdw_radius(element_vdw_radius(element))
        if (spec%get_vdw_radius() < m_zero) then
          call spec%set_vdw_radius(m_zero)
          call messages_write("The default vdW radius for species '"//trim(spec%get_label())//"' is not defined.", &
            new_line = .true.)
          call messages_write("You can specify the vdW radius in %Species block.")
          call messages_warning(namespace=namespace)
        end if
        call messages_write('Info: default vdW radius for species '//trim(spec%get_label())//':')
        call messages_write(spec%get_vdw_radius())
        call messages_write(' [b]')
        call messages_info(namespace=namespace)
      end if
 
      call element_end(element)
 
 
    end subroutine check_real_atom_species
  end function read_from_block
 
end module species_factory_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: