TABLE OF CONTENTS
ABINIT/mkkin [ Functions ]
NAME
mkkin
FUNCTION
compute elements of kinetic energy operator in reciprocal space at a given k point
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (DCA, XG, GMR, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$) idir1 = 1st direction of the derivative (if 1 <= idir1 <= 3, not used otherwise) idir2 = 2st direction of the derivative (if 1 <= idir1,idir2 <= 3, not used otherwise)) kg(3,npw)=integer coordinates of planewaves in basis sphere. kpt(3)=reduced coordinates of k point npw=number of plane waves at kpt.
OUTPUT
kinpw(npw)=(modified) kinetic energy (or derivative) for each plane wave (Hartree)
SIDE EFFECTS
NOTES
Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $ However, the present implementation allows for a modification of this kinetic energy, in order to obtain smooth total energy curves with respect to the cut-off energy or the cell size and shape. Thus the usual expression is kept if it is lower then ecut-ecutsm, zero is returned beyond ecut, and in between, the kinetic energy is DIVIDED by a smearing factor (to make it infinite at the cut-off energy). The smearing factor is $x^2 (3-2x)$, where x = (ecut- unmodified energy)/ecutsm. This smearing factor is also used to derived a modified kinetic contribution to stress, in another routine (forstrnps.f) Also, in order to break slightly the symmetry between axes, that causes sometimes a degeneracy of eigenvalues and do not allow to obtain the same results on different machines, there is a modification by one part over 1.0e12 of the metric tensor elements (1,1) and (3,3)
PARENTS
calc_vhxc_me,d2frnl,dfpt_nsteltwf,dfpt_nstpaw,dfpt_rhofermi,energy getgh1c,ks_ddiago,m_io_kss,m_vkbr,mkffnl,vtorho
CHILDREN
SOURCE
58 #if defined HAVE_CONFIG_H 59 #include "config.h" 60 #endif 61 62 #include "abi_common.h" 63 64 65 subroutine mkkin (ecut,ecutsm,effmass_free,gmet,kg,kinpw,kpt,npw,idir1,idir2) 66 67 use defs_basis 68 use m_profiling_abi 69 70 !This section has been created automatically by the script Abilint (TD). 71 !Do not modify the following lines by hand. 72 #undef ABI_FUNC 73 #define ABI_FUNC 'mkkin' 74 !End of the abilint section 75 76 implicit none 77 78 !Arguments ------------------------------------ 79 !scalars 80 integer,intent(in) :: npw 81 integer,intent(in) :: idir1,idir2 82 real(dp),intent(in) :: ecut,ecutsm,effmass_free 83 84 !arrays 85 integer,intent(in) :: kg(3,npw) 86 real(dp),intent(in) :: gmet(3,3),kpt(3) 87 real(dp),intent(out) :: kinpw(npw) 88 89 !Local variables------------------------------- 90 !scalars 91 integer :: ig,order 92 real(dp),parameter :: break_symm=1.0d-11 93 real(dp) :: ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq,kinetic,kpg2,dkpg2,xx 94 real(dp) :: d1kpg2,d2kpg2,ddfsm, dfsm 95 !arrays 96 real(dp) :: gmet_break(3,3) 97 98 ! ************************************************************************* 99 ! 100 !htpisq is (1/2) (2 Pi) **2: 101 htpisq=0.5_dp*(two_pi)**2 102 103 ecutsm_inv=0.0_dp 104 if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm 105 106 gmet_break(:,:)=gmet(:,:) 107 gmet_break(1,1)=(1.0_dp+break_symm)*gmet(1,1) 108 gmet_break(3,3)=(1.0_dp-break_symm)*gmet(3,3) 109 110 order=0 ! Compute the kinetic operator 111 if (idir1>0.and.idir1<4) then 112 order=1 ! Compute the 1st derivative of the kinetic operator 113 if (idir2>0.and.idir2<4) then 114 order=2 ! Compute the 2nd derivative of the kinetic operator 115 end if 116 end if 117 118 !$OMP PARALLEL DO PRIVATE(dkpg2,d1kpg2,d2kpg2,gpk1,gpk2,gpk3,ig,kinetic,kpg2,xx,fsm,dfsm,ddfsm) & 119 !$OMP SHARED(kinpw,ecut,ecutsm,ecutsm_inv) & 120 !$OMP SHARED(gmet_break,htpisq,idir1,idir2,kg,kpt,npw) 121 do ig=1,npw 122 gpk1=dble(kg(1,ig))+kpt(1) 123 gpk2=dble(kg(2,ig))+kpt(2) 124 gpk3=dble(kg(3,ig))+kpt(3) 125 kpg2=htpisq*& 126 & ( gmet_break(1,1)*gpk1**2+ & 127 & gmet_break(2,2)*gpk2**2+ & 128 & gmet_break(3,3)*gpk3**2 & 129 & +2.0_dp*(gpk1*gmet_break(1,2)*gpk2+ & 130 & gpk1*gmet_break(1,3)*gpk3+ & 131 & gpk2*gmet_break(2,3)*gpk3 ) ) 132 select case (order) 133 case(0) 134 kinetic=kpg2 135 case(1) 136 dkpg2=htpisq*2.0_dp*& 137 & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3) 138 kinetic=dkpg2 139 case(2) 140 dkpg2=htpisq*2.0_dp*gmet_break(idir1,idir2) 141 kinetic=dkpg2 142 end select 143 if(kpg2>ecut-ecutsm)then 144 if(kpg2>ecut-tol12)then 145 if(order==0) then 146 ! Will filter the wavefunction, based on this value, in cgwf.f, getghc.f and precon.f 147 kinetic=huge(0.0_dp)*1.d-10 148 else 149 ! The wavefunction has been filtered : no derivative 150 kinetic=0 151 end if 152 else 153 if(order==0) then 154 xx=max( (ecut-kpg2)*ecutsm_inv , 1.0d-20) 155 else 156 xx=(ecut-kpg2)*ecutsm_inv 157 end if 158 if(order==2) then 159 d1kpg2=htpisq*2.0_dp*& 160 & (gmet_break(idir1,1)*gpk1+gmet_break(idir1,2)*gpk2+gmet_break(idir1,3)*gpk3) 161 d2kpg2=htpisq*2.0_dp*& 162 & (gmet_break(idir2,1)*gpk1+gmet_break(idir2,2)*gpk2+gmet_break(idir2,3)*gpk3) 163 end if 164 ! This kinetic cutoff smoothing function and its xx derivatives 165 ! were produced with Mathematica and the fortran code has been 166 ! numerically checked against Mathematica. 167 fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx)))) 168 if(order>0) dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2 169 if(order>1) ddfsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*(-144+45*xx))))))*fsm**3 170 select case (order) 171 case(0) 172 kinetic=kpg2*fsm 173 case(1) 174 kinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm) 175 case(2) 176 kinetic=dkpg2*fsm& 177 & -2.0_dp*d1kpg2*dfsm*ecutsm_inv*d2kpg2& 178 & +kpg2*ddfsm*(ecutsm_inv**2)*d1kpg2*d2kpg2& 179 & -kpg2*dfsm*ecutsm_inv*dkpg2 180 end select 181 end if 182 end if 183 kinpw(ig)=kinetic/effmass_free 184 end do 185 !$OMP END PARALLEL DO 186 187 end subroutine mkkin