./56_recipspace/kpgstr.F90

ABINIT/kpgstr

ABINIT/kpgstr [ Functions ]

NAME

 kpgstr

FUNCTION

 Compute elements of the derivative the kinetic energy operator in reciprocal
 space at given k point wrt a single cartesian strain component

 Copyright (C) 1999-2017 ABINIT group (DRH, 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 .

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 (Bohr**-2)
  gprimd(3,3)=reciprocal space dimensional primitive translations
  istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21
  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

  dkinpw(npw)=d/deps(istr) ( (1/2)*(2 pi)**2 * (k+G)**2 )

NOTES

  Src_6response/kpg3.f

PARENTS

      dfpt_nsteltwf,dfpt_nstpaw,dfpt_rhofermi,getgh1c

CHILDREN

SOURCE

 40 #if defined HAVE_CONFIG_H
 41 #include "config.h"
 42 #endif
 43 
 44 #include "abi_common.h"
 45 
 46 subroutine kpgstr(dkinpw,ecut,ecutsm,effmass_free,gmet,gprimd,istr,kg,kpt,npw)
 47 
 48  use defs_basis
 49  use m_profiling_abi
 50  use m_errors
 51 
 52 !This section has been created automatically by the script Abilint (TD).
 53 !Do not modify the following lines by hand.
 54 #undef ABI_FUNC
 55 #define ABI_FUNC 'kpgstr'
 56 !End of the abilint section
 57 
 58  implicit none
 59 
 60 !Arguments -------------------------------
 61 !scalars
 62  integer,intent(in) :: istr,npw
 63  real(dp),intent(in) :: ecut,ecutsm,effmass_free
 64 !arrays
 65  integer,intent(in) :: kg(3,npw)
 66  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt(3)
 67  real(dp),intent(out) :: dkinpw(npw)
 68 
 69 !Local variables -------------------------
 70 !scalars
 71  integer :: ig,ii,ka,kb
 72  real(dp) :: dfsm,dkinetic,dkpg2,ecutsm_inv,fsm,gpk1,gpk2,gpk3,htpisq
 73 ! real(dp) :: d2fsm ! used in commented section below
 74  real(dp) :: kpg2,xx
 75  character(len=500) :: message
 76 !arrays
 77  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 78  real(dp) :: dgmetds(3,3)
 79 
 80 ! *********************************************************************
 81 
 82 !htpisq is (1/2) (2 Pi) **2:
 83  htpisq=0.5_dp*(two_pi)**2
 84 
 85  ecutsm_inv=0.0_dp
 86  if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
 87 
 88 !Compute derivative of metric tensor wrt strain component istr
 89  if(istr<1 .or. istr>6)then
 90    write(message, '(a,i10,a,a,a)' )&
 91 &   'Input istr=',istr,' not allowed.',ch10,&
 92 &   'Possible values are 1,2,3,4,5,6 only.'
 93    MSG_BUG(message)
 94  end if
 95 
 96  ka=idx(2*istr-1);kb=idx(2*istr)
 97  do ii = 1,3
 98    dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
 99  end do
100 !For historical reasons:
101  dgmetds(:,:)=0.5_dp*dgmetds(:,:)
102 
103  do ig=1,npw
104    gpk1=dble(kg(1,ig))+kpt(1)
105    gpk2=dble(kg(2,ig))+kpt(2)
106    gpk3=dble(kg(3,ig))+kpt(3)
107    kpg2=htpisq*&
108 &   ( gmet(1,1)*gpk1**2+         &
109 &   gmet(2,2)*gpk2**2+         &
110 &   gmet(3,3)*gpk3**2          &
111 &   +2.0_dp*(gpk1*gmet(1,2)*gpk2+  &
112 &   gpk1*gmet(1,3)*gpk3+  &
113 &   gpk2*gmet(2,3)*gpk3 )  )
114    dkpg2=htpisq*2.0_dp*&
115 &   (gpk1*(dgmetds(1,1)*gpk1+dgmetds(1,2)*gpk2+dgmetds(1,3)*gpk3)+  &
116 &   gpk2*(dgmetds(2,1)*gpk1+dgmetds(2,2)*gpk2+dgmetds(2,3)*gpk3)+  &
117 &   gpk3*(dgmetds(3,1)*gpk1+dgmetds(3,2)*gpk2+dgmetds(3,3)*gpk3) )
118    dkinetic=dkpg2
119    if(kpg2>ecut-ecutsm)then
120      if(kpg2>ecut-tol12)then
121 !      The wavefunction has been filtered : no derivative
122        dkinetic=0.0_dp
123      else
124        xx=(ecut-kpg2)*ecutsm_inv
125 !      This kinetic cutoff smoothing function and its xx derivatives
126 !      were produced with Mathematica and the fortran code has been
127 !      numerically checked against Mathematica.
128        fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
129        dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
130 !      d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
131 !      &                        (-144+45*xx))))))*fsm**3
132        dkinetic=dkpg2*(fsm-ecutsm_inv*kpg2*dfsm)
133      end if
134    end if
135    dkinpw(ig)=dkinetic/effmass_free
136  end do
137 
138 end subroutine kpgstr

TABLE OF CONTENTS

ABINIT/kpgstr [ Functions ]