curv_gygi.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2025 S. Ohlmann
!!
!! 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 curv_gygi_oct_m
  use coordinate_system_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use root_solver_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                   &
    curv_gygi_t,              &
    curv_gygi_copy,           &
    curv_gygi_hessian
 
  type, extends(coordinate_system_t) :: curv_gygi_t
    private
    real(real64), public :: A
    real(real64), public :: alpha
    real(real64), public :: beta
    real(real64), allocatable :: pos(:, :)
    integer :: npos
    type(root_solver_t) :: rs
  contains
    procedure :: to_cartesian => curv_gygi_to_cartesian
    procedure :: from_cartesian => curv_gygi_from_cartesian
    procedure :: write_info => curv_gygi_write_info
    procedure :: surface_element => curv_gygi_surface_element
    procedure :: jacobian => curv_gygi_jacobian
    procedure :: jacobian_inverse => curv_gygi_jacobian_inverse
    procedure :: trace_hessian => curv_gygi_trace_hessian
    final :: curv_gygi_finalize
  end type curv_gygi_t
 
  interface curv_gygi_t
    procedure curv_gygi_constructor
  end interface curv_gygi_t
 
  ! Auxiliary variables for the root solver.
  class(curv_gygi_t), pointer  :: gygi_p
  type(curv_gygi_t), target :: gygi_global
  real(real64), allocatable :: chi_p(:)
 
