TABLE OF CONTENTS
ABINIT/getcut [ Functions ]
NAME
getcut
FUNCTION
For input kpt, fft box dim ngfft(1:3), recip space metric gmet, and kinetic energy cutoff ecut (hartree), COMPUTES: if iboxcut==0: gsqcut : cut-off on G^2 for "large sphere" of radius double that of the basis sphere corresponding to ecut boxcut : where boxcut == gcut(box)/gcut(sphere). boxcut >=2 for no aliasing. boxcut < 1 is wrong and halts subroutine. if iboxcut==1: gsqcut : cut-off on G^2 for "large sphere" containing the whole fft box boxcut : no meaning (zero)
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (DCA, XG, GMR, MT) 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
ecut=kinetic energy cutoff for planewave sphere (hartree) gmet(3,3)=reciprocal space metric (bohr^-2) iboxcut=0: compute gsqcut and boxcut with boxcut>=1 1: compute gsqcut for boxcut=1 (sphere_cutoff=box_cutoff) iout=unit number for output file kpt(3)=input k vector (reduced coordinates--in terms of reciprocal lattice primitive translations) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
OUTPUT
boxcut=defined above (dimensionless), ratio of basis sphere diameter to fft box length (smallest value) gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere--appropriate for charge density rho(G), Hartree potential, and pseudopotentials
NOTES
2*gcut arises from rho(g)=sum g prime (psi(g primt)*psi(g prime+g)) where psi(g) is only nonzero for |g| <= gcut). ecut (currently in hartree) is proportional to gcut(sphere)**2.
PARENTS
dfpt_looppert,dfpt_scfcv,m_pawfgr,nonlinear,respfn,scfcv,setup1 setup_positron
CHILDREN
bound,wrtout
SOURCE
58 #if defined HAVE_CONFIG_H 59 #include "config.h" 60 #endif 61 62 #include "abi_common.h" 63 64 65 subroutine getcut(boxcut,ecut,gmet,gsqcut,iboxcut,iout,kpt,ngfft) 66 67 use defs_basis 68 use m_profiling_abi 69 use m_errors 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'getcut' 75 use interfaces_14_hidewrite 76 use interfaces_52_fft_mpi_noabirule 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 integer,intent(in) :: iboxcut,iout 84 real(dp),intent(in) :: ecut 85 real(dp),intent(out) :: boxcut,gsqcut 86 !arrays 87 integer,intent(in) :: ngfft(18) 88 real(dp),intent(in) :: gmet(3,3),kpt(3) 89 90 !Local variables------------------------------- 91 !scalars 92 integer :: plane 93 real(dp) :: boxsq,cutrad,ecut_pw,effcut,largesq,sphsq 94 character(len=500) :: message 95 !arrays 96 integer :: gbound(3) 97 98 ! ************************************************************************* 99 100 !This is to treat the case where ecut has not been initialized, 101 !for wavelet computations. The default for ecut is -1.0 , allowed 102 !only for wavelets calculations 103 ecut_pw=ecut 104 if(ecut<-tol8)ecut_pw=ten 105 106 !gcut(box)**2=boxsq; gcut(sphere)**2=sphsq 107 !get min. d**2 to boundary of fft box: 108 !(gmet sets dimensions: bohr**-2) 109 !ecut(sphere)=0.5*(2 pi)**2 * sphsq: 110 call bound(largesq,boxsq,gbound,gmet,kpt,ngfft,plane) 111 effcut=0.5_dp * (two_pi)**2 * boxsq 112 sphsq=2._dp*ecut_pw/two_pi**2 113 114 if (iboxcut/=0) then 115 boxcut=10._dp 116 gsqcut=(largesq/sphsq)*(2.0_dp*ecut)/two_pi**2 117 118 write(message, '(a,a,3f8.4,a,3i4,a,a,f11.3,a,a)' ) ch10,& 119 & ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,& 120 & ' ecut(hartree)=',ecut_pw+tol8,ch10,'=> whole FFT box selected' 121 if(iout/=std_out) then 122 call wrtout(iout,message,'COLL') 123 end if 124 call wrtout(std_out,message,'COLL') 125 126 else 127 128 ! Get G^2 cutoff for sphere of double radius of basis sphere 129 ! for selecting G s for rho(G), V_Hartree(G), and V_psp(G)-- 130 ! cut off at fft box boundary or double basis sphere radius, whichever 131 ! is smaller. If boxcut were 2, then relation would be 132 !$ecut_eff = (1/2) * (2 Pi Gsmall)^2 and gsqcut=4*Gsmall^2$. 133 boxcut = sqrt(boxsq/sphsq) 134 cutrad = min(2.0_dp,boxcut) 135 gsqcut = (cutrad**2)*(2.0_dp*ecut_pw)/two_pi**2 136 137 if(ecut>-tol8)then 138 139 write(message, '(a,a,3f8.4,a,3i4,a,a,f11.3,3x,a,f10.5)' ) ch10,& 140 & ' getcut: wavevector=',kpt,' ngfft=',ngfft(1:3),ch10,& 141 & ' ecut(hartree)=',ecut+tol8,'=> boxcut(ratio)=',boxcut+tol8 142 if(iout/=std_out) then 143 call wrtout(iout,message,'COLL') 144 end if 145 call wrtout(std_out,message,'COLL') 146 147 if (boxcut<1.0_dp) then 148 write(message, '(a,a,a,a,a,a,a,a,a,f12.6)' )& 149 & ' Choice of acell, ngfft, and ecut',ch10,& 150 & ' ===> basis sphere extends BEYOND fft box !',ch10,& 151 & ' Recall that boxcut=Gcut(box)/Gcut(sphere) must be > 1.',ch10,& 152 & ' Action : try larger ngfft or smaller ecut.',ch10,& 153 & ' Note that ecut=effcut/boxcut**2 and effcut=',effcut+tol8 154 if(iout/=std_out) then 155 call wrtout(iout,message,'COLL') 156 end if 157 MSG_ERROR(message) 158 end if 159 160 if (boxcut>2.2_dp) then 161 write(message, '(a,a,a,a,a,a,a,a,a,a,a,f12.6,a,a)' ) ch10,& 162 & ' getcut : COMMENT -',ch10,& 163 & ' Note that boxcut > 2.2 ; recall that',' boxcut=Gcut(box)/Gcut(sphere) = 2',ch10,& 164 & ' is sufficient for exact treatment of convolution.',ch10,& 165 & ' Such a large boxcut is a waste : you could raise ecut',ch10,& 166 & ' e.g. ecut=',effcut*0.25_dp+tol8,' Hartrees makes boxcut=2',ch10 167 if(iout/=std_out) then 168 call wrtout(iout,message,'COLL') 169 end if 170 call wrtout(std_out,message,'COLL') 171 end if 172 173 if (boxcut<1.5_dp) then 174 write(message, '(a,a,a,a,a,a,a,a,a,a)' ) ch10,& 175 & ' getcut : WARNING -',ch10,& 176 & ' Note that boxcut < 1.5; this usually means',ch10,& 177 & ' that the forces are being fairly strongly affected by',' the smallness of the fft box.',ch10,& 178 & ' Be sure to test with larger ngfft(1:3) values.',ch10 179 if(iout/=std_out) then 180 call wrtout(iout,message,'COLL') 181 end if 182 call wrtout(std_out,message,'COLL') 183 end if 184 185 end if 186 187 end if ! iboxcut 188 189 end subroutine getcut