src/ggamt.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: ggamt
   8 ! !INTERFACE:
   9 subroutine ggamt(is,ia,grhomt,gupmt,gdnmt,g2upmt,g2dnmt,g3rhomt,g3upmt,g3dnmt)
  10 ! !USES:
  11 use modmain
  12 ! !INPUT/OUTPUT PARAMETERS:
  13 !   is      : species number (in,integer)
  14 !   ia      : atom number (in,integer)
  15 !   grhomt  : |grad rho| (out,real(lmmaxvr,nrmtmax))
  16 !   gupmt   : |grad rhoup| (out,real(lmmaxvr,nrmtmax))
  17 !   gdnmt   : |grad rhodn| (out,real(lmmaxvr,nrmtmax))
  18 !   g2upmt  : grad^2 rhoup (out,real(lmmaxvr,nrmtmax))
  19 !   g2dnmt  : grad^2 rhodn (out,real(lmmaxvr,nrmtmax))
  20 !   g3rhomt : (grad rho).(grad |grad rho|) (out,real(lmmaxvr,nrmtmax))
  21 !   g3upmt  : (grad rhoup).(grad |grad rhoup|) (out,real(lmmaxvr,nrmtmax))
  22 !   g3dnmt  : (grad rhodn).(grad |grad rhodn|) (out,real(lmmaxvr,nrmtmax))
  23 ! !DESCRIPTION:
  24 !   Computes $|\nabla\rho|$, $|\nabla\rho^{\uparrow}|$,
  25 !   $|\nabla\rho^{\downarrow}|$, $\nabla^2\rho^{\uparrow}$,
  26 !   $\nabla^2\rho^{\downarrow}$, $\nabla\rho\cdot(\nabla|\nabla\rho|)$,
  27 !   $\nabla\rho^{\uparrow}\cdot(\nabla|\nabla\rho^{\uparrow}|)$ and
  28 !   $\nabla\rho^{\downarrow}\cdot(\nabla|\nabla\rho^{\downarrow}|)$
  29 !   for a muffin-tin charge density, as required by the generalised gradient
  30 !   approximation functional for spin-polarised densities. In the case of spin
  31 !   unpolarised calculations, $|\nabla\rho|$, $\nabla^2\rho$ and
  32 !   $\nabla\rho\cdot(\nabla|\nabla\rho|)$ are returned in the arrays
  33 !   {\tt gupmt}, {\tt g2upmt} and {\tt g3upmt}, respectively, while
  34 !   {\tt grhomt}, {\tt gdnmt}, {\tt g2dnmt}, {\tt g3rhomt} and {\tt g3dnmt} are
  35 !   not referenced. The input densities are in terms of real spherical harmonic
  36 !   expansions but the returned functions are in spherical coordinates. See
  37 !   routines {\tt potxc}, {\tt modxcifc}, {\tt gradrfmt}, {\tt genrlm} and
  38 !   {\tt genshtmat}.
  39 !
  40 ! !REVISION HISTORY:
  41 !   Created April 2004 (JKD)
  42 !EOP
  43 !BOC
  44 implicit none
  45 ! arguments
  46 integer, intent(in) :: is
  47 integer, intent(in) :: ia
  48 real(8), intent(out) :: grhomt(lmmaxvr,nrmtmax)
  49 real(8), intent(out) :: gupmt(lmmaxvr,nrmtmax)
  50 real(8), intent(out) :: gdnmt(lmmaxvr,nrmtmax)
  51 real(8), intent(out) :: g2upmt(lmmaxvr,nrmtmax)
  52 real(8), intent(out) :: g2dnmt(lmmaxvr,nrmtmax)
  53 real(8), intent(out) :: g3rhomt(lmmaxvr,nrmtmax)
  54 real(8), intent(out) :: g3upmt(lmmaxvr,nrmtmax)
  55 real(8), intent(out) :: g3dnmt(lmmaxvr,nrmtmax)
  56 ! local variables
  57 integer ias,nr,ir,i,itp
  58 ! allocatable arrays
  59 real(8), allocatable :: rfmt1(:,:)
  60 real(8), allocatable :: rftp1(:,:,:)
  61 real(8), allocatable :: rftp2(:,:,:)
  62 real(8), allocatable :: rftp3(:)
  63 real(8), allocatable :: grfmt1(:,:,:)
  64 real(8), allocatable :: grfmt2(:,:,:)
  65 allocate(rfmt1(lmmaxvr,nrmtmax))
  66 allocate(rftp1(lmmaxvr,nrmtmax,3))
  67 allocate(rftp2(lmmaxvr,nrmtmax,3))
  68 allocate(rftp3(lmmaxvr))
  69 allocate(grfmt1(lmmaxvr,nrmtmax,3))
  70 allocate(grfmt2(lmmaxvr,nrmtmax,3))
  71 ias=idxas(ia,is)
  72 nr=nrmt(is)
  73 if (spinpol) then
  74 ! rhoup for spin-polarised case
  75   rfmt1(:,1:nr)=0.5d0*(rhomt(:,1:nr,ias)+magmt(:,1:nr,ias,ndmag))
  76 else
  77 ! rho for spin-unpolarised case
  78   rfmt1(:,1:nr)=rhomt(:,1:nr,ias)
  79 end if
  80 ! |grad rhoup|
  81 call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,rfmt1,grfmt1)
  82 do ir=1,nr
  83   do i=1,3
  84     call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt1(1,ir,i),1, &
  85      0.d0,rftp1(1,ir,i),1)
  86   end do
  87   do itp=1,lmmaxvr
  88     gupmt(itp,ir)=sqrt(rftp1(itp,ir,1)**2+rftp1(itp,ir,2)**2+rftp1(itp,ir,3)**2)
  89   end do
  90 end do
  91 ! grad^2 rhoup
  92 g2upmt(:,1:nr)=0.d0
  93 do i=1,3
  94   call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,grfmt1(1,1,i),grfmt2)
  95   do ir=1,nr
  96     call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt2(1,ir,i),1, &
  97      1.d0,g2upmt(1,ir),1)
  98   end do
  99 end do
 100 ! (grad rhoup).(grad |grad rhoup|)
 101 do ir=1,nr
 102   call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rfshtvr,lmmaxvr,gupmt(1,ir),1,0.d0, &
 103    rfmt1(1,ir),1)
 104 end do
 105 call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,rfmt1,grfmt1)
 106 g3upmt(:,1:nr)=0.d0
 107 do i=1,3
 108   do ir=1,nr
 109     call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt1(1,ir,i),1, &
 110      0.d0,rftp3,1)
 111     do itp=1,lmmaxvr
 112       g3upmt(itp,ir)=g3upmt(itp,ir)+rftp1(itp,ir,i)*rftp3(itp)
 113     end do
 114   end do
 115 end do
 116 if (spinpol) then
 117 ! rhodn
 118   rfmt1(:,1:nr)=0.5d0*(rhomt(:,1:nr,ias)-magmt(:,1:nr,ias,ndmag))
 119 ! |grad rhodn|
 120   call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,rfmt1,grfmt1)
 121   do ir=1,nr
 122     do i=1,3
 123       call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt1(1,ir,i),1, &
 124        0.d0,rftp2(1,ir,i),1)
 125     end do
 126     do itp=1,lmmaxvr
 127       gdnmt(itp,ir)=sqrt(rftp2(itp,ir,1)**2+rftp2(itp,ir,2)**2 &
 128        +rftp2(itp,ir,3)**2)
 129     end do
 130   end do
 131 ! grad^2 rhodn
 132   g2dnmt(:,1:nr)=0.d0
 133   do i=1,3
 134     call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,grfmt1(1,1,i),grfmt2)
 135     do ir=1,nr
 136       call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt2(1,ir,i),1, &
 137        1.d0,g2dnmt(1,ir),1)
 138     end do
 139   end do
 140 ! (grad rhodn).(grad |grad rhodn|)
 141   do ir=1,nr
 142     call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rfshtvr,lmmaxvr,gdnmt(1,ir),1,0.d0, &
 143      rfmt1(1,ir),1)
 144   end do
 145   call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,rfmt1,grfmt1)
 146   g3dnmt(:,1:nr)=0.d0
 147   do i=1,3
 148     do ir=1,nr
 149       call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt1(1,ir,i),1, &
 150        0.d0,rftp3,1)
 151       do itp=1,lmmaxvr
 152         g3dnmt(itp,ir)=g3dnmt(itp,ir)+rftp2(itp,ir,i)*rftp3(itp)
 153       end do
 154     end do
 155   end do
 156 ! |grad rho|
 157   do ir=1,nr
 158     do itp=1,lmmaxvr
 159       grhomt(itp,ir)=sqrt((rftp1(itp,ir,1)+rftp2(itp,ir,1))**2 &
 160                          +(rftp1(itp,ir,2)+rftp2(itp,ir,2))**2 &
 161                          +(rftp1(itp,ir,3)+rftp2(itp,ir,3))**2)
 162     end do
 163   end do
 164 ! (grad rho).(grad |grad rho|)
 165   do ir=1,nr
 166     call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rfshtvr,lmmaxvr,grhomt(1,ir),1,0.d0, &
 167      rfmt1(1,ir),1)
 168   end do
 169   call gradrfmt(lmaxvr,nr,spr(1,is),lmmaxvr,nrmtmax,rfmt1,grfmt1)
 170   g3rhomt(:,1:nr)=0.d0
 171   do i=1,3
 172     do ir=1,nr
 173       call dgemv('N',lmmaxvr,lmmaxvr,1.d0,rbshtapw,lmmaxapw,grfmt1(1,ir,i),1, &
 174        0.d0,rftp3,1)
 175       do itp=1,lmmaxvr
 176         g3rhomt(itp,ir)=g3rhomt(itp,ir) &
 177          +(rftp1(itp,ir,i)+rftp2(itp,ir,i))*rftp3(itp)
 178       end do
 179     end do
 180   end do
 181 end if
 182 deallocate(rfmt1,rftp1,rftp2,rftp3,grfmt1,grfmt2)
 183 return
 184 end subroutine
 185 !EOC