curv_modine.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_modine_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 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_modine_t,              &
    curv_modine_copy
 
  type, extends(coordinate_system_t) :: curv_modine_t
    private
    real(real64), allocatable :: lsize(:)
    real(real64)       :: xbar
    real(real64)       :: Jbar
 
    real(real64), allocatable :: Jlocal(:)
    real(real64), allocatable :: Jrange(:)
 
    real(real64), allocatable :: xi(:,:)
    real(real64), allocatable :: Q(:,:,:)
 
    integer :: npos
    type(root_solver_t) :: rs, rs_init
  contains
    procedure :: to_cartesian => curv_modine_to_cartesian
    procedure :: from_cartesian => curv_modine_from_cartesian
    procedure :: jacobian => curv_modine_jacobian
    procedure :: jacobian_inverse => curv_modine_jacobian_inverse
    procedure :: write_info => curv_modine_write_info
    procedure :: surface_element => curv_modine_surface_element
    final :: curv_modine_finalize
  end type curv_modine_t
 
  interface curv_modine_t
    procedure curv_modine_constructor
  end interface curv_modine_t
 
  ! Auxiliary variables for the root solver. They are threadprivate, as
  ! from_cartesian is called from OpenMP regions. They must not be local
  ! variables of from_cartesian, as getf would then need a trampoline on
  ! the stack, which marks the stack as executable.
  class(curv_modine_t), pointer :: modine_p
  real(real64), allocatable :: x_p(:)
  !$omp threadprivate(modine_p, x_p)
 
