orbitalbasis.F90

Go to the documentation of this file.

!! Copyright (C) 2016 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 orbitalbasis_oct_m
  use accel_oct_m
  use allelectron_oct_m
  use atomic_orbital_oct_m
  use debug_oct_m
  use distributed_oct_m
  use global_oct_m
  use io_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use namespace_oct_m
  use orbitalset_oct_m
  use orbitalset_utils_oct_m
  use parser_oct_m
  use profiling_oct_m
  use species_oct_m
  use species_factory_oct_m
  use submesh_oct_m
  use types_oct_m
 
  implicit none
 
  private
 
  public ::                         &
    orbitalbasis_t,                 &
    orbitalbasis_init,              &
    orbitalbasis_end,               &
    dorbitalbasis_build,            &
    zorbitalbasis_build,            &
    dorbitalbasis_build_empty,      &
    zorbitalbasis_build_empty
 
  type orbitalbasis_t
    private
    type(orbitalset_t), allocatable, public :: orbsets(:)
 
    integer,            public :: norbsets = 0
    integer,            public :: maxnorbs = 0
    integer,            public :: max_np   = 0
    integer,            public :: size     = 0
    integer, allocatable, public :: global2os(:,:)
    integer, allocatable         :: os2global(:,:)
 
    integer(int64)              :: truncation
    real(real64)               :: threshold = 0.01_real64 
 
    logical,            public :: normalize = .true.          
    logical,            public :: use_submesh   = .false.     
    logical,            public :: orthogonalization = .false. 
 
    character(len=256), public :: debugdir
 
    logical,            public :: combine_j_orbitals
  end type orbitalbasis_t
 
