TABLE OF CONTENTS
ABINIT/invcb [ 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