TABLE OF CONTENTS


ABINIT/getcut [ Functions ]

[ Top ] [ 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