Commit ce379ccf authored by Daniel Wortmann's avatar Daniel Wortmann

Added postprocessing code for libxc, major refactoring in x_tofrom_grid routines needed as well

parent b7565d63
......@@ -39,7 +39,9 @@ MODULE m_types_xcpot
REAL,ALLOCATABLE :: gggrd(:),grgru(:),grgrd(:)
!These are the contracted Gradients used in libxc
REAL,ALLOCATABLE :: sigma(:,:)
REAL,ALLOCATABLE :: vsigma(:,:)
REAL,ALLOCATABLE :: gr(:,:,:)
REAL,ALLOCATABLE :: laplace(:,:)
END TYPE t_gradients
CONTAINS
......@@ -64,15 +66,14 @@ CONTAINS
a_ex=-1
END FUNCTION xcpot_get_exchange_weight
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh,vxc,vx,grad,drdsigma)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh,vxc,vx,grad)
CLASS(t_xcpot),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
!--> charge density
REAL,INTENT (IN) :: rh(:,:)
!---> xc potential
REAL, INTENT (OUT) :: vxc (:,:),vx(:,:)
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:)
TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
END SUBROUTINE xcpot_get_vxc
......
......@@ -135,7 +135,7 @@ CONTAINS
END FUNCTION xcpot_get_exchange_weight
!***********************************************************************
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad,drdsigma)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
!***********************************************************************
!
USE m_xcxal, ONLY : vxcxal
......@@ -164,8 +164,7 @@ CONTAINS
REAL, INTENT (OUT) :: vxc(:,:)
! optional arguments for GGA
TYPE(t_gradients),INTENT(IN),OPTIONAL::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:) !This will not be allocated
TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad
!c
!c ---> local scalars
INTEGER :: ngrid
......
......@@ -95,7 +95,7 @@ CONTAINS
END FUNCTION xcpot_get_exchange_weight
!***********************************************************************
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad,drdsigma)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
!***********************************************************************
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
......@@ -104,22 +104,21 @@ CONTAINS
REAL, INTENT (OUT) :: vx (:,:) !points,spin
REAL, INTENT (OUT ) :: vxc(:,:) !
! optional arguments for GGA
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:)
TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
#ifdef CPP_LIBXC
REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:)
REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
!libxc uses the spin as a first index, hence we have to transpose....
ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
IF (xcpot%is_gga()) THEN
CALL judft_error("libxc GGA not implemented yet")
IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
ALLOCATE(drdsigma(SIZE(grad%sigma)))
CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,drdsigma)
ALLOCATE(grad%vsigma,mold=grad%sigma)
CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,grad%vsigma)
IF (xcpot%func_id_c>0) THEN
ALLOCATE(vsigma(SIZE(grad%sigma)))
ALLOCATE(vsigma,mold=grad%sigma)
CALL xc_f03_gga_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,vsigma)
drdsigma=drdsigma+vsigma
grad%vsigma=grad%vsigma+vsigma
vxc_tmp=vxc_tmp+vx_tmp
ENDIF
ELSE !LDA potentials
......
......@@ -2,8 +2,6 @@ set(fleur_F77 ${fleur_F77}
vgen/dylm3.f
vgen/fft3dxc.f
vgen/grdrsvac.f
vgen/mkgxyz3.f
vgen/mkgylm.f
vgen/mkgz.f
vgen/modcyli.f
vgen/modcylk.f
......@@ -17,6 +15,8 @@ vgen/visp5_z.f
)
set(fleur_F90 ${fleur_F90}
vgen/b_field.F90
vgen/mkgylm.f90
vgen/mkgxyz3.f90
vgen/convol.f90
vgen/grdrsis.f90
vgen/mt_tofrom_grid.F90
......
......@@ -3,75 +3,85 @@
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mkgxyz3
c.....------------------------------------------------------------------
c by use of cartesian x,y,z components of charge density gradients,
c make the quantities
cc agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
cc used to calculate gradient contribution to xc potential and
cc energy.
c.....------------------------------------------------------------------
CONTAINS
SUBROUTINE mkgxyz3(
> ndm,jsdm,ng3,jspins,vl,
> dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,
< grad)
USE m_types
IMPLICIT NONE
INTEGER, INTENT (IN) :: ndm,ng3,jsdm,jspins
REAL, INTENT (IN) :: vl(ndm,jsdm)
REAL, INTENT (IN) :: dvx(ndm,jsdm),dvy(ndm,jsdm),dvz(ndm,jsdm)
REAL, INTENT (IN) :: dvxx(ndm,jsdm),dvyy(ndm,jsdm),dvzz(ndm,jsdm)
REAL, INTENT (IN) :: dvyz(ndm,jsdm),dvzx(ndm,jsdm),dvxy(ndm,jsdm)
TYPE(t_gradients),INTENT(INOUT)::grad
REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,
& vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,
& vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,
& dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,
+ dagrzu,dzdx,dzdy,dzdz,
+ sml
INTEGER i
sml = 1.e-14
IF(ALLOCATED(grad%sigma)) THEN
!Use only contracted gradients for libxc
if (jspins==1) THEN
DO i=1,ng3
grad%sigma(1,i)=
+ dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
ENDDO
ELSE
DO i=1,ng3
grad%sigma(1,i)=
+ dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
grad%sigma(2,i)=
+ dvx(i,1)*dvx(i,2)+dvy(i,1)*dvy(i,2)+dvz(i,1)*dvz(i,2)
grad%sigma(3,i)=
+ dvx(i,2)*dvx(i,2)+dvy(i,2)*dvy(i,2)+dvz(i,2)*dvz(i,2)
ENDDO
ENDIF
RETURN
ENDIF
DO i = 1,ndm
grad%agrt(i) = 0.0
grad%agru(i) = 0.0
grad%agrd(i) = 0.0
grad%gggrt(i) = 0.0
grad%gggru(i) = 0.0
grad%gggrd(i) = 0.0
grad%gzgr(i) = 0.0
grad%g2rt(i) = 0.0
grad%g2ru(i) = 0.0
grad%g2rd(i) = 0.0
ENDDO
IF (jspins.eq.1) THEN
DO 10 i = 1,ng3
MODULE m_mkgxyz3
!.....------------------------------------------------------------------
!c by use of cartesian x,y,z components of charge density gradients,
!c make the quantities
!cc agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
!cc used to calculate gradient contribution to xc potential and
!cc energy.
!c.....------------------------------------------------------------------
CONTAINS
SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,grad)
USE m_types
IMPLICIT NONE
REAL, INTENT (IN) :: vl(:,:)
REAL, INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
REAL, INTENT (IN) :: dvyz(:,:),dvzx(:,:),dvxy(:,:)
TYPE(t_gradients),INTENT(INOUT)::grad
REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,&
vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,&
vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,&
dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
dagrzu,dzdx,dzdy,dzdz,sml
INTEGER i,js,jspins,nsp
nsp=SIZE(dvx,1)
jspins=SIZE(dvx,2)
sml = 1.e-14
IF (ALLOCATED(grad%gr)) THEN
! Gradients for libxc
DO js=1,jspins
DO i=1,nsp
grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
ENDDO
END DO
IF(ALLOCATED(grad%sigma)) THEN
!Use only contracted gradients for libxc
IF (jspins==1) THEN
DO i=1,nsp
grad%sigma(1,i)= dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
ENDDO
ELSE
DO i=1,nsp
grad%sigma(1,i)= dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
grad%sigma(2,i)= dvx(i,1)*dvx(i,2)+dvy(i,1)*dvy(i,2)+dvz(i,1)*dvz(i,2)
grad%sigma(3,i)= dvx(i,2)*dvx(i,2)+dvy(i,2)*dvy(i,2)+dvz(i,2)*dvz(i,2)
ENDDO
ENDIF
END IF
IF(ALLOCATED(grad%laplace)) THEN
DO js=1,jspins
DO i=1,nsp
grad%laplace(i,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
ENDDO
ENDDO
ENDIF
RETURN
ENDIF
IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
DO i = 1,size(grad%agrt)
grad%agrt(i) = 0.0
grad%agru(i) = 0.0
grad%agrd(i) = 0.0
grad%gggrt(i) = 0.0
grad%gggru(i) = 0.0
grad%gggrd(i) = 0.0
grad%gzgr(i) = 0.0
grad%g2rt(i) = 0.0
grad%g2ru(i) = 0.0
grad%g2rd(i) = 0.0
ENDDO
IF (jspins.eq.1) THEN
DO 10 i = 1,nsp
vlu=max(vl(i,1)/2,sml)
dvxu=dvx(i,1)/2
......@@ -108,7 +118,7 @@ c.....------------------------------------------------------------------
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
c agr: abs(grad(ro)), t,u,d for total, up and down.
! agr: abs(grad(ro)), t,u,d for total, up and down.
grad%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml)
......@@ -130,28 +140,28 @@ c agr: abs(grad(ro)), t,u,d for total, up and down.
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
c dzdx=d(zeta)/dx,..
! dzdx=d(zeta)/dx,..
dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
c gzgr=grad(zeta)*grad(ro).
! gzgr=grad(zeta)*grad(ro).
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
c g2r: grad(grad(ro))
! g2r: grad(grad(ro))
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
10 ENDDO
10 ENDDO
ELSE
ELSE
DO 20 i = 1,ng3
DO 20 i = 1,nsp
vlu = max(vl(i,1),sml)
dvxu=dvx(i,1)
......@@ -187,7 +197,7 @@ c g2r: grad(grad(ro))
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
c agr: abs(grad(ro)), t,u,d for total, up and down.
!c agr: abs(grad(ro)), t,u,d for total, up and down.
grad%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml)
......@@ -209,25 +219,25 @@ c agr: abs(grad(ro)), t,u,d for total, up and down.
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
c dzdx=d(zeta)/dx,..
!c dzdx=d(zeta)/dx,..
dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
c gzgr=grad(zeta)*grad(ro).
!c gzgr=grad(zeta)*grad(ro).
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
c g2r: grad(grad(ro))
!c g2r: grad(grad(ro))
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
20 ENDDO
20 ENDDO
ENDIF
END SUBROUTINE mkgxyz3
END MODULE m_mkgxyz3
ENDIF
END SUBROUTINE mkgxyz3
END MODULE m_mkgxyz3
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mkgylm
CONTAINS
SUBROUTINE mkgylm(
> jspins,rv,thet,nsp,
> rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff,
> rhdtf,rhdrt,rhdrf,grad,kt)
c.....------------------------------------------------------------------
c by use of charge density and its polar coord. gradient components
c c calculate agr and others used to evaluate gradient
c c contributions to potential and energy. t.a. 1996.
c.....------------------------------------------------------------------
c ro=sum(ro*ylh), rdr=sum(drr*ylh), drdr=sum(ddrr*ylh),
c c rdt=sum(ro*ylht1), rdtt=sum(ro*ylht2), ...
c c rdf=sum(ro*ylhf1), rdff=sum(ro*ylhf2), ...
c c rdtf=sum(ro*ylhtf), rdff=sum(ro*ylhf2), ...
c c rdrt=sum(drr*ylht1),rdrf=sum(drr*ylhf1),
c agr: abs(grad(ro)), g2r: laplacian(ro),
c c gggr: grad(ro)*grad(agr),
c c grgru,d: grad(ro)*grad(rou),for rod., gzgr: grad(zeta)*grad(ro).
c dagrr,-t,-f: d(agr)/dr, d(agr)/dth/r, d(agr)/dfi/r/sint.
c.....------------------------------------------------------------------
USE m_types
IMPLICIT NONE
C ..
C .. Arguments ..
REAL, INTENT (IN) :: rv
REAL, INTENT (IN) :: thet(:)
REAL, INTENT (IN) :: rh(:,:),rhdr(:,:)
REAL, INTENT (IN) :: rhdf(:,:),rhdrr(:,:)
REAL, INTENT (IN) :: rhdtt(:,:),rhdff(:,:)
REAL, INTENT (IN) :: rhdtf(:,:),rhdrt(:,:)
REAL, INTENT (IN) :: rhdrf(:,:),rhdt(:,:)
TYPE(t_gradients),INTENT(INOUT) :: grad
INTEGER,INTENT(IN) :: jspins,nsp
INTEGER,INTENT(IN) :: kt !index of first point to use in gradients
C ..
C .. Locals ..
INTEGER i
REAL chsml, dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt,
+ dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr,
+ grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu,
+ rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu,
+ rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu,
+ ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1
chsml = 1.e-10
rv1 = rv
rv2 = rv1**2
rv3 = rv1**3
IF (ALLOCATED(grad%sigma)) THEN
! Contracted gradients for libxc
IF (jspins==1) THEN
DO i=1,nsp
grad%gr(:,kt+i,1)=(/rhdr(i,1),rhdt(i,1),rhdf(i,1)/)
grad%sigma(1,kt+i)=
+ dot_product(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
ENDDO
ELSE
DO i=1,nsp
grad%gr(:,kt+i,1)=(/rhdr(i,1),rhdt(i,1),rhdf(i,1)/)
grad%gr(:,kt+i,2)=(/rhdr(i,2),rhdt(i,2),rhdf(i,2)/)
grad%sigma(1,kt+i)=
+ dot_product(grad%gr(:,i,1),grad%gr(:,i,1))
grad%sigma(2,kt+i)=
+ dot_product(grad%gr(:,i,1),grad%gr(:,i,2))
grad%sigma(3,kt+i)=
+ dot_product(grad%gr(:,i,2),grad%gr(:,i,2))
ENDDO
ENDIF
ELSE
! Old code for in-build xcpots
IF (jspins.EQ.1) THEN
points_1 : DO i = 1,nsp
grad%agrt(kt+i) = 0.0
grad%agru(kt+i) = 0.0
grad%agrd(kt+i) = 0.0
grad%g2rt(kt+i) = 0.0
grad%g2ru(kt+i) = 0.0
grad%g2rd(kt+i) = 0.0
grad%gggrt(kt+i) = 0.0
grad%gggru(kt+i) = 0.0
grad%gggrd(kt+i) = 0.0
grad%grgru(kt+i) = 0.0
grad%grgrd(kt+i) = 0.0
grad%gzgr(kt+i) = 0.0
ro = rh(i,1)
IF (ro.LT.chsml) CYCLE points_1
sint1 = sin(thet(i))
sint2 = sint1**2
tant1 = tan(thet(i))
rvsin1 = rv1*sint1
rou = ro/2
rdru = rhdr(i,1)/2
rdtu = rhdt(i,1)/2
rdfu = rhdf(i,1)/2
rdrru = rhdrr(i,1)/2
rdttu = rhdtt(i,1)/2
rdffu = rhdff(i,1)/2
rdtfu = rhdtf(i,1)/2
rdrtu = rhdrt(i,1)/2
rdrfu = rhdrf(i,1)/2
rod = rou
rdrd = rdru
rdtd = rdtu
rdfd = rdfu
rdrrd = rdrru
rdttd = rdttu
rdffd = rdffu
rdtfd = rdtfu
rdrtd = rdrtu
rdrfd = rdrfu
rdr = rdru + rdrd
rdt = rdtu + rdtd
rdf = rdfu + rdfd
drdr = rdrru + rdrrd
rdtt = rdttu + rdttd
rdff = rdffu + rdffd
rdtf = rdtfu + rdtfd
rdrt = rdrtu + rdrtd
rdrf = rdrfu + rdrfd
ro2 = ro**2
grr = rdr
grt = rdt/rv1
grf = rdf/rvsin1
grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
IF (grad%agrt(kt+i).LT.chsml) CYCLE points_1
dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+
+ rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/
+ (grad%agrt(kt+i)*rv3)
dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/
+ (grad%agrt(kt+i)*rv3*sint1)
grad%g2rt(kt+i)=drdr+2.0*rdr/rv1+
+ (rdtt+rdt/tant1+rdff/sint2)/rv2
dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
!-- > dzdtr,dzdfs vanish by definition.
dzdtr = 0.0
dzdfs = 0.0
grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
grru = rdru
grtu = rdtu/rv1
grfu = rdfu/rvsin1
grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+
+ rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+
+ rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/
+ (grad%agru(kt+i)*rv3*sint1)
grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 +
+ (rdttu+rdtu/tant1+rdffu/sint2)/rv2
grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
grrd = rdrd
grtd = rdtd/rv1
grfd = rdfd/rvsin1
grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+
+ rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+
+ rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/
+ (grad%agrd(kt+i)*rv3*sint1)
grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 +
+ (rdttd+rdtd/tant1+rdffd/sint2)/rv2
grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
ENDDO points_1
ELSE
points_2 : DO i = 1,nsp
grad%agrt(kt+i) = 0.0
grad%agru(kt+i) = 0.0
grad%agrd(kt+i) = 0.0
grad%g2rt(kt+i) = 0.0
grad%g2ru(kt+i) = 0.0
grad%g2rd(kt+i) = 0.0
grad%gggrt(kt+i) = 0.0
grad%gggru(kt+i) = 0.0
grad%gggrd(kt+i) = 0.0
grad%grgru(kt+i) = 0.0
grad%grgrd(kt+i) = 0.0
grad%gzgr(kt+i) = 0.0
ro = rh(i,1) + rh(i,jspins)
IF (ro.LT.chsml) CYCLE points_2
sint1 = sin(thet(i))
sint2 = sint1**2
tant1 = tan(thet(i))
rvsin1 = rv1*sint1
rou = rh(i,1)
rdru = rhdr(i,1)
rdtu = rhdt(i,1)
rdfu = rhdf(i,1)
rdrru = rhdrr(i,1)
rdttu = rhdtt(i,1)
rdffu = rhdff(i,1)