contains
 
  subroutine orbitalbasis_init(this, namespace, periodic_dim)
    type(orbitalbasis_t),    intent(inout) :: this
    type(namespace_t),       intent(in)    :: namespace
    integer,                 intent(in)    :: periodic_dim
 
    push_sub(orbitalbasis_init)
 
    !%Variable AOTruncation
    !%Type flag
    !%Default ao_full
    !%Section Atomic Orbitals
    !%Description
    !% This option determines how Octopus will truncate the orbitals used for DFT+U.
    !% Except for the full method, the other options are only there to get a quick idea.
    !%Option ao_full bit(0)
    !% The full size of the orbitals used. The radius is controled by variable AOThreshold.
    !%Option ao_box bit(1)
    !% The radius of the orbitals are restricted to the size of the simulation box.
    !% This reduces the number of points used to discretize the orbitals.
    !% This is mostly a debug option, and one should be aware that changing the size of the simulation box
    !% will affect the result of the calculation. It is recommended to use ao_nlradius instead.
    !%Option ao_nlradius bit(2)
    !% The radius of the orbitals are restricted to the radius of the non-local part of the pseudopotential
    !% of the corresponding atom.
    !%End
    call parse_variable(namespace, 'AOTruncation', option__aotruncation__ao_full, this%truncation)
    call messages_print_var_option('AOTruncation', this%truncation, namespace=namespace)
 
    !%Variable AOThreshold
    !%Type float
    !%Default 0.01
    !%Section Atomic Orbitals
    !%Description
    !% Determines the threshold used to compute the radius of the atomic orbitals for DFT+U and for Wannier90.
    !% This radius is computed by making sure that the
    !% absolute value of the radial part of the atomic orbital is below the specified threshold.
    !% This value should be converged to be sure that results do not depend on this value.
    !% However increasing this value increases the number of grid points covered by the orbitals and directly affect performances.
    !%End
    call parse_variable(namespace, 'AOThreshold', 0.01_real64, this%threshold)
    if (this%threshold <= m_zero) call messages_input_error(namespace, 'AOThreshold')
    call messages_print_var_value('AOThreshold', this%threshold, namespace=namespace)
 
    !%Variable AONormalize
    !%Type logical
    !%Default yes
    !%Section Atomic Orbitals
    !%Description
    !% If set to yes, Octopus will normalize the atomic orbitals individually.
    !% This variable is ignored is <tt>AOLoewdin</tt> is set to yes.
    !%End
    call parse_variable(namespace, 'AONormalize', .true., this%normalize)
    call messages_print_var_value('AONormalize', this%normalize, namespace=namespace)
 
    !%Variable AOSubmesh
    !%Type logical
    !%Section Atomic Orbitals
    !%Description
    !% If set to yes, Octopus will use submeshes to internally store the orbitals with
    !% their phase instead of storing them on the mesh. This is usually slower for small
    !% periodic systems, but becomes advantageous for large supercells.
    !% Submeshes are not compatible with Loewdin orthogonalization.
    !% For periodic systems, the default is set to no, whereas it is set to yes for isolated systems.
    !%End
    if (periodic_dim > 0) then
      call parse_variable(namespace, 'AOSubmesh', .false., this%use_submesh)
    else
      call parse_variable(namespace, 'AOSubmesh', .true., this%use_submesh)
    end if
    call messages_print_var_value('AOSubmesh', this%use_submesh, namespace=namespace)
 
    !%Variable AOLoewdin
    !%Type logical
    !%Default no
    !%Section Atomic Orbitals
    !%Description
    !% This option determines if the atomic orbital basis is orthogonalized or not.
    !% This is done for using the Loewdin orthogonalization scheme.
    !% The default is set to no for the moment as this option is
    !% not yet implemented for isolated systems, and seems to lead to important egg-box effect
    !%End
    call parse_variable(namespace, 'AOLoewdin', .false., this%orthogonalization)
    call messages_print_var_value('AOLoewdin', this%orthogonalization, namespace=namespace)
    if (this%orthogonalization) call messages_experimental("AOLoewdin")
    if (this%orthogonalization) this%normalize = .false.
 
    !%Variable AOCombineJOrbitals
    !%Type logical
    !%Default no
    !%Section Atomic Orbitals
    !%Description
    !% By default, Octopus creates two sets of atomic orbitals for j-dependent pseudopotentials.
    !% When doing ACBN0 functional, this leads to a Hubbard U for j=l-1/2 and one for j=l+1/2, unless
    !% hubbard_j is specified in the species block.
    !% By setting this variable to yes, one can instead define a single set of atomic orbitals that
    !% contains both the j=l-1/2 and j=l+1/2 orbitals.
    !% This is only relevant for spinor calculations with j-dependent pseudopotentials.
    !%End
    call parse_variable(namespace, 'AOCombineJOrbitals', .false., this%combine_j_orbitals)
    call messages_print_var_value('AOCombineJOrbitals', this%combine_j_orbitals, namespace=namespace)
 
 
    if (this%orthogonalization .and. this%use_submesh) then
      call messages_not_implemented("AOLoewdin=yes with AOSubmesh=yes.")
    end if
 
    if (debug%info) then
      write(this%debugdir, '(a)') 'debug/ao_basis'
      call io_mkdir(this%debugdir, namespace)
    end if
 
    pop_sub(orbitalbasis_init)
  end subroutine orbitalbasis_init
 
 
  subroutine orbitalbasis_end(this)
    type(orbitalbasis_t), intent(inout) :: this
 
    integer :: ios
 
    push_sub(orbitalbasis_end)
 
    do ios = 1, this%norbsets
      call orbitalset_end(this%orbsets(ios))
    end do
    safe_deallocate_a(this%orbsets)
 
    safe_deallocate_a(this%global2os)
    safe_deallocate_a(this%os2global)
 
    pop_sub(orbitalbasis_end)
  end subroutine orbitalbasis_end
 
  ! ---------------------------------------------------------
  subroutine orbitalbasis_set_nlmj_radius(this, mesh, os, hubbard_l, hubbard_j, norb)
    type(orbitalbasis_t),    intent(in)    :: this
    type(mesh_t),            intent(in)    :: mesh
    type(orbitalset_t),      intent(inout) :: os
    integer,                 intent(in)    :: hubbard_l
    real(real64),            intent(in)    :: hubbard_j
    integer,                 intent(out)   :: norb
 
    integer :: iorb, i, l, m, n
    real(real64) :: j
 
    norb = 0
 
    push_sub(orbitalbasis_set_nlmj_radius)
 
    do iorb = 1, os%spec%get_niwfs()
      call os%spec%get_iwf_ilm(iorb, 1, i, l, m)
      call os%spec%get_iwf_n(iorb, 1, n)
      call os%spec%get_iwf_j(iorb, j)
      if (l .eq. hubbard_l .and. is_close(hubbard_j, j)) then
        norb = norb + 1
        call orbitalset_set_jln(os, j, hubbard_l, n)
        os%ii = i
        if (os%spec%is_full()) os%ii = n
        os%radius = atomic_orbital_get_radius(os%spec, mesh, iorb, 1, &
          this%truncation, this%threshold)
      end if
    end do
    pop_sub(orbitalbasis_set_nlmj_radius)
  end subroutine orbitalbasis_set_nlmj_radius
 
 
#include "undef.F90"
#include "real.F90"
#include "orbitalbasis_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "orbitalbasis_inc.F90"
 
end module orbitalbasis_oct_m