logrid.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 logrid_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::         &
    logrid_t,       &
    logrid_init,    &
    logrid_end,     &
    logrid_copy,    &
    logrid_index,   &
    logrid_find_parameters
 
  integer, parameter, public :: &
    LOGRID_PSF  = 1, & !< log grid used in Troullier-Martins code
    logrid_cpi  = 2    
 
  type logrid_t
    ! Components are public by default
    integer  :: flavor
 
    real(real64)    :: a, b
    integer  :: nrval
 
    real(real64), allocatable :: rofi(:)
    real(real64), allocatable :: r2ofi(:)
    real(real64), allocatable :: drdi(:)
    real(real64), allocatable :: s(:)
  end type logrid_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine logrid_init(grid, flavor, aa, bb, nrval)
    type(logrid_t), intent(out) :: grid
    integer,        intent(in)  :: flavor
    real(real64),   intent(in)  :: aa, bb
    integer,        intent(in)  :: nrval
 
    real(real64) :: rpb, ea
    integer  :: ir
 
    push_sub(logrid_init)
 
    assert(flavor == logrid_psf .or. flavor == logrid_cpi)
    assert(nrval >= 2)
    assert(aa > m_zero)
    assert(bb > m_zero)
 
    grid%flavor = flavor
    grid%a = aa
    grid%b = bb
    grid%nrval = nrval
 
    safe_allocate(grid%rofi(1:nrval))
    safe_allocate(grid%r2ofi(1:nrval))
    safe_allocate(grid%drdi(1:nrval))
    safe_allocate(grid%s(1:nrval))
 
    select case (grid%flavor)
    case (logrid_psf)
      rpb = bb
      ea  = exp(aa)
      do ir = 1, nrval
        grid%drdi(ir) = aa*rpb
        rpb           = rpb*ea
        grid%rofi(ir) = bb*expm1(aa*(ir-1))
      end do
 
    case (logrid_cpi)
      grid%rofi(1) = m_zero
      grid%drdi(1) = m_zero
 
      rpb = log(aa)
      grid%rofi(2) = bb
      grid%drdi(2) = bb*rpb
      do ir = 3, grid%nrval
        grid%rofi(ir) = grid%rofi(ir-1)*aa
        grid%drdi(ir) = grid%rofi(ir)*rpb
      end do
    end select
 
    ! calculate sqrt(drdi)
    do ir = 1, grid%nrval
      grid%s(ir)     = sqrt(grid%drdi(ir))
      grid%r2ofi(ir) = grid%rofi(ir)**2
    end do
 
    pop_sub(logrid_init)
  end subroutine logrid_init
 
 
  ! ---------------------------------------------------------
  subroutine logrid_end(grid)
    type(logrid_t), intent(inout) :: grid
 
    push_sub(logrid_end)
 
    safe_deallocate_a(grid%rofi)
    safe_deallocate_a(grid%r2ofi)
    safe_deallocate_a(grid%drdi)
    safe_deallocate_a(grid%s)
 
    pop_sub(logrid_end)
  end subroutine logrid_end
 
 
  ! ---------------------------------------------------------
  subroutine logrid_copy(grid_in, grid_out)
    type(logrid_t), intent(in)    :: grid_in
    type(logrid_t), intent(inout) :: grid_out
 
    push_sub(logrid_copy)
 
    call logrid_end(grid_out)
 
    grid_out%flavor = grid_in%flavor
    grid_out%a      = grid_in%a
    grid_out%b      = grid_in%b
    grid_out%nrval  = grid_in%nrval
 
    safe_allocate(grid_out%rofi (1:grid_out%nrval))
    safe_allocate(grid_out%r2ofi(1:grid_out%nrval))
    safe_allocate(grid_out%drdi (1:grid_out%nrval))
    safe_allocate(grid_out%s    (1:grid_out%nrval))
 
    grid_out%rofi(:)  = grid_in%rofi(:)
    grid_out%r2ofi(:) = grid_in%r2ofi(:)
    grid_out%drdi(:)  = grid_in%drdi(:)
    grid_out%s(:)     = grid_in%s(:)
 
    pop_sub(logrid_copy)
  end subroutine logrid_copy
 
 
  ! ---------------------------------------------------------
  integer function logrid_index(grid, rofi) result(ii)
    type(logrid_t), intent(in) :: grid
    real(real64),   intent(in) :: rofi
 
    integer :: ir
 
    push_sub(logrid_index)
 
    ii = 0
    do ir = 1, grid%nrval-1
 
      if (rofi >= grid%rofi(ir).and.rofi < grid%rofi(ir+1)) then
        if (abs(rofi-grid%rofi(ir)) < abs(rofi-grid%rofi(ir+1))) then
          ii = ir
        else
          ii = ir + 1
        end if
        exit
      end if
 
    end do
 
    pop_sub(logrid_index)
  end function logrid_index
 
  ! ----------------------------------------------------------
  subroutine logrid_find_parameters(namespace, zz, aa, bb, np)
    type(namespace_t),  intent(in)  :: namespace
    integer,            intent(in)  :: zz
    real(real64),       intent(out) :: aa, bb
    integer,            intent(out) :: np
 
    real(real64) :: xmin, xmax, a1, a2, f1, fm
 
    push_sub(logrid_find_parameters)
 
    ! Initializes the logarithmic grid.
    ! Parameters are obtained using the default values for the first non-zero point xmin,
    ! the last point xmax, and the number of points np
    ! These values have a default value obtained from the atomic number
    ! Adapted from APE
    xmin = sqrt(real(zz, real64) )*1e-5_real64
    xmax = sqrt(real(zz, real64) )*30.0_real64
    np = floor(sqrt(real(zz, real64) )*200_real64)
    ! The code wants np to be an odd number
    np = floor(np/m_two)*2+1
 
    a1 = 1e-8_real64
    f1 = func(xmin, xmax, real(np, real64) , a1)
    a2 = m_one
    do
      aa = (a2 + a1)*m_half
      fm = func(xmin, xmax, real(np, real64) , aa)
      if (m_half*abs(a1 - a2) < 1.0e-16_real64) exit
      if (fm*f1 > m_zero) then
        a1 = aa
        f1 = fm
      else
        a2 = aa
      end if
    end do
 
    bb = xmin/expm1(aa)
 
    write(message(1), '(a,es13.6,a,es13.6,a,i4)') 'Debug: Log grid parameters: a = ', aa, &
      ' b = ', bb, ' np = ', np
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(logrid_find_parameters)
  contains
    real(real64) function func(r1, rn, n, a)
      real(real64), intent(in) :: r1, rn, a, n
      if((n-m_one)*a < m_max_exp_arg) then ! To avoid FPE
        func = r1*expm1((n-m_one)*a) - rn*expm1(a)
      else
        func = m_huge*r1 - m_one*r1 - rn*expm1(a)
      end if
    end function func
  end subroutine logrid_find_parameters
 
end module logrid_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: