TABLE OF CONTENTS


ABINIT/invcb [ Functions ]

[ Top ] [ Functions ]

NAME

 invcb

FUNCTION

 Compute a set of inverse cubic roots as fast as possible :
 rspts(:)=rhoarr(:)$^\frac{-1}{3}$

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  npts=number of real space points on which density is provided
  rhoarr(npts)=input data

OUTPUT

  rspts(npts)=inverse cubic root of rhoarr

PARENTS

      drivexc,gammapositron,xchcth,xcpbe,xcpositron,xctfw

CHILDREN

SOURCE

 31 #if defined HAVE_CONFIG_H
 32 #include "config.h"
 33 #endif
 34 
 35 #include "abi_common.h"
 36 
 37 
 38 subroutine invcb(rhoarr,rspts,npts)
 39 
 40  use defs_basis
 41  use m_profiling_abi
 42 
 43 !This section has been created automatically by the script Abilint (TD).
 44 !Do not modify the following lines by hand.
 45 #undef ABI_FUNC
 46 #define ABI_FUNC 'invcb'
 47 !End of the abilint section
 48 
 49  implicit none
 50 
 51 !Arguments ------------------------------------
 52 !scalars
 53  integer,intent(in) :: npts
 54 !arrays
 55  real(dp),intent(in) :: rhoarr(npts)
 56  real(dp),intent(out) :: rspts(npts)
 57 
 58 !Local variables-------------------------------
 59 !scalars
 60  integer :: ii,ipts
 61  real(dp),parameter :: c2_27=2.0e0_dp/27.0e0_dp,c5_9=5.0e0_dp/9.0e0_dp
 62  real(dp),parameter :: c8_9=8.0e0_dp/9.0e0_dp,m1thrd=-third
 63  real(dp) :: del,prod,rho,rhom1,rhomtrd
 64  logical :: test
 65 !character(len=500) :: message
 66 
 67 ! *************************************************************************
 68 
 69 !Loop over points : here, brute force algorithm
 70 !do ipts=1,npts
 71 !rspts(ipts)=sign( (abs(rhoarr(ipts)))**m1thrd,rhoarr(ipts))
 72 !end do
 73 !
 74 
 75 !write(std_out,*)' invcb : rhoarr, rspts'
 76 
 77  rhomtrd=sign( (abs(rhoarr(1)))**m1thrd, rhoarr(1) )
 78  rhom1=one/rhoarr(1)
 79  rspts(1)=rhomtrd
 80  do ipts=2,npts
 81 !  write(std_out,*)
 82 !  write(std_out,*)rhoarr(ipts),rspts(ipts)
 83    rho=rhoarr(ipts)
 84    prod=rho*rhom1
 85 !  If the previous point is too far ...
 86    if(prod < 0.01_dp .or. prod > 10._dp )then
 87      rhomtrd=sign( (abs(rho))**m1thrd , rho )
 88      rhom1=one/rho
 89    else
 90      del=prod-one
 91      do ii=1,5
 92 !      Choose one of the two next lines, the last one is more accurate
 93 !      rhomtrd=((one+third*del)/(one+two_thirds*del))*rhomtrd
 94        rhomtrd=((one+c5_9*del)/(one+del*(c8_9+c2_27*del)))*rhomtrd
 95        rhom1=rhomtrd*rhomtrd*rhomtrd
 96        del=rho*rhom1-one
 97 !      write(std_out,*)rhomtrd,del
 98        test = del*del < 1.0e-24_dp
 99        if(test) exit
100      end do
101      if( .not. test) then
102        rhomtrd=sign( (abs(rho))**m1thrd , rho )
103      end if
104    end if
105    rspts(ipts)=rhomtrd
106  end do
107 
108 end subroutine invcb