contains
 
  ! ---------------------------------------------------------
  function curv_modine_constructor(namespace, dim, npos, pos, lsize, spacing) result(modine)
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: dim
    integer,           intent(in)  :: npos
    real(real64),      intent(in)  :: pos(1:dim,1:npos)
    real(real64),      intent(in)  :: lsize(1:dim)
    real(real64),      intent(in)  :: spacing(1:dim)
    class(curv_modine_t), pointer :: modine
 
    integer :: iatom, n
    real(real64) :: kappa(1:npos)
    type(block_t) :: blk
 
    push_sub(curv_modine_constructor)
 
    safe_allocate(modine)
 
    modine%dim = dim
    modine%local_basis = .true.
    modine%orthogonal = .true. ! This needs to be checked.
 
    modine%npos = npos
 
    safe_allocate(modine%lsize(1:dim))
    modine%lsize = lsize
 
    !%Variable CurvModineXBar
    !%Type float
    !%Default 0.25
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Size of central flat region (in units of <tt>2*Lsize</tt>). Must be between 0 and 1.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineXBar', m_one/m_four, modine%xbar)
 
    !%Variable CurvModineJBar
    !%Type float
    !%Default 0.5
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Increase in density of points is inverse of this parameter.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    call parse_variable(namespace, 'CurvModineJBar', m_half, modine%Jbar)
 
    if (modine%xbar < m_zero .or. modine%xbar > m_one) then
      message(1) = 'The parameter "CurvModineXBar" must lie between 0 and 1.'
      call messages_fatal(1, namespace=namespace)
    end if
    if (modine%Jbar <= modine%xbar) then
      message(1) = 'The parameter "CurvModineJBar" must be larger than the parameter "CurvModineXBar".'
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(modine%Jlocal(1:npos))
    safe_allocate(modine%Jrange(1:npos))
 
    !%Variable CurvModineJlocal
    !%Type block
    !%Default 0.125
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Local refinement around the atoms for each atom. If only one line is specified, the value is taken
    !% for all atoms. Must be between 0 and 1.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    if (parse_block(namespace, "CurvModineJlocal", blk) == 0) then
      n = parse_block_n(blk)
      if (n == 1) then
        call parse_block_float(blk, 0, 0, modine%Jlocal(1))
        modine%Jlocal(:) = modine%Jlocal(1)
      else if (n == npos) then
        do iatom = 1, npos
          call parse_block_float(blk, iatom - 1, 0, modine%Jlocal(iatom))
        end do
      else
        write(message(1), '(3a,i2)') 'Error in block CurvModineJlocal: number of lines must be 1 or ', npos
        call messages_fatal(1, namespace=namespace)
      end if
    else
      ! for default value, see page 55, second paragraph
      modine%Jlocal(:) = 0.125_real64
    end if
 
    do iatom = 1, npos
      if (modine%Jlocal(iatom)<m_zero.or.modine%Jlocal(iatom)>m_one) then
        message(1) = 'The parameter "CurvModineJlocal" must lie between 0 and 1.'
        call messages_fatal(1, namespace=namespace)
      end if
    end do
 
    !%Variable CurvModineJrange
    !%Type block
    !%Default 1.2 b
    !%Section Mesh::Curvilinear::Modine
    !%Description
    !% Local refinement range (a length) for each atom. If only one line is specified, the value is
    !% taken for all atoms.
    !% See N. A. Modine, G. Zumbach, and E. Kaxiras, <i>Phys. Rev. B</i> <b>55</b>, 10289-10301 (1997).
    !%End
    if (parse_block(namespace, "CurvModineJrange", blk) == 0) then
      n = parse_block_n(blk)
      if (n == 1) then
        call parse_block_float(blk, 0, 0, kappa(1), units_inp%length)
        kappa(:) = kappa(1)
      else if (n == npos) then
        do iatom = 1, npos
          call parse_block_float(blk, iatom - 1, 0, kappa(iatom), units_inp%length)
        end do
      else
        write(message(1), '(3a,i2)') 'Error in block CurvModineJrange: number of lines must be 1 or ', npos
        call messages_fatal(1, namespace=namespace)
      end if
    else
      ! for default value, see page 55, second paragraph
      kappa(:) = 1.2_real64
    end if
 
    ! solve eq. 18
    do iatom = 1, npos
      if (modine%Jlocal(iatom) > m_one-1e-5_real64) then
        modine%Jrange(iatom) = kappa(iatom)
      else
        modine%Jrange(iatom) = kappa(iatom) / &
          sqrt(m_one - product_log(-m_half*(((m_one+modine%Jlocal(1))/m_two)**m_third/(m_one-modine%Jlocal(iatom)**m_third)-m_one)&
          *sqrt(m_e)))
      end if
    end do
 
    ! initialize root solver for the inversion of the coordinates
    call root_solver_init(modine%rs, namespace, dim, &
      solver_type = root_newton, maxiter = 500, abs_tolerance = 1.0e-13_real64)
 
    safe_allocate(modine%xi(1:modine%dim, 1:modine%npos))
    safe_allocate(modine%Q(1:modine%dim, 1:modine%dim, 1:modine%npos))
    call optimize()
 
    modine%min_mesh_scaling_product = modine%Jbar**modine%dim
 
    pop_sub(curv_modine_constructor)
  contains
 
    ! product log function, i.e., the inverse of x exp(x)
    real(real64) pure function product_log(x)
      real(real64), intent(in) :: x
 
      ! no push_sub/pop_sub, pure function
 
      ! Third-order Taylor expansion around 0
      product_log = x * (m_one + x * (-m_one + m_three*m_half*x))
    end function product_log
 
    subroutine optimize()
      integer :: iatom, idim
      real(real64), allocatable :: xi_tmp(:, :), q_tmp(:, :, :)
      real(real64) :: residue, rr2, dd
      integer :: iter, nu, mu, idim1, idim2
      integer, parameter :: maxiter = 1000
      real(real64), parameter :: tolerance = 1.0e-10_real64
 
      push_sub(curv_modine_constructor.optimize)
 
      safe_allocate(xi_tmp(1:modine%dim, 1:modine%npos))
      safe_allocate(q_tmp(1:modine%dim, 1:modine%dim, 1:modine%npos))
 
      ! initial guess: xi = R, Q = delta (1-J**1/3)
      modine%Q(:, :, :) = m_zero
      do iatom = 1, modine%npos
        modine%xi(:, iatom) = pos(:, iatom)
        do idim = 1, modine%dim
          modine%Q(idim, idim, iatom) = m_one-modine%Jlocal(iatom)**(m_one/real(modine%dim, real64))
        end do
      end do
 
      ! relaxation method to solve eq. 16 and 17
      do iter = 1, maxiter
        ! copy current values to temporary arrays
        xi_tmp(:, :) = modine%xi(:, :)
        q_tmp(:, :, :) = modine%Q(:, :, :)
 
        modine%Q(:, :, :) = m_zero
        do mu = 1, modine%npos
          ! right-hand side
          modine%xi(:, mu) = pos(:, mu)
          do idim = 1, modine%dim
            ! right-hand side and diagonal part
            modine%Q(idim, idim, mu) = m_one-modine%Jlocal(mu)**(m_one/real(modine%dim, real64))
          end do
          do nu = 1, modine%npos
            if (mu == nu) cycle
            rr2 = sum((xi_tmp(:, mu) - xi_tmp(:, nu))**2)
            dd = exp(-rr2/(m_two*modine%Jrange(nu)**2))
 
            modine%xi(:, mu) = modine%xi(:, mu) + &
              matmul(modine%Q(:, :, mu), xi_tmp(:, mu) - xi_tmp(:, nu)) * dd
 
            do idim1 = 1, modine%dim
              do idim2 = 1, modine%dim
                modine%Q(idim1, idim2, mu) = modine%Q(idim1, idim2, mu) - q_tmp(idim1, idim2, nu) * dd
                modine%Q(idim1, idim2, mu) = modine%Q(idim1, idim2, mu) + &
                  sum(q_tmp(idim1, :, nu)*(xi_tmp(:, mu) - xi_tmp(:, nu)))*(xi_tmp(idim2, mu) - xi_tmp(idim2, nu)) * &
                  m_one/modine%Jrange(nu) * dd
              end do
            end do
          end do
        end do
 
        ! compute residue
        residue = norm2(xi_tmp(:, :) - modine%xi(:, :)) / real(modine%dim*modine%npos, real64)
        residue = residue + norm2(q_tmp(:, :, :) - modine%Q(:, :, :)) / real(modine%dim*modine%dim*modine%npos, real64)
        if (residue <= tolerance) exit
      end do
 
      ! check for convergence
      if (residue > tolerance) then
        message(1) = "During the construction of the adaptive grid, the relaxation"
        message(2) = "method did not converge."
        message(3) = "This might be due to values of CurvModineJrange that are too large"
        message(4) = "or CurvModineJlocal that are too small."
        call messages_fatal(4, namespace=namespace)
      end if
 
 
      safe_deallocate_a(xi_tmp)
      safe_deallocate_a(q_tmp)
 
      pop_sub(curv_modine_constructor.optimize)
    end subroutine optimize
 
  end function curv_modine_constructor
 
  ! ---------------------------------------------------------
  subroutine curv_modine_copy(this_out, this_in)
    type(curv_modine_t), intent(inout) :: this_out
    type(curv_modine_t), intent(in)    :: this_in
 
    push_sub(curv_modine_copy)
 
    this_out%dim = this_in%dim
    this_out%local_basis = this_in%local_basis
    safe_allocate_source_a(this_out%lsize, this_in%lsize)
    this_out%xbar = this_in%xbar
    this_out%Jbar = this_in%Jbar
    safe_allocate_source_a(this_out%Jlocal, this_in%Jlocal)
    safe_allocate_source_a(this_out%Jrange, this_in%Jrange)
    safe_allocate_source_a(this_out%xi, this_in%xi)
    safe_allocate_source_a(this_out%Q, this_in%Q)
    this_out%npos = this_in%npos
 
    pop_sub(curv_modine_copy)
  end subroutine curv_modine_copy
 
  ! ---------------------------------------------------------
  subroutine curv_modine_finalize(this)
    type(curv_modine_t), intent(inout) :: this
 
    push_sub(curv_modine_finalize)
 
    safe_deallocate_a(this%lsize)
    safe_deallocate_a(this%Jlocal)
    safe_deallocate_a(this%Jrange)
    safe_deallocate_a(this%xi)
    safe_deallocate_a(this%Q)
 
    pop_sub(curv_modine_finalize)
  end subroutine curv_modine_finalize
 
  ! ---------------------------------------------------------
  function curv_modine_to_cartesian(this, chi) result(xx)
    class(curv_modine_t), target, intent(in)  :: this
    real(real64),                 intent(in)  :: chi(:)
    real(real64) :: xx(1:this%dim)
 
    real(real64) :: chi2(this%dim), rr2, dd
    integer :: iatom
 
    ! no PUSH_SUB, called too often
 
    ! eq. 14
    call curv_modine_chi2chi2(this, chi, chi2)
 
    ! eq. 15
    xx(:) = chi2(:)
    do iatom = 1, this%npos
      rr2 = sum((chi2(:) - this%xi(:, iatom))**2)
      dd = exp(-rr2/(m_two*this%Jrange(iatom)**2))
 
      xx(:) = xx(:) - matmul(this%Q(:, :, iatom), chi2(:) - this%xi(:, iatom)) * dd
    end do
 
  end function curv_modine_to_cartesian
 
  ! ---------------------------------------------------------
  function curv_modine_from_cartesian(this, xx) result(chi)
    class(curv_modine_t), target, intent(in)  :: this
    real(real64),                 intent(in)  :: xx(:)
    real(real64) :: chi(1:this%dim)
 
    logical :: conv
    real(real64) :: start(1:this%dim)
    integer :: i
 
    ! no PUSH_SUB, called too often
 
    modine_p  => this
    safe_allocate(x_p(1:this%dim))
    x_p(:) = xx(:)
 
    ! get a better starting point for the global backdrop
    do i = 1, this%dim
      if (abs(xx(i)) > this%lsize(i)) then
        start(i) = xx(i)
      else if (abs(xx(i)) > this%xbar*this%lsize(i)) then
        ! linear interpolation between the two cases
        start(i) = sign(m_one, xx(i)) * (abs(xx(i))/this%Jbar + abs(xx(i)) *(m_one - m_one/this%Jbar) * &
          (abs(xx(i)) - this%xbar*this%lsize(i)) / (this%lsize(i) - this%xbar*this%lsize(i)))
        ! ensure that we stay in the bounds
        if (abs(start(i)) >= this%lsize(i)) then
          start(i) = sign(this%lsize(i), xx(i))
        end if
      else
        start(i) = xx(i)/this%Jbar
      end if
    end do
 
    call droot_solver_run(this%rs, getf, chi, conv, startval=start)
 
    safe_deallocate_a(x_p)
    nullify(modine_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),'(3f14.6)') xx(1:this%dim)
      call messages_fatal(3)
    end if
 
  end function curv_modine_from_cartesian
 
  ! ---------------------------------------------------------
  subroutine curv_modine_write_info(this, iunit, namespace)
    class(curv_modine_t),         intent(in) :: this
    integer,            optional, intent(in) :: iunit
    type(namespace_t),  optional, intent(in) :: namespace
 
    push_sub(curv_modine_write_info)
 
    write(message(1), '(a)') '  Curvilinear Method = modine'
    call messages_info(1, iunit=iunit, namespace=namespace)
 
    pop_sub(curv_modine_write_info)
  end subroutine curv_modine_write_info
 
  ! ---------------------------------------------------------
  real(real64) function curv_modine_surface_element(this, idir) result(ds)
    class(curv_modine_t), intent(in) :: this
    integer,              intent(in) :: idir
 
    ds = m_zero
    message(1) = 'Surface element with modine curvilinear coordinates not implemented'
    call messages_fatal(1)
 
  end function curv_modine_surface_element
 
  ! ---------------------------------------------------------
  pure subroutine curv_modine_chi2chi2(this, chi_, chi2, Jac)
    class(curv_modine_t), intent(in) :: this
    real(real64), intent(in) :: chi_(:)
    real(real64), contiguous, intent(out) :: chi2(:)
    real(real64), contiguous, optional, intent(out) :: jac(:)
 
    integer, parameter :: qq = 3
    real(real64) :: chibar(this%dim), rr, chi
    logical :: neg
    integer :: i
 
    ! no PUSH_SUB, called too often
 
    chibar = this%xbar*this%lsize/this%Jbar
 
    do i = 1, this%dim
      neg = (chi_(i) < 0)
      chi = abs(chi_(i))
 
      chi2(i)  = this%Jbar * chi
      if (present(jac)) jac(i) = this%Jbar
 
      if (chi > this%lsize(i)) then
        ! needed for boundary points
        chi2(i) = chi
        if (present(jac)) jac(i) = m_one
      else if (chi > chibar(i)) then
        rr = (chi - chibar(i))/(this%lsize(i) - chibar(i))
 
        chi2(i)  = chi2(i) + this%lsize(i)/m_two*(1 - this%Jbar) * rr**qq * &
          (qq + m_one - (qq - m_one)*rr)
 
        if (present(jac)) then
          jac(i) = jac(i) + this%lsize(i)/m_two*(m_one - this%Jbar) * rr**(qq - 1)/(this%lsize(i) - chibar(i)) * &
            (qq*(qq + m_one) - (qq**2 - m_one)*rr)
        end if
      end if
 
      if (neg) chi2(i) = -chi2(i)
      ! CHECK if Jacobian does not have to be negated!
    end do
 
  end subroutine curv_modine_chi2chi2
 
  ! ---------------------------------------------------------
  ! Jacobian for the coordinate transformation, i.e. for eqs. 14 and 15
  pure subroutine curv_modine_jacobian_inv(this, chi, xx, Jac)
    type(curv_modine_t), intent(in)  :: this
    real(real64),        intent(in)  :: chi(:)
    real(real64),        intent(out) :: xx(:)
    real(real64),        intent(out) :: jac(:, :)
 
    real(real64) :: chi2(this%dim), rr2, dd, j2(this%dim)
    integer :: iatom, idim, idim2
 
    ! no PUSH_SUB, called too often
 
    call curv_modine_chi2chi2(this, chi, chi2, j2)
 
    ! initialize both xx and the Jacobian
    xx(:) = chi2(:)
    jac(:,:) = m_zero
    do idim = 1, this%dim
      jac(idim, idim) = m_one
    end do
 
    do iatom = 1, this%npos
      rr2 = sum((chi2(:) - this%xi(:, iatom))**2)
      dd = exp(-rr2/(m_two*this%Jrange(iatom)**2))
 
      xx(:) = xx(:) - matmul(this%Q(:, :, iatom), chi2(:) - this%xi(:, iatom)) * dd
 
      do idim = 1, this%dim
        do idim2 = 1, this%dim
          jac(idim, idim2) = jac(idim, idim2) - this%Q(idim, idim2, iatom) * dd
          jac(idim, idim2) = jac(idim, idim2) + &
            sum(this%Q(idim, :, iatom)*(chi2(:) - this%xi(:, iatom)))*(chi2(idim2) - this%xi(idim2, iatom)) * &
            m_one/this%Jrange(iatom) * dd
        end do
      end do
    end do
 
    do idim = 1, this%dim
      jac(idim, :) = jac(idim, :) * j2(:)
    end do
 
  end subroutine curv_modine_jacobian_inv
 
  ! ---------------------------------------------------------
  pure subroutine getf(yy, ff, jf)
    real(real64), intent(in)  :: yy(:)
    real(real64), intent(out) :: ff(:), jf(:, :)
 
    ! no PUSH_SUB, called too often
 
    call curv_modine_jacobian_inv(modine_p, yy, ff, jf)
    ff(:) = ff(:) - x_p(:)
 
  end subroutine getf
 
  ! ---------------------------------------------------------
  function curv_modine_jacobian(this, chi) result(jacobian)
    class(curv_modine_t), intent(in)  :: this
    real(real64),         intent(in)  :: chi(:)
    real(real64) :: jacobian(1:this%dim, 1:this%dim)
 
    real(real64) :: dummy(this%dim)
 
    call curv_modine_jacobian_inv(this, chi, dummy, jacobian)
 
  end function curv_modine_jacobian
 
  ! ---------------------------------------------------------
  function curv_modine_jacobian_inverse(this, chi) result(jacobian_inverse)
    class(curv_modine_t), intent(in)  :: this
    real(real64),         intent(in)  :: chi(:)
    real(real64) :: jacobian_inverse(1:this%dim, 1:this%dim)
 
    real(real64) :: dummy(this%dim)
 
    call curv_modine_jacobian_inv(this, chi, dummy, jacobian_inverse)
    call lalg_inverse(this%dim, jacobian_inverse, "dir")
 
  end function curv_modine_jacobian_inverse
 
end module curv_modine_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: