Commit 9b23d5db authored by Alexander Neukirchen's avatar Alexander Neukirchen

Gradient also seems to be working just fine. Some beautification in old code.

parent 1b2ec199
......@@ -101,7 +101,7 @@ CONTAINS
TYPE(t_wann) :: wann
TYPE(t_potden) :: vTot, vx, vCoul, vTemp, vxcForPlotting
TYPE(t_potden) :: inDen, outDen, EnergyDen, dummyDen
TYPE(t_potden), DIMENSION(3) :: testDen
TYPE(t_potden), DIMENSION(3) :: testDen, testGrad
CLASS(t_xcpot), ALLOCATABLE :: xcpot
CLASS(t_forcetheo), ALLOCATABLE :: forcetheo
......@@ -517,14 +517,20 @@ CONTAINS
! Test: Build a field, for which the theoretical divergence etc. are known and
! compare with the result of the routine.
!CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeBxc(stars,atoms,sphhar,vacuum,input,noco,vTot,testDen)
CALL matrixsplit(stars, atoms, sphhar, vacuum, input, noco, 1.0, inDen, dummyDen, testDen(1), testDen(2), testDen(3))
CALL checkplotinp()
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', dummyDen, testDen(1), testDen(2), testDen(3))
!CALL matrixsplit(stars, atoms, sphhar, vacuum, input, noco, 1.0, inDen, dummyDen, testDen(1), testDen(2), testDen(3))
!CALL checkplotinp()
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', dummyDen, testDen(1), testDen(2), testDen(3))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDeny ', testDen(2))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDenz ', testDen(3))
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,testDen)
!CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,testDen)
DO i=1,3
CALL testGrad(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(testGrad(i)%pw_w,mold=testGrad(i)%pw)
ENDDO
CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,testDen(2),testGrad)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testGrad ', testGrad(1), testGrad(1), testGrad(2), testGrad(3))
CALL add_usage_data("Iterations",iter)
......
......@@ -67,6 +67,9 @@ MODULE m_types_xcpot
REAL,ALLOCATABLE :: vsigma(:,:)
REAL,ALLOCATABLE :: gr(:,:,:)
REAL,ALLOCATABLE :: laplace(:,:)
!For sourcefree calculations, the cartesian components of the gradient
!on the spherical grid are needed (no spin dependency needed though).
REAL,ALLOCATABLE :: grxyz(:)
END TYPE t_gradients
CONTAINS
SUBROUTINE kED_alloc_mt(kED,nsp_x_jmtd, jspins, n_start, n_types, n_stride)
......
......@@ -4,243 +4,239 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mkgxyz3
USE m_judft
!.....------------------------------------------------------------------
!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.....------------------------------------------------------------------
USE m_judft
!-----------------------------------------------------------------------------
! Using the cartesian components and derivatives of a charge density rho on
! the real space grid, make the following quantaties:
!
! gr(js): grad(rho_js)
! sigma: |grad(rho)|^2 for jspins==1, otherwise three components, namely
! |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2
! laplace(js): laplace(rho_js)
!
! and these older components:
! agrt/u/d: |grad(rho)| for total density/spin-up/spin-down
! g2rt/u/d: laplace(rho)
! gggrt/u/d: (grad(rho))*(grad(|grad(rho)|)) [scalar product]
! gzgr: (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho
!
! which are used to calculate gradient contributions to the xc potential and
! energy.
!
! vl is rho, dv[i][j] are the partial derivatives along one/two directions.
!
! Modified so only allocated old quantities require calculations. A.N. 2019
!
! Quantities fo libxc are calculated as well. D.W./M.R. 2018
!
! Original script by T.A. 1996
!-----------------------------------------------------------------------------
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
REAL, PARAMETER :: sml = 1.e-14
INTEGER i,js,jspins,nsp
nsp=SIZE(dvx,1)
jspins=SIZE(dvx,2)
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")
IF(ALLOCATED(grad%agrt)) THEN
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
END IF
IF (jspins.eq.1) THEN
DO i = 1,nsp
vlu = max(vl(i,1)/2,sml)
dvxu = dvx(i,1)/2
dvyu = dvy(i,1)/2
dvzu = dvz(i,1)/2
dvxxu = dvxx(i,1)/2
dvyyu = dvyy(i,1)/2
dvzzu = dvzz(i,1)/2
dvyzu = dvyz(i,1)/2
dvzxu = dvzx(i,1)/2
dvxyu = dvxy(i,1)/2
vld = vlu
dvxd = dvxu
dvyd = dvyu
dvzd = dvzu
dvxxd = dvxxu
dvyyd = dvyyu
dvzzd = dvzzu
dvyzd = dvyzu
dvzxd = dvzxu
dvxyd = dvxyu
vlt = vlu + vld
dvxt = dvxu + dvxd
dvyt = dvyu + dvyd
dvzt = dvzu + dvzd
dvxxt = dvxxu + dvxxd
dvyyt = dvyyu + dvyyd
dvzzt = dvzzu + dvzzd
dvyzt = dvyzu + dvyzd
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
IF(ALLOCATED(grad%agrt)) THEN
! 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)
grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
! 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
! gzgr=grad(zeta)*grad(ro).
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
! g2r: grad(grad(ro))
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
END IF
ENDDO
ELSE
DO i = 1,nsp
vlu = max(vl(i,1),sml)
dvxu = dvx(i,1)
dvyu = dvy(i,1)
dvzu = dvz(i,1)
dvxxu = dvxx(i,1)
dvyyu = dvyy(i,1)
dvzzu = dvzz(i,1)
dvyzu = dvyz(i,1)
dvzxu = dvzx(i,1)
dvxyu = dvxy(i,1)
vld = max(vl(i,jspins),sml)
dvxd = dvx(i,jspins)
dvyd = dvy(i,jspins)
dvzd = dvz(i,jspins)
dvxxd = dvxx(i,jspins)
dvyyd = dvyy(i,jspins)
dvzzd = dvzz(i,jspins)
dvyzd = dvyz(i,jspins)
dvzxd = dvzx(i,jspins)
dvxyd = dvxy(i,jspins)
vlt = vlu + vld
dvxt = dvxu + dvxd
dvyt = dvyu + dvyd
dvzt = dvzu + dvzd
dvxxt = dvxxu + dvxxd
dvyyt = dvyyu + dvyyd
dvzzt = dvzzu + dvzzd
dvyzt = dvyzu + dvyzd
dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd
IF(ALLOCATED(grad%agrt)) THEN
!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)
grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
!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).
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
!c g2r: grad(grad(ro))
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
END IF
ENDDO
ENDIF
END SUBROUTINE mkgxyz3
SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvxz,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(:,:),dvxz(:,:),dvxy(:,:)
TYPE(t_gradients), INTENT(INOUT) :: grad
REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvxzt,dvxyt, &
vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvxzu,dvxyu, &
vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvxzd,dvxyd, &
dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
dagrzu,dzdx,dzdy,dzdz
REAL, PARAMETER :: sml = 1.e-14
INTEGER i,js,jspins,nsp
nsp=SIZE(dvx,1)
jspins=SIZE(dvx,2)
! Gradients for libxc and sourcefree calculations.
IF (ALLOCATED(grad%gr)) THEN
DO js=1,jspins
DO i=1,nsp
grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
END DO
END DO
! Use contracted gradients only for libxc.
IF(ALLOCATED(grad%sigma)) THEN
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)
END DO
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)
END DO
END IF
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)
END DO
END DO
END IF
RETURN ! Do not calculate arrays for in-build GGA.
END IF
IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
IF(ALLOCATED(grad%agrt)) THEN
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
END DO
IF (jspins.eq.1) THEN
DO i = 1,nsp
vlu = max(vl(i,1)/2,sml)
dvxu = dvx(i,1)/2
dvyu = dvy(i,1)/2
dvzu = dvz(i,1)/2
dvxxu = dvxx(i,1)/2
dvyyu = dvyy(i,1)/2
dvzzu = dvzz(i,1)/2
dvyzu = dvyz(i,1)/2
dvxzu = dvxz(i,1)/2
dvxyu = dvxy(i,1)/2
vld = vlu
dvxd = dvxu
dvyd = dvyu
dvzd = dvzu
dvxxd = dvxxu
dvyyd = dvyyu
dvzzd = dvzzu
dvyzd = dvyzu
dvxzd = dvxzu
dvxyd = dvxyu
vlt = vlu + vld
dvxt = dvxu + dvxd
dvyt = dvyu + dvyd
dvzt = dvzu + dvzd
dvxxt = dvxxu + dvxxd
dvyyt = dvyyu + dvyyd
dvzzt = dvzzu + dvzzd
dvyzt = dvyzu + dvyzd
dvxzt = dvxzu + dvxzd
dvxyt = dvxyu + dvxyd
grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(i)
dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(i)
dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(i)
dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
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
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
ENDDO
ELSE
DO i = 1,nsp
vlu = max(vl(i,1),sml)
dvxu = dvx(i,1)
dvyu = dvy(i,1)
dvzu = dvz(i,1)
dvxxu = dvxx(i,1)
dvyyu = dvyy(i,1)
dvzzu = dvzz(i,1)
dvyzu = dvyz(i,1)
dvxzu = dvxz(i,1)
dvxyu = dvxy(i,1)
vld = max(vl(i,jspins),sml)
dvxd = dvx(i,jspins)
dvyd = dvy(i,jspins)
dvzd = dvz(i,jspins)
dvxxd = dvxx(i,jspins)
dvyyd = dvyy(i,jspins)
dvzzd = dvzz(i,jspins)
dvyzd = dvyz(i,jspins)
dvxzd = dvxz(i,jspins)
dvxyd = dvxy(i,jspins)
vlt = vlu + vld
dvxt = dvxu + dvxd
dvyt = dvyu + dvyd
dvzt = dvzu + dvzd
dvxxt = dvxxu + dvxxd
dvyyt = dvyyu + dvyyd
dvzzt = dvzzu + dvzzd
dvyzt = dvyzu + dvyzd
dvxzt = dvxzu + dvxzd
dvxyt = dvxyu + dvxyd
grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(i)
dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(i)
dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(i)
dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
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
grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
grad%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd
END DO
END IF
END IF
END SUBROUTINE mkgxyz3
END MODULE m_mkgxyz3
......@@ -4,76 +4,81 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mkgylm
!-----------------------------------------------------------------------------
! Using the components and derivatives of a charge density rho in spherical
! coordinates on the real space grid, make the following quantaties:
!
! gr(js): [partial_r (rho),partial_theta (rho),partial_phi (rho)]
! sigma: |grad(rho)|^2 for jspins==1, otherwise three components, namely
! |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2
! laplace(js): laplace(rho_js)
!
! and these older components:
! agrt/u/d: |grad(rho)| for total density/spin-up/spin-down
! g2rt/u/d: laplace(rho)
! gggrt/u/d: (grad(rho))*(grad(|grad(rho)|)) [scalar product]
! gzgr: (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho
!
! which are used to calculate gradient contributions to the xc potential and
! energy.
!
! rh is rho, rhd[i][j] are the partial derivatives along one/two directions.
! rv and thet are the radius and polar angles and t/f in derivatives stand for
! theta/phi.
!
! Modified so only allocated old quantities require calculations. A.N. 2019
!
! Quantities fo libxc are calculated as well. D.W./M.R. 2018
!
! Original script by T.A. 1996
!-----------------------------------------------------------------------------
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