contains
 
  ! ---------------------------------------------------------
  function curv_gygi_constructor(namespace, dim, npos, pos) result(gygi)
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: dim
    integer,           intent(in)  :: npos
    real(real64),      intent(in)  :: pos(1:dim,1:npos)
    class(curv_gygi_t), pointer :: gygi
 
    push_sub(curv_gygi_constructor)
 
    safe_allocate(gygi)
 
    gygi%dim = dim
    gygi%local_basis = .true.
    gygi%orthogonal = .true. ! This needs to be checked.
 
    gygi%npos = npos
    safe_allocate(gygi%pos(1:dim, 1:gygi%npos))
    gygi%pos(:, :) = pos(:, :)
 
    !%Variable CurvGygiA
    !%Type float
    !%Default 0.5
    !%Section Mesh::Curvilinear::Gygi
    !%Description
    !% The grid spacing is reduced locally around each atom, and the reduction is
    !% given by 1/(1+<i>A</i>), where <i>A</i> is specified by this variable. So, if
    !% <i>A</i>=1/2 (the default), the grid spacing is reduced to two thirds = 1/(1+1/2).
    !% [This is the <math>A_{\alpha}</math> variable in Eq. 2 of F. Gygi and G. Galli, <i>Phys.
    !% Rev. B</i> <b>52</b>, R2229 (1995)]. It must be larger than zero.
    !% It must not be larger than two, as the coordinate transformation is not bijective
    !% anymore for larger values of <i>A</i>. This caps the local refinement that can be
    !% obtained at 1/(1+2), i.e. a factor of three, which is the governing limitation when
    !% attempting all-electron calculations on heavy atoms.
    !%End
    call parse_variable(namespace, 'CurvGygiA', m_half, gygi%A)
 
    !%Variable CurvGygiAlpha
    !%Type float
    !%Default 2.0 a.u.
    !%Section Mesh::Curvilinear::Gygi
    !%Description
    !% This number determines the region over which the grid is enhanced (range of
    !% enhancement of the resolution). That is, the grid is enhanced on a sphere
    !% around each atom, whose radius is given by this variable. [This is the <math>a_{\alpha}</math>
    !% variable in Eq. 2 of F. Gygi and G. Galli, <i>Phys. Rev. B</i> <b>52</b>, R2229 (1995)].
    !% It must be larger than zero.
    !%End
 
    call parse_variable(namespace, 'CurvGygiAlpha', m_two, gygi%alpha, units_inp%length)
    !%Variable CurvGygiBeta
    !%Type float
    !%Default 4.0 a.u.
    !%Section Mesh::Curvilinear::Gygi
    !%Description
    !% This number determines the distance over which Euclidean coordinates are
    !% recovered. [This is the <math>b_{\alpha}</math> variable in Eq. 2 of F. Gygi and G. Galli,
    !% <i>Phys. Rev. B</i> <b>52</b>, R2229 (1995)]. It must be larger than zero.
    !%End
    call parse_variable(namespace, 'CurvGygiBeta', m_four, gygi%beta, units_inp%length)
 
    if (gygi%a <= m_zero)     call messages_input_error(namespace, 'CurvGygiA')
    if (gygi%a > m_two)       call messages_input_error(namespace, 'CurvGygiA')
    if (gygi%alpha <= m_zero) call messages_input_error(namespace, 'CurvGygiAlpha')
    if (gygi%beta <= m_zero)  call messages_input_error(namespace, 'CurvGygiBeta')
 
    gygi%min_mesh_scaling_product = (m_one / (m_one + gygi%A))**gygi%dim
 
    ! initialize root solver
    call root_solver_init(gygi%rs, namespace, dim, solver_type = root_newton, maxiter = 500, abs_tolerance = 1.0e-10_real64)
 
    pop_sub(curv_gygi_constructor)
  end function curv_gygi_constructor
 
  ! ---------------------------------------------------------
  subroutine curv_gygi_copy(this_out, this_in)
    type(curv_gygi_t), intent(inout) :: this_out
    type(curv_gygi_t), intent(in)    :: this_in
 
    push_sub(curv_gygi_copy)
 
    this_out%A = this_in%A
    this_out%alpha = this_in%alpha
    this_out%beta = this_in%beta
    safe_allocate_source_a(this_out%pos, this_in%pos)
    this_out%pos = this_in%pos
    this_out%npos = this_in%npos
    this_out%dim = this_in%dim
    this_out%local_basis = this_in%local_basis
    this_out%orthogonal = this_in%orthogonal
    call root_solver_init(this_out%rs, global_namespace, this_in%dim, solver_type = root_newton, &
      maxiter = 500, abs_tolerance = 1.0e-10_real64)
 
    pop_sub(curv_gygi_copy)
  end subroutine curv_gygi_copy
 
  ! ---------------------------------------------------------
  subroutine curv_gygi_finalize(this)
    type(curv_gygi_t), intent(inout) :: this
 
    push_sub(curv_gygi_finalize)
 
    safe_deallocate_a(this%pos)
 
    pop_sub(curv_gygi_finalize)
  end subroutine curv_gygi_finalize
 
  ! ---------------------------------------------------------
  real(real64) pure function gygi_f(this, r)
    class(curv_gygi_t), intent(in) :: this
    real(real64),       intent(in) :: r
 
    ! no PUSH_SUB, called too often
 
    if (r < 1.0e-4_real64) then
      ! below 1e-4, the relative deviation of this Taylor expansion from the function is less than 1e-15
      gygi_f = this%A * (m_one + (-m_one/(m_three*this%alpha**2) - m_one/this%beta**2)*r**2)
    else if (-(r/this%beta)**2 <= m_min_exp_arg) then
      gygi_f = m_zero
    else
      gygi_f = this%A * this%alpha/r * tanh(r/this%alpha) * exp(-(r/this%beta)**2)
    end if
  end function gygi_f
 
  ! ---------------------------------------------------------
  ! absorb division by r into the function to get the correct expansion
  real(real64) pure function gygi_dfdr_over_r(this, r)
    class(curv_gygi_t), intent(in) :: this
    real(real64),       intent(in) :: r
 
    ! no PUSH_SUB, called too often
 
    if (r < 3.0e-3_real64) then
      ! below 3e-3, the expansion is closer than 1e-13 to the function and does not suffer from
      ! rounding errors
      gygi_dfdr_over_r = this%A * ((-m_two/(m_three*this%alpha**2) - m_two/this%beta**2) + &
        (8._real64/(15._real64*this%alpha**4) + m_four/(m_three*this%alpha**2*this%beta**2) &
        + m_two/this%beta**4) * r**2)
    else if (-(r/this%beta)**2 <= m_min_exp_arg) then
      gygi_dfdr_over_r = m_zero
    else
      gygi_dfdr_over_r = this%A * (this%beta**2*r/cosh(r/this%alpha)**2 - &
        this%alpha*(this%beta**2+m_two*r**2) * tanh(r/this%alpha)) &
        / (this%beta**2*r**3) * exp(-(r/this%beta)**2)
    end if
  end function gygi_dfdr_over_r
 
  ! ---------------------------------------------------------
  ! only a combination of the second derivative is needed, namely
  ! (d2f/d2r - 1/r df/dr)/r**2
  real(real64) pure function gygi_d2fdr2_combination(this, r)
    class(curv_gygi_t), intent(in) :: this
    real(real64),       intent(in) :: r
 
    ! no PUSH_SUB, called too often
 
    if (r < 3.0e-3_real64) then
      ! below 3e-3, the expansion is closer than 1e-5 to the function and does not suffer from
      ! rounding errors
      gygi_d2fdr2_combination = this%A * ( &
        m_two*(8._real64/(15._real64*this%alpha**4) + m_four/(m_three*this%alpha**2*this%beta**2) &
        + m_two/this%beta**4) &
        - m_four*(34._real64/(105._real64*this%alpha**6) + m_four/(m_five*this%alpha**4*this%beta**2) &
        + m_one/(this%alpha**2*this%beta**4) + m_one/this%beta**6) * r**2 &
        )
    else if (-(r/this%beta)**2 <= m_min_exp_arg) then
      gygi_d2fdr2_combination = m_zero
    else
      gygi_d2fdr2_combination = this%A * exp(-(r/this%beta)**2)*( &
        -m_two*tanh(r/this%alpha)/(cosh(r/this%alpha)**2*this%alpha*r) &
        + tanh(r/this%alpha)*(m_three*this%alpha/r**3+m_four*this%alpha/(this%beta**2*r)+m_four*this%alpha*r/this%beta**4) &
        + m_one/cosh(r/this%alpha)**2*(-m_three/r**2-m_four/this%beta**2))/r**2
    end if
  end function gygi_d2fdr2_combination
 
 
  ! ---------------------------------------------------------
  function curv_gygi_to_cartesian(this, chi) result(xx)
    class(curv_gygi_t), target, intent(in)  :: this
    real(real64),               intent(in)  :: chi(:)
    real(real64) :: xx(1:this%dim), xx_start(1:this%dim)
 
    logical :: conv
    integer :: i_conv, n_conv
 
    ! no PUSH_SUB, called too often
 
    gygi_p => this
    safe_allocate(chi_p(1:this%dim))
    chi_p(:) = chi(:)
 
    call droot_solver_run(this%rs, getf, xx, conv, startval = chi)
    nullify(gygi_p)
 
    if (.not. conv) then
      call curv_gygi_copy(gygi_global, this)
      gygi_p => gygi_global
      ! try to converge with decreasing A
      n_conv = 8
      do i_conv = 1, n_conv
        gygi_p%A = gygi_p%A / m_two
        call droot_solver_run(gygi_p%rs, getf, xx, conv, startval = chi)
        if (conv) then
          exit
        end if
      end do
      if (conv) then
        ! increase A again and take previous xx as starting point
        n_conv = i_conv
        do i_conv = n_conv, 1, -1
          gygi_p%A = gygi_p%A * m_two
          xx_start = xx
          call droot_solver_run(gygi_p%rs, getf, xx, conv, startval = xx_start)
        end do
      end if
      nullify(gygi_p)
      call curv_gygi_finalize(gygi_global)
    end if
 
    safe_deallocate_a(chi_p)
 
    if (.not. conv) then
      message(1) = "During the construction of the adaptive grid, the Newton-Raphson"
      message(2) = "method did not converge for point:"
      write(message(3),'(9f14.6)') xx(1:this%dim)
      message(4) = "Try varying the Gygi parameters -- usually reducing CurvGygiA or"
      message(5) = "CurvGygiAlpha (or both) solves the problem."
      call messages_fatal(5)
    end if
 
  end function curv_gygi_to_cartesian
 
  ! ---------------------------------------------------------
  pure function curv_gygi_from_cartesian(this, xx) result(chi)
    class(curv_gygi_t), target, intent(in)  :: this
    real(real64),               intent(in)  :: xx(:)
    real(real64) :: chi(1:this%dim)
 
    integer :: i, ia
    real(real64) :: diff(this%dim), f
 
    ! no PUSH_SUB, called too often
 
    chi(1:this%dim) = xx(1:this%dim)
    do ia = 1, this%npos
      diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
      f = gygi_f(this, norm2(diff))
      do i = 1, this%dim
        chi(i) = chi(i) + diff(i) * f
      end do
    end do
 
  end function curv_gygi_from_cartesian
 
  ! ---------------------------------------------------------
  subroutine curv_gygi_write_info(this, iunit, namespace)
    class(curv_gygi_t),           intent(in) :: this
    integer,            optional, intent(in) :: iunit
    type(namespace_t),  optional, intent(in) :: namespace
 
    push_sub(curv_gygi_write_info)
 
    write(message(1), '(a)')  '  Curvilinear Method = gygi'
    write(message(2), '(a)')  '  Gygi Parameters:'
    write(message(3), '(4x,a,f6.3)')  'A = ', this%a
    write(message(4), '(4x,3a,f6.3)') 'alpha [', trim(units_abbrev(units_out%length)), '] = ', &
      units_from_atomic(units_out%length, this%alpha)
    write(message(5), '(4x,3a,f6.3)') 'beta  [', trim(units_abbrev(units_out%length)), '] = ', &
      units_from_atomic(units_out%length, this%beta)
    call messages_info(5, iunit=iunit, namespace=namespace)
 
    pop_sub(curv_gygi_write_info)
  end subroutine curv_gygi_write_info
 
  ! ---------------------------------------------------------
  real(real64) function curv_gygi_surface_element(this, idir) result(ds)
    class(curv_gygi_t), intent(in) :: this
    integer,            intent(in) :: idir
 
    ds = m_zero
    message(1) = 'Surface element with gygi curvilinear coordinates not implemented'
    call messages_fatal(1)
 
  end function curv_gygi_surface_element
 
  ! ---------------------------------------------------------
  function curv_gygi_jacobian(this, chi) result(jacobian)
    class(curv_gygi_t), intent(in)  :: this
    real(real64),       intent(in)  :: chi(:)
    real(real64) :: jacobian(1:this%dim, 1:this%dim)
 
    jacobian = this%jacobian_inverse(chi)
    call lalg_inverse(this%dim, jacobian, "dir")
  end function curv_gygi_jacobian
 
  ! ---------------------------------------------------------
  function curv_gygi_jacobian_inverse(this, chi) result(jacobian_inverse)
    class(curv_gygi_t), intent(in)  :: this
    real(real64),       intent(in)  :: chi(:)
    real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
 
    real(real64) :: xx(1:this%dim)
 
    ! no PUSH_SUB, called too often
 
    xx(:) = this%to_cartesian(chi)
    jacobian_inverse = curv_gygi_jacobian_inverse_cartesian(this, xx)
 
  end function curv_gygi_jacobian_inverse
 
  ! ---------------------------------------------------------
  function curv_gygi_jacobian_inverse_cartesian(this, xx) result(jacobian_inverse)
    class(curv_gygi_t), intent(in)  :: this
    real(real64),       intent(in)  :: xx(:)
    real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
 
    integer :: i, ix, iy
    real(real64) :: r, diff(1:this%dim), dfdr_over_r, f
 
    ! no PUSH_SUB, called too often
 
    jacobian_inverse(1:this%dim, 1:this%dim) = diagonal_matrix(this%dim, m_one)
 
    do i = 1, this%npos
      diff = xx(1:this%dim) - this%pos(1:this%dim, i)
      r = norm2(diff)
 
      f = gygi_f(this, r)
      dfdr_over_r = gygi_dfdr_over_r(this, r)
 
      do ix = 1, this%dim
        jacobian_inverse(ix, ix) = jacobian_inverse(ix, ix) + f
        do iy = 1, this%dim
          jacobian_inverse(ix, iy) = jacobian_inverse(ix, iy) + diff(ix)*diff(iy)*dfdr_over_r
        end do
      end do
    end do
 
  end function curv_gygi_jacobian_inverse_cartesian
 
 
  ! ---------------------------------------------------------
  pure subroutine curv_gygi_hessian(this, xx, hessian, natoms)
    class(curv_gygi_t), intent(in)  :: this
    real(real64),       intent(in)  :: xx(:)
    real(real64),       intent(out) :: hessian(:, :, :)
    integer,  optional, intent(in)  :: natoms
 
    integer :: ia, i, j, k, natoms_
    real(real64) :: r, diff(this%dim), dfdr_over_r, d2fdr2_combination
 
    ! no PUSH_SUB, called too often
 
    hessian(1:this%dim, 1:this%dim, 1:this%dim) = m_zero
 
    natoms_ = this%npos
    if (present(natoms)) natoms_ = natoms
 
    do ia = 1, natoms_
      diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
      r = norm2(diff)
 
      dfdr_over_r = gygi_dfdr_over_r(this, r)
      d2fdr2_combination = gygi_d2fdr2_combination(this, r)
 
      do i = 1, this%dim
        do j = 1, this%dim
          do k = 1, this%dim
            if (i == j) then
              hessian(i, j, k) = hessian(i, j, k) + diff(k)*dfdr_over_r
            end if
            if (i == k) then
              hessian(i, j, k) = hessian(i, j, k) + diff(j)*dfdr_over_r
            end if
            if (j == k) then
              hessian(i, j, k) = hessian(i, j, k) + diff(i)*dfdr_over_r
            end if
            hessian(i, j, k) = hessian(i, j, k) + diff(i)*diff(j)*diff(k)*d2fdr2_combination
          end do
        end do
      end do
    end do
 
  end subroutine curv_gygi_hessian
 
  ! ---------------------------------------------------------
  function curv_gygi_trace_hessian(this, chi) result(trace_hessian)
    class(curv_gygi_t), intent(in)  :: this
    real(real64),       intent(in)  :: chi(:)
    real(real64) :: trace_hessian(1:this%dim)
 
    integer :: ia, i, k
    real(real64) :: r, xx(this%dim), diff(this%dim), dfdr_over_r, d2fdr2_combination
 
    ! no PUSH_SUB, called too often
 
    xx(:) = this%to_cartesian(chi)
    trace_hessian(:) = m_zero
 
    do ia = 1, this%npos
      diff = xx(1:this%dim) - this%pos(1:this%dim, ia)
      r = norm2(diff)
 
      dfdr_over_r = gygi_dfdr_over_r(this, r)
      d2fdr2_combination = gygi_d2fdr2_combination(this, r)
 
      do i = 1, this%dim
        do k = 1, this%dim
          if (i == k) then
            trace_hessian(i) = trace_hessian(i) + m_two*diff(k)*dfdr_over_r
          end if
          trace_hessian(i) = trace_hessian(i) + diff(i)*dfdr_over_r
          trace_hessian(i) = trace_hessian(i) + diff(i)*diff(k)*diff(k)*d2fdr2_combination
        end do
      end do
    end do
 
  end function curv_gygi_trace_hessian
 
  ! ---------------------------------------------------------
  subroutine getf(y, f, jf)
    real(real64), intent(in)    :: y(:)
    real(real64), intent(out)   :: f(:), jf(:, :)
 
    ! no PUSH_SUB, called too often
 
    jf = curv_gygi_jacobian_inverse_cartesian(gygi_p, y)
    f = gygi_p%from_cartesian(y)
    f(1:gygi_p%dim) = f(1:gygi_p%dim) - chi_p(1:gygi_p%dim)
 
  end subroutine getf
 
end module curv_gygi_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: