TABLE OF CONTENTS


ABINIT/mkkin [ Functions ]

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