atomic_orbital.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 atomic_orbital_oct_m
  use allelectron_oct_m
  use box_cylinder_oct_m
  use box_minimum_oct_m
  use box_sphere_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use jellium_oct_m
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use orbitalset_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use submesh_oct_m
 
  implicit none
 
  private
 
  public ::                           &
    atomic_orbital_get_radius,        &
    datomic_orbital_get_submesh,      &
    zatomic_orbital_get_submesh,      &
    datomic_orbital_get_submesh_safe, &
    zatomic_orbital_get_submesh_safe, &
    atomic_orbital_init_submesh,      &
    dget_atomic_orbital,              &
    zget_atomic_orbital,              &
    l_notation
 
  character(len=1), parameter :: &
    l_notation(0:3) = (/ 's', 'p', 'd', 'f' /)
 
  integer, public, parameter :: MAX_L = 4
 
contains
 
  ! ---------------------------------------------------------
  real(real64) function atomic_orbital_get_radius(species, mesh, iorb, ispin, truncation, threshold) result(radius)
    class(species_t),         intent(in)   :: species
    type(mesh_t),             intent(in)   :: mesh
    integer,                  intent(in)   :: iorb, ispin
    integer(int64),           intent(in)   :: truncation
    real(real64),             intent(in)   :: threshold
 
    integer :: ii, ll, mm
 
    push_sub(atomic_orbital_get_radius)
 
    call species%get_iwf_ilm(iorb, ispin, ii, ll, mm)
 
    if (truncation == option__aotruncation__ao_full) then
      radius = species%get_iwf_radius(ii, ispin, threshold)
    else
      radius = species%get_iwf_radius(ii, ispin)
 
      if (truncation == option__aotruncation__ao_box) then
        ! if the orbital is larger than the size of the box, we restrict it to this size,
        ! otherwise the orbital will overlap more than one time with the simulation box.
        ! This would induces phase problem if the complete mesh is used instead of the sphere
        radius = min(radius, minval(mesh%box%bounding_box_l(1:mesh%box%dim)-mesh%spacing(1:mesh%box%dim)*1.01_real64))
      else
        !If asked, we truncate the orbital to the radius on the projector spheres
        !of the NL part of the pseudopotential.
        !This is a way to garanty no overlap between orbitals of different atoms.
        select type(species)
        type is(pseudopotential_t)
          radius = min(radius, species%get_radius())
        end select
      end if
 
    end if
    ! make sure that if the spacing is too large, the orbitals fit in a few points at least
    radius = max(radius, m_two*maxval(mesh%spacing))
 
    pop_sub(atomic_orbital_get_radius)
  end function atomic_orbital_get_radius
 
  ! ---------------------------------------------------------
  subroutine atomic_orbital_init_submesh (namespace, space, latt, center, mesh, sm, radius)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(lattice_vectors_t),  intent(in)    :: latt
    real(real64),             intent(in)    :: center(:)
    class(mesh_t),            intent(in)    :: mesh
    type(submesh_t),          intent(inout) :: sm
    real(real64),             intent(in)    :: radius
 
    push_sub(atomic_orbital_init_submesh)
 
    if (sm%np == -1) then
 
      select type (box => mesh%box)
      type is (box_minimum_t)
        if (radius > box%radius) then
          message(1) = "The radius of an orbital set is bigger than the radius of the simulation box."
          message(2) = "Increase the value of Radius or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
 
      type is (box_sphere_t)
        if (radius > box%radius) then
          message(1) = "The radius of an orbital set is bigger than the radius of the simulation box."
          message(2) = "Increase the value of Radius or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
        if (norm2(center - box%center) + radius > box%radius) then
          message(1) = "An orbital set has points outside of the simulation box."
          message(2) = "Increase the value of Radius or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
 
      type is (box_cylinder_t)
        if (radius > box%radius) then
          message(1) = "The radius of an orbital set is bigger than the radius of the simulation box."
          message(2) = "Increase the value of Radius or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
        if (radius > box%half_length) then
          message(1) = "The radius of an orbital set is bigger than the length of the cylinder box."
          message(2) = "Increase the value of XLength or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
 
        if (norm2(center(2:mesh%box%dim) - box%center(2:mesh%box%dim)) + radius > box%radius) then
          message(1) = "An orbital set has points outside of the simulation box."
          message(2) = "Increase the value of Radius or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
        if (abs(center(1) - box%center(1)) + radius > box%half_length) then
          message(1) = "An orbital set has points outside of the simulation box."
          message(2) = "Increase the value of Xlength or decrease the value of AOThreshold."
          write(message(3),'(a,f8.5,a,i5,a)') 'The value of the radius is ', radius, ' Bohr.'
          call messages_fatal(3, namespace=namespace)
        end if
 
      end select
 
      !We initialise the submesh corresponding to the orbital
      call submesh_init(sm, space, mesh, latt, center, radius)
 
    end if
 
    pop_sub(atomic_orbital_init_submesh)
 
  end subroutine atomic_orbital_init_submesh
 
 
 
#include "undef.F90"
#include "real.F90"
#include "atomic_orbital_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "atomic_orbital_inc.F90"
 
end module atomic_orbital_oct_m