epot.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 epot_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use ion_interaction_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use projector_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_oct_m
  use species_pot_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use spline_filter_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use tdfunction_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                        &
    epot_t,                        &
    epot_init,                     &
    epot_end,                      &
    epot_bind_poisson_solver,      &
    epot_generate,                 &
    epot_local_potential,          &
    epot_precalc_local_potential,  &
    epot_have_external_potentials, &
    local_potential_has_density
 
  integer, public, parameter ::   &
    NOREL      = 0,               &
    spin_orbit = 1,               &
    scalar_relativistic_zora = 2, &
    fully_relativistic_zora  = 3
 
  type epot_t
    ! Components are public by default
 
    ! Ions
    real(real64),      allocatable :: vpsl(:)
    !                                               !< plus the potential from static electric fields
    type(projector_t), allocatable :: proj(:)
    logical                        :: non_local
    integer                        :: natoms
 
    ! External e-m fields
    real(real64), allocatable     :: e_field(:)
    real(real64), allocatable     :: v_ext(:)
    real(real64), allocatable     :: b_field(:)
    real(real64), allocatable     :: a_static(:,:)
    integer                :: reltype
    integer                :: proj_reltype
 
    real(real64) :: gyromagnetic_ratio
 
    real(real64) :: so_strength
 
    real(real64)       :: eii
    real(real64), allocatable :: fii(:, :)
    real(real64), allocatable :: vdw_forces(:, :)
    real(real64), allocatable :: photon_forces(:)
 
    real(real64) :: vdw_stress(3, 3)
 
    real(real64), allocatable, private :: local_potential(:,:)
    logical,            private :: local_potential_precalculated
 
    logical,                  private :: have_density
    type(poisson_t), pointer, private :: poisson_solver
 
    logical, private :: proj_check_done = .false.
 
    logical              :: nlcc = .false.   
  end type epot_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine epot_init(ep, namespace, gr, ions, psolver, ispin, xc_family, kpoints)
    type(epot_t),                       intent(out)   :: ep
    type(namespace_t),                  intent(in)    :: namespace
    type(grid_t),                       intent(in)    :: gr
    type(ions_t),                       intent(inout) :: ions
    type(poisson_t),  target,           intent(in)    :: psolver
    integer,                            intent(in)    :: ispin
    integer,                            intent(in)    :: xc_family
    type(kpoints_t),                    intent(in)    :: kpoints
 
 
    integer :: ispec, ia
    integer :: filter
 
    push_sub(epot_init)
 
    !%Variable FilterPotentials
    !%Type integer
    !%Default filter_ts
    !%Section Hamiltonian
    !%Description
    !% <tt>Octopus</tt> can filter the pseudopotentials so that they no
    !% longer contain Fourier components larger than the mesh itself. This is
    !% very useful to decrease the egg-box effect, and so should be used in
    !% all instances where atoms move (<i>e.g.</i> geometry optimization,
    !% molecular dynamics, and vibrational modes).
    !% No filtering is applied to HGH pseudopotentials, as these are already
    !% smooth, analytical potentials in both real and Fourier space.
    !%Option filter_none 0
    !% Do not filter.
    !%Option filter_TS 2
    !% The filter of M. Tafipolsky and R. Schmid, <i>J. Chem. Phys.</i> <b>124</b>, 174102 (2006).
    !%Option filter_BSB 3
    !% The filter of E. L. Briggs, D. J. Sullivan, and J. Bernholc, <i>Phys. Rev. B</i> <b>54</b>, 14362 (1996).
    !%End
    call parse_variable(namespace, 'FilterPotentials', ps_filter_ts, filter)
    if (.not. varinfo_valid_option('FilterPotentials', filter)) call messages_input_error(namespace, 'FilterPotentials')
    call messages_print_var_option("FilterPotentials", filter, namespace=namespace)
 
    if (family_is_mgga(xc_family) .and. filter /= ps_filter_none) then
      call messages_not_implemented("FilterPotentials different from filter_none with MGGA", namespace=namespace)
    end if
 
    if (filter == ps_filter_ts) call spline_filter_mask_init()
    do ispec = 1, ions%nspecies
      call ions%species(ispec)%s%init_potential(namespace, mesh_gcutoff(gr), filter)
    end do
 
    safe_allocate(ep%vpsl(1:gr%np))
 
    ep%vpsl(1:gr%np) = m_zero
 
    ! No more "UserDefinedTDPotential" from this version on.
    call messages_obsolete_variable(namespace, 'UserDefinedTDPotential', 'TDExternalFields')
 
    call messages_obsolete_variable(namespace, 'ClassicalPotential')
 
    !%Variable GyromagneticRatio
    !%Type float
    !%Default 2.0023193043768
    !%Section Hamiltonian
    !%Description
    !% The gyromagnetic ratio of the electron. This is of course a physical
    !% constant, and the default value is the exact one that you should not
    !% touch, unless:
    !% (i)  You want to disconnect the anomalous Zeeman term in the Hamiltonian
    !% (then set it to zero; this number only affects that term);
    !% (ii) You are using an effective Hamiltonian, as is the case when
    !% you calculate a 2D electron gas, in which case you have an effective
    !% gyromagnetic factor that depends on the material.
    !%End
    call parse_variable(namespace, 'GyromagneticRatio', p_g, ep%gyromagnetic_ratio)
 
    !%Variable RelativisticCorrection
    !%Type integer
    !%Default non_relativistic
    !%Section Hamiltonian
    !%Description
    !% The default value means that <i>no</i> relativistic correction is used. To
    !% include spin-orbit coupling turn <tt>RelativisticCorrection</tt> to <tt>spin_orbit</tt>
    !% (this will only work if <tt>SpinComponents</tt> has been set to <tt>non_collinear</tt>, which ensures
    !% the use of spinors).
    !%
    !% The two ZORA options are meant for all-electron species (see <tt>AllElectronType</tt>).
    !% For pseudopotentials the relativistic effects are already contained in the pseudopotential
    !% itself, so <tt>fully_relativistic_zora</tt> is rejected for pseudopotential species and
    !% <tt>scalar_relativistic_zora</tt> only issues a warning; use <tt>spin_orbit</tt> instead.
    !%
    !% Three properties of the ZORA implementation are worth knowing:
    !%
    !% <ul>
    !% <li> The potential <i>V</i> entering the ZORA prefactor <math>c^2/(2mc^2-V)</math> is the
    !% ionic plus external potential only; the Hartree and exchange-correlation potential are
    !% <i>not</i> included. This is a deliberate approximation, which makes this ZORA comparable
    !% to a one-electron X2C rather than to a full ZORA. </li>
    !% <li> As a consequence the ZORA prefactor carries no spin dependence: it is built from a
    !% single, spin-independent potential and applied identically to both spin channels. </li>
    !% <li> The forces receive no contribution from the ZORA terms, and are therefore only
    !% approximate; a warning is printed if they are requested. The stress tensor is not
    !% available with ZORA at all. </li>
    !% </ul>
    !%Option non_relativistic 0
    !% No relativistic corrections.
    !%Option spin_orbit 1
    !% Spin-orbit, from the j-dependent projectors of a fully-relativistic pseudopotential.
    !%Option scalar_relativistic_zora 2
    !% Scalar relativistic ZORA Hamiltonian
    !%Option fully_relativistic_zora 3
    !% Fully relativistic spin-orbit ZORA Hamiltonian including SR and SO terms
    !%End
    call parse_variable(namespace, 'RelativisticCorrection', norel, ep%reltype)
    if (.not. varinfo_valid_option('RelativisticCorrection', ep%reltype)) then
      call messages_input_error(namespace, 'RelativisticCorrection')
    end if
    if (ispin /= spinors .and. ( ep%reltype == spin_orbit .or. ep%reltype == fully_relativistic_zora ) ) then
      message(1) = "The spin-orbit term can only be applied when using spinors."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if((ep%reltype == spin_orbit .or. ep%reltype == fully_relativistic_zora) .and. kpoints%use_symmetries) then
      call messages_not_implemented("Spin-orbit coupling and k-point symmetries", namespace=namespace)
    end if
 
    if (ep%reltype == scalar_relativistic_zora .or. ep%reltype == fully_relativistic_zora) then
      do ia = 1, ions%nspecies
        if (.not. ions%species(ia)%s%is_ps()) cycle
        message(1) = "ZORA is meant for all-electron species, but species '" &
          //trim(ions%species(ia)%s%get_label())//"' is a pseudopotential,"
        message(2) = "which already contains the relativistic effects."
        if (ep%reltype == fully_relativistic_zora) then
          message(3) = "Use RelativisticCorrection = spin_orbit for pseudopotentials, or an"
          message(4) = "all-electron species (see AllElectronType) for ZORA."
          call messages_fatal(4, namespace=namespace)
        else
          message(3) = "Scalar relativistic effects will therefore be counted twice."
          call messages_warning(3, namespace=namespace)
        end if
      end do
    end if
 
    ep%proj_reltype = norel
    if (ep%reltype == spin_orbit) ep%proj_reltype = spin_orbit
 
    call messages_print_var_option("RelativisticCorrection", ep%reltype, namespace=namespace)
 
    !%Variable SOStrength
    !%Type float
    !%Default 1.0
    !%Section Hamiltonian
    !%Description
    !% Tuning of the spin-orbit coupling strength: setting this value to zero turns off spin-orbit terms in
    !% the Hamiltonian, and setting it to one corresponds to full spin-orbit.
    !%
    !% This variable is parsed both for <tt>RelativisticCorrection = spin_orbit</tt> and for
    !% <tt>fully_relativistic_zora</tt>. In the latter case it scales the ZORA spin-orbit term only,
    !% so <tt>SOStrength = 0.0</tt> reproduces <tt>scalar_relativistic_zora</tt> exactly. This is
    !% the cleanest way to separate the scalar-relativistic and the spin-orbit ZORA contributions
    !% within a single code path.
    !%End
    if (ep%reltype == spin_orbit .or. ep%reltype == fully_relativistic_zora) then
      call parse_variable(namespace, 'SOStrength', m_one, ep%so_strength)
    else
      ep%so_strength = m_one
    end if
 
    safe_allocate(ep%proj(1:ions%natoms))
 
    ep%natoms = ions%natoms
    ep%non_local = .false.
 
    ep%eii = m_zero
    safe_allocate(ep%fii(1:ions%space%dim, 1:ions%natoms))
    ep%fii = m_zero
 
    safe_allocate(ep%vdw_forces(1:ions%space%dim, 1:ions%natoms))
    ep%vdw_forces = m_zero
 
    safe_allocate(ep%photon_forces(1:ions%space%dim))
    ep%photon_forces = m_zero
 
    ep%local_potential_precalculated = .false.
 
 
    ep%have_density = .false.
    do ia = 1, ions%nspecies
      if (local_potential_has_density(ions%space, ions%species(ia)%s)) then
        ep%have_density = .true.
        exit
      end if
    end do
 
    if (ep%have_density) then
      ep%poisson_solver => psolver
    else
      nullify(ep%poisson_solver)
    end if
 
    ! find out if we need non-local core corrections
    ep%nlcc = .false.
    do ia = 1, ions%nspecies
      ep%nlcc = (ep%nlcc .or. ions%species(ia)%s%is_ps_with_nlcc())
    end do
 
    pop_sub(epot_init)
  end subroutine epot_init
 
  ! ---------------------------------------------------------
  subroutine epot_end(ep)
    type(epot_t), intent(inout) :: ep
 
    integer :: iproj
 
    push_sub(epot_end)
 
    if (ep%have_density) then
      nullify(ep%poisson_solver)
    end if
 
    safe_deallocate_a(ep%local_potential)
    safe_deallocate_a(ep%fii)
    safe_deallocate_a(ep%vdw_forces)
    safe_deallocate_a(ep%vpsl)
    safe_deallocate_a(ep%photon_forces)
 
    ! the macroscopic fields
    safe_deallocate_a(ep%e_field)
    safe_deallocate_a(ep%v_ext)
    safe_deallocate_a(ep%b_field)
    safe_deallocate_a(ep%a_static)
 
    do iproj = 1, ep%natoms
      if (projector_is_null(ep%proj(iproj))) cycle
      call projector_end(ep%proj(iproj))
    end do
 
    assert(allocated(ep%proj))
    safe_deallocate_a(ep%proj)
 
    pop_sub(epot_end)
 
  end subroutine epot_end
 
  ! ---------------------------------------------------------
  subroutine epot_bind_poisson_solver(ep, psolver)
    type(epot_t),             intent(inout) :: ep
    type(poisson_t), target,  intent(in)    :: psolver
 
    push_sub(epot_bind_poisson_solver)
 
    if (ep%have_density) then
      ep%poisson_solver => psolver
    else
      nullify(ep%poisson_solver)
    end if
 
    pop_sub(epot_bind_poisson_solver)
  end subroutine epot_bind_poisson_solver
 
  ! ---------------------------------------------------------
  subroutine epot_generate(ep, namespace, mesh, ions, st_d)
    type(epot_t),             intent(inout) :: ep
    type(namespace_t),        intent(in)    :: namespace
    class(mesh_t),    target, intent(in)    :: mesh
    type(ions_t),     target, intent(inout) :: ions
    type(states_elec_dim_t),  intent(inout) :: st_d
 
    integer :: ia
    type(ps_t), pointer :: ps
    logical, allocatable :: spec_checked(:)
 
    call profiling_in("EPOT_GENERATE")
    push_sub(epot_generate)
 
    ! Local part
    ep%vpsl = m_zero
 
    ! we assume that we need to recalculate the ion-ion energy
    call ion_interaction_calculate(ions%ion_interaction, ions%space, ions%latt, ions%atom, &
      ions%natoms, ions%pos, mesh%box%bounding_box_l, ep%eii, ep%fii)
 
    ! the pseudopotential part.
    do ia = 1, ions%natoms
      select type(spec=>ions%atom(ia)%species)
      type is(pseudopotential_t)
        call projector_end(ep%proj(ia))
        call projector_init(ep%proj(ia), spec, namespace, st_d%dim, ep%proj_reltype)
      end select
    end do
 
    do ia = ions%atoms_dist%start, ions%atoms_dist%end
      if (ep%proj(ia)%type == proj_none) cycle
      select type(spec=>ions%atom(ia)%species)
      type is(pseudopotential_t)
        ps => spec%ps
        call submesh_init(ep%proj(ia)%sphere, ions%space, mesh, ions%latt, ions%pos(:, ia), ps%rc_max)
      end select
    end do
 
    if (ions%atoms_dist%parallel) then
      do ia = 1, ions%natoms
        if (ep%proj(ia)%type == proj_none) cycle
        select type(spec=>ions%atom(ia)%species)
        type is(pseudopotential_t)
          ps => spec%ps
          call submesh_broadcast(ep%proj(ia)%sphere, ions%space, mesh, ions%pos(:, ia), ps%rc_max, &
            ions%atoms_dist%process(ia), ions%atoms_dist%mpi_grp)
        end select
      end do
    end if
 
    do ia = 1, ions%natoms
      select type(spec=>ions%atom(ia)%species)
      type is(pseudopotential_t)
        call projector_build(ep%proj(ia), spec, ep%so_strength)
        if (.not. projector_is(ep%proj(ia), proj_none)) ep%non_local = .true.
      end select
    end do
 
    ! Check once, for one atom of each species, that the projectors are properly
    ! resolved on the grid, to warn about possible spurious (ghost) states
    if (.not. ep%proj_check_done) then
      safe_allocate(spec_checked(1:ions%nspecies))
      spec_checked = .false.
      do ia = 1, ions%natoms
        if (ep%proj(ia)%type == proj_none) cycle
        if (spec_checked(ions%atom(ia)%species%get_index())) cycle
        select type(spec=>ions%atom(ia)%species)
        type is(pseudopotential_t)
          call projector_check_discretization(spec, ep%proj(ia)%sphere, namespace)
          spec_checked(ions%atom(ia)%species%get_index()) = .true.
        end select
      end do
      safe_deallocate_a(spec_checked)
      ep%proj_check_done = .true.
    end if
 
    pop_sub(epot_generate)
    call profiling_out("EPOT_GENERATE")
  end subroutine epot_generate
 
  ! ---------------------------------------------------------
 
  logical pure function local_potential_has_density(space, species) result(has_density)
    class(space_t),       intent(in)    :: space
    class(species_t),     intent(in)    :: species
 
    has_density = species%has_density .or. (species%is_ps() .and. space%is_periodic())
 
  end function local_potential_has_density
 
  ! ---------------------------------------------------------
  subroutine epot_local_potential(ep, namespace, space, latt, mesh, species, pos, iatom, vpsl)
    type(epot_t),             intent(in)    :: ep
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(lattice_vectors_t),  intent(in)    :: latt
    class(mesh_t),            intent(in)    :: mesh
    class(species_t), target, intent(in)    :: species
    real(real64),             intent(in)    :: pos(1:space%dim)
    integer,                  intent(in)    :: iatom
    real(real64), contiguous,        intent(inout) :: vpsl(:)
 
    integer :: ip
    real(real64) :: radius
    real(real64), allocatable :: vl(:), rho(:)
    type(submesh_t)  :: sphere
    type(ps_t), pointer :: ps
 
    push_sub(epot_local_potential)
    call profiling_in("EPOT_LOCAL")
 
    if (ep%local_potential_precalculated) then
 
      call lalg_axpy(mesh%np, m_one, ep%local_potential(:, iatom), vpsl)
 
    else
 
      !Local potential, we can get it by solving the Poisson equation
      !(for all-electron species or pseudopotentials in periodic
      !systems) or by applying it directly to the grid
      safe_allocate(vl(1:mesh%np))
 
      if (local_potential_has_density(space, species)) then
        safe_allocate(rho(1:mesh%np))
 
        call species_get_long_range_density(species, namespace, space, latt, pos, mesh, rho, sphere)
 
        call dpoisson_solve(ep%poisson_solver, namespace, vl, rho, all_nodes = .false.)
 
        safe_deallocate_a(rho)
 
      else
 
        call species_get_local(species, namespace, space, latt, pos, mesh, vl)
 
      end if
 
      call lalg_axpy(mesh%np, m_one, vl, vpsl)
      safe_deallocate_a(vl)
 
      !the localized part
      select type(species)
      class is(pseudopotential_t)
 
        ps => species%ps
 
        radius = min(ps%vl%x_threshold*1.05_real64, spline_range_max(ps%vl))
        if (.not. submesh_compatible(sphere, radius, pos, minval(mesh%spacing(1:space%dim)))) then
          call submesh_end(sphere)
          call submesh_init(sphere, space, mesh, latt, pos, radius)
        end if
        safe_allocate(vl(1:sphere%np))
        vl = m_zero
 
        do ip = 1, sphere%np
          if(sphere%r(ip) <= radius) then
            vl(ip) = spline_eval(ps%vl, sphere%r(ip))
          end if
        end do
 
        call submesh_add_to_mesh(sphere, vl, vpsl)
 
        safe_deallocate_a(vl)
        nullify(ps)
 
      end select
      call submesh_end(sphere)
 
    end if
 
    call profiling_out("EPOT_LOCAL")
    pop_sub(epot_local_potential)
  end subroutine epot_local_potential
 
  ! ---------------------------------------------------------
  subroutine epot_precalc_local_potential(ep, namespace, gr, ions)
    type(epot_t),      intent(inout) :: ep
    type(namespace_t), intent(in)    :: namespace
    type(grid_t),      intent(in)    :: gr
    type(ions_t),      intent(in)    :: ions
 
    integer :: iatom
 
    push_sub(epot_precalc_local_potential)
 
    if (.not. allocated(ep%local_potential)) then
      safe_allocate(ep%local_potential(1:gr%np, 1:ions%natoms))
    end if
 
    ep%local_potential_precalculated = .false.
 
    do iatom = 1, ions%natoms
      ep%local_potential(1:gr%np, iatom) = m_zero
      call epot_local_potential(ep, namespace, ions%space, ions%latt, gr, ions%atom(iatom)%species, &
        ions%pos(:, iatom), iatom, ep%local_potential(1:gr%np, iatom))!, time)
    end do
    ep%local_potential_precalculated = .true.
 
    pop_sub(epot_precalc_local_potential)
  end subroutine epot_precalc_local_potential
 
  ! ---------------------------------------------------------
 
  logical function epot_have_external_potentials(ep)
    type(epot_t), intent(in)  :: ep
 
    push_sub(epot_have_external_potentials)
 
    epot_have_external_potentials = allocated(ep%e_field)
 
    pop_sub(epot_have_external_potentials)
 
  end function epot_have_external_potentials
 
end module epot_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: