Commit 1018338f authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'xc-pot-refactor' into 'develop'

Merge libxc branch as GGA seems to be approximately working...

See merge request fleur/fleur!9
parents 295f78fc c941f55d
...@@ -82,6 +82,7 @@ ...@@ -82,6 +82,7 @@
! arltv(i) : length of reciprical lattice vector along ! arltv(i) : length of reciprical lattice vector along
! direction (i) ! direction (i)
! !
IF (.NOT.xcpot%is_gga()) xcpot%gmaxxc=stars%gmax
WRITE (6,'('' gmaxxc should be: 2*kmax <= gmaxxc <= gmax '')') WRITE (6,'('' gmaxxc should be: 2*kmax <= gmaxxc <= gmax '')')
IF ( abs( xcpot%gmaxxc - stars%gmax ) .le. 10.0**(-6) ) THEN IF ( abs( xcpot%gmaxxc - stars%gmax ) .le. 10.0**(-6) ) THEN
WRITE (6,'('' concerning memory, you may want to choose'',& WRITE (6,'('' concerning memory, you may want to choose'',&
......
...@@ -20,7 +20,7 @@ CONTAINS ...@@ -20,7 +20,7 @@ CONTAINS
!! TE_EXC : charge density-ex-corr.energy density integral !! TE_EXC : charge density-ex-corr.energy density integral
SUBROUTINE vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,& SUBROUTINE vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,&
obsolete,cell,oneD,sliceplot,mpi,results,noco,den,vTot,vXC,vCoul) obsolete,cell,oneD,sliceplot,mpi,results,noco,den,vTot,vx,vCoul)
USE m_rotate_int_den_to_local USE m_rotate_int_den_to_local
USE m_bfield USE m_bfield
...@@ -33,6 +33,7 @@ CONTAINS ...@@ -33,6 +33,7 @@ CONTAINS
#endif #endif
IMPLICIT NONE IMPLICIT NONE
TYPE(t_results), INTENT(INOUT) :: results TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_hybrid), INTENT(IN) :: hybrid TYPE(t_hybrid), INTENT(IN) :: hybrid
...@@ -51,26 +52,29 @@ CONTAINS ...@@ -51,26 +52,29 @@ CONTAINS
TYPE(t_sphhar), INTENT(IN) :: sphhar TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_potden), INTENT(INOUT) :: den TYPE(t_potden), INTENT(INOUT) :: den
TYPE(t_potden), INTENT(INOUT) :: vTot,vXC,vCoul TYPE(t_potden), INTENT(INOUT) :: vTot,vx,vCoul
TYPE(t_potden) :: workden,denRot TYPE(t_potden) :: workden,denRot
if (mpi%irank==0) WRITE (6,FMT=8000) if (mpi%irank==0) WRITE (6,FMT=8000)
8000 FORMAT (/,/,t10,' p o t e n t i a l g e n e r a t o r',/) 8000 FORMAT (/,/,t10,' p o t e n t i a l g e n e r a t o r',/)
CALL vTot%resetPotDen() CALL vTot%resetPotDen()
CALL vCoul%resetPotDen() CALL vCoul%resetPotDen()
CALL vXC%resetPotDen() CALL vx%resetPotDen()
ALLOCATE(vXC%pw_w,vTot%pw_w,mold=vTot%pw) ALLOCATE(vx%pw_w,vTot%pw_w,mold=vTot%pw)
ALLOCATE(vCoul%pw_w(SIZE(den%pw,1),1)) ALLOCATE(vCoul%pw_w(SIZE(den%pw,1),1))
CALL workDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,0) CALL workDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,0)
!sum up both spins in den into workden !sum up both spins in den into workden
CALL den%sum_both_spin(workden) CALL den%sum_both_spin(workden)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,workden,vCoul,results) CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,workden,vCoul,results)
CALL vCoul%copy_both_spin(vTot) CALL vCoul%copy_both_spin(vTot)
IF (noco%l_noco) THEN IF (noco%l_noco) THEN
...@@ -79,19 +83,21 @@ CONTAINS ...@@ -79,19 +83,21 @@ CONTAINS
CALL rotate_int_den_to_local(dimension,sym,stars,atoms,sphhar,vacuum,cell,input,noco,oneD,denRot) CALL rotate_int_den_to_local(dimension,sym,stars,atoms,sphhar,vacuum,cell,input,noco,oneD,denRot)
ENDIF ENDIF
CALL vgen_xcpot(hybrid,input,xcpot,dimension,atoms,sphhar,stars,vacuum,sym,& CALL vgen_xcpot(hybrid,input,xcpot,dimension,atoms,sphhar,stars,vacuum,sym,&
obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vXC,results) obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vx,results)
!ToDo, check if this is needed for more potentials as well... !ToDo, check if this is needed for more potentials as well...
CALL vgen_finalize(atoms,stars,vacuum,sym,noco,input,vTot,denRot) CALL vgen_finalize(atoms,stars,vacuum,sym,noco,input,vTot,denRot)
DEALLOCATE(vcoul%pw_w,vXC%pw_w) DEALLOCATE(vcoul%pw_w,vx%pw_w)
CALL bfield(input,noco,atoms,field,vTot) CALL bfield(input,noco,atoms,field,vTot)
#ifdef CPP_MPI #ifdef CPP_MPI
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vTot) CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vTot)
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vCoul) CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vCoul)
CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vXC) CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vx)
#endif #endif
END SUBROUTINE vgen END SUBROUTINE vgen
......
...@@ -39,6 +39,9 @@ MODULE m_types_xcpot ...@@ -39,6 +39,9 @@ MODULE m_types_xcpot
REAL,ALLOCATABLE :: gggrd(:),grgru(:),grgrd(:) REAL,ALLOCATABLE :: gggrd(:),grgru(:),grgrd(:)
!These are the contracted Gradients used in libxc !These are the contracted Gradients used in libxc
REAL,ALLOCATABLE :: sigma(:,:) REAL,ALLOCATABLE :: sigma(:,:)
REAL,ALLOCATABLE :: vsigma(:,:)
REAL,ALLOCATABLE :: gr(:,:,:)
REAL,ALLOCATABLE :: laplace(:,:)
END TYPE t_gradients END TYPE t_gradients
CONTAINS CONTAINS
...@@ -63,15 +66,14 @@ CONTAINS ...@@ -63,15 +66,14 @@ CONTAINS
a_ex=-1 a_ex=-1
END FUNCTION xcpot_get_exchange_weight 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 CLASS(t_xcpot),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins INTEGER, INTENT (IN) :: jspins
!--> charge density !--> charge density
REAL,INTENT (IN) :: rh(:,:) REAL,INTENT (IN) :: rh(:,:)
!---> xc potential !---> xc potential
REAL, INTENT (OUT) :: vxc (:,:),vx(:,:) REAL, INTENT (OUT) :: vxc (:,:),vx(:,:)
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:)
END SUBROUTINE xcpot_get_vxc END SUBROUTINE xcpot_get_vxc
...@@ -87,8 +89,14 @@ CONTAINS ...@@ -87,8 +89,14 @@ CONTAINS
SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad) SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
INTEGER, INTENT (IN) :: jspins,ngrid INTEGER, INTENT (IN) :: jspins,ngrid
TYPE(t_gradients),INTENT(OUT):: grad TYPE(t_gradients),INTENT(INOUT):: grad
IF (allocated(grad%agrt)) THEN
DEALLOCATE(grad%agrt,grad%agru,grad%agrd)
DEALLOCATE(grad%g2ru,grad%g2rd,grad%gggrt)
DEALLOCATE(grad%gggru,grad%gzgr,grad%g2rt)
DEALLOCATE(grad%gggrd,grad%grgru,grad%grgrd)
ENDIF
!For the in-build xc-pots !For the in-build xc-pots
ALLOCATE(grad%agrt(ngrid),grad%agru(ngrid),grad%agrd(ngrid)) ALLOCATE(grad%agrt(ngrid),grad%agru(ngrid),grad%agrd(ngrid))
ALLOCATE(grad%g2ru(ngrid),grad%g2rd(ngrid),grad%gggrt(ngrid)) ALLOCATE(grad%g2ru(ngrid),grad%g2rd(ngrid),grad%gggrt(ngrid))
......
...@@ -135,7 +135,7 @@ CONTAINS ...@@ -135,7 +135,7 @@ CONTAINS
END FUNCTION xcpot_get_exchange_weight 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 USE m_xcxal, ONLY : vxcxal
...@@ -164,8 +164,7 @@ CONTAINS ...@@ -164,8 +164,7 @@ CONTAINS
REAL, INTENT (OUT) :: vxc(:,:) REAL, INTENT (OUT) :: vxc(:,:)
! optional arguments for GGA ! optional arguments for GGA
TYPE(t_gradients),INTENT(IN),OPTIONAL::grad TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:) !This will not be allocated
!c !c
!c ---> local scalars !c ---> local scalars
INTEGER :: ngrid INTEGER :: ngrid
......
...@@ -44,7 +44,15 @@ CONTAINS ...@@ -44,7 +44,15 @@ CONTAINS
xcpot%jspins=jspins xcpot%jspins=jspins
xcpot%func_id_x=id_x xcpot%func_id_x=id_x
xcpot%func_id_c=id_c xcpot%func_id_c=id_c
if(xcpot%func_id_x == 0 .or. xcpot%func_id_c == 0) then
CALL judft_error("LibXC exchange- and correlation-function indicies need to be set"&
,hint='Try this: ' // ACHAR(10) //&
'<xcFunctional name="libxc" relativisticCorrections="F">' // ACHAR(10) //&
' <libXC exchange="1" correlation="1" /> ' // ACHAR(10) //&
'</xcFunctional> ')
endif
IF (jspins==1) THEN IF (jspins==1) THEN
CALL xc_f03_func_init(xcpot%xc_func_x, xcpot%func_id_x, XC_UNPOLARIZED) CALL xc_f03_func_init(xcpot%xc_func_x, xcpot%func_id_x, XC_UNPOLARIZED)
IF (xcpot%func_id_c>0) CALL xc_f03_func_init(xcpot%xc_func_c, xcpot%func_id_c, XC_UNPOLARIZED) IF (xcpot%func_id_c>0) CALL xc_f03_func_init(xcpot%xc_func_c, xcpot%func_id_c, XC_UNPOLARIZED)
...@@ -95,7 +103,7 @@ CONTAINS ...@@ -95,7 +103,7 @@ CONTAINS
END FUNCTION xcpot_get_exchange_weight 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 IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
...@@ -104,23 +112,23 @@ CONTAINS ...@@ -104,23 +112,23 @@ CONTAINS
REAL, INTENT (OUT) :: vx (:,:) !points,spin REAL, INTENT (OUT) :: vx (:,:) !points,spin
REAL, INTENT (OUT ) :: vxc(:,:) ! REAL, INTENT (OUT ) :: vxc(:,:) !
! optional arguments for GGA ! optional arguments for GGA
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
REAL,ALLOCATABLE,OPTIONAL,INTENT(OUT)::drdsigma(:)
#ifdef CPP_LIBXC #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.... !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(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 ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
IF (xcpot%is_gga()) THEN 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") IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
ALLOCATE(drdsigma(SIZE(grad%sigma))) ALLOCATE(vsigma,mold=grad%vsigma)
CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,drdsigma) !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,vsigma)
IF (xcpot%func_id_c>0) THEN IF (xcpot%func_id_c>0) THEN
ALLOCATE(vsigma(SIZE(grad%sigma))) CALL xc_f03_gga_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma)
CALL xc_f03_gga_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,vsigma) grad%vsigma=grad%vsigma+vsigma
drdsigma=drdsigma+vsigma
vxc_tmp=vxc_tmp+vx_tmp vxc_tmp=vxc_tmp+vx_tmp
ELSE
vxc_tmp=vx_tmp
ENDIF ENDIF
ELSE !LDA potentials ELSE !LDA potentials
CALL xc_f03_lda_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp) CALL xc_f03_lda_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
...@@ -169,9 +177,14 @@ CONTAINS ...@@ -169,9 +177,14 @@ CONTAINS
SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad) SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
INTEGER, INTENT (IN) :: jspins,ngrid INTEGER, INTENT (IN) :: jspins,ngrid
TYPE(t_gradients),INTENT(OUT):: grad TYPE(t_gradients),INTENT(INOUT):: grad
!For libxc we only need the sigma array... !For libxc we only need the sigma array...
IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid)) ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
ALLOCATE(grad%gr(3,ngrid,jspins))
ALLOCATE(grad%laplace(ngrid,jspins))
ALLOCATE(grad%vsigma,mold=grad%sigma)
END SUBROUTINE xcpot_alloc_gradients END SUBROUTINE xcpot_alloc_gradients
......
...@@ -2,8 +2,6 @@ set(fleur_F77 ${fleur_F77} ...@@ -2,8 +2,6 @@ set(fleur_F77 ${fleur_F77}
vgen/dylm3.f vgen/dylm3.f
vgen/fft3dxc.f vgen/fft3dxc.f
vgen/grdrsvac.f vgen/grdrsvac.f
vgen/mkgxyz3.f
vgen/mkgylm.f
vgen/mkgz.f vgen/mkgz.f
vgen/modcyli.f vgen/modcyli.f
vgen/modcylk.f vgen/modcylk.f
...@@ -17,8 +15,11 @@ vgen/visp5_z.f ...@@ -17,8 +15,11 @@ vgen/visp5_z.f
) )
set(fleur_F90 ${fleur_F90} set(fleur_F90 ${fleur_F90}
vgen/b_field.F90 vgen/b_field.F90
vgen/mkgylm.f90
vgen/mkgxyz3.f90
vgen/convol.f90 vgen/convol.f90
vgen/grdrsis.f90 vgen/grdrsis.f90
vgen/mt_tofrom_grid.F90
vgen/int_nv.F90 vgen/int_nv.F90
vgen/lhglptg.f90 vgen/lhglptg.f90
vgen/lhglpts.f90 vgen/lhglpts.f90
...@@ -32,12 +33,11 @@ vgen/prp_xcfft_map.f90 ...@@ -32,12 +33,11 @@ vgen/prp_xcfft_map.f90
vgen/psqpw.F90 vgen/psqpw.F90
vgen/rotate_int_den_to_local.F90 vgen/rotate_int_den_to_local.F90
vgen/vintcz.f90 vgen/vintcz.f90
vgen/visxc.f90 vgen/vis_xc.F90
vgen/visxcg.f90 vgen/pw_tofrom_grid.F90
vgen/vmatgen.f90 vgen/vmatgen.f90
vgen/vmts.F90 vgen/vmts.F90
vgen/vmtxc.f90 vgen/vmt_xc.F90
vgen/vmtxcg.F90
vgen/vvac.f90 vgen/vvac.f90
vgen/vvacis.f90 vgen/vvacis.f90
vgen/vvacxc.f90 vgen/vvacxc.f90
......
This diff is collapsed.
...@@ -3,75 +3,85 @@ ...@@ -3,75 +3,85 @@
! This file is part of FLEUR and available as free software under the conditions ! 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. ! of the MIT license as expressed in the LICENSE file in more detail.
!-------------------------------------------------------------------------------- !--------------------------------------------------------------------------------
MODULE m_mkgxyz3 MODULE m_mkgxyz3
c.....------------------------------------------------------------------ !.....------------------------------------------------------------------
c by use of cartesian x,y,z components of charge density gradients, !c by use of cartesian x,y,z components of charge density gradients,
c make the quantities !c make the quantities
cc agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr !cc agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
cc used to calculate gradient contribution to xc potential and !cc used to calculate gradient contribution to xc potential and
cc energy. !cc energy.
c.....------------------------------------------------------------------ !c.....------------------------------------------------------------------
CONTAINS CONTAINS
SUBROUTINE mkgxyz3( SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,grad)
> ndm,jsdm,ng3,jspins,vl, USE m_types
> dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy, IMPLICIT NONE
< grad) REAL, INTENT (IN) :: vl(:,:)
USE m_types REAL, INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
IMPLICIT NONE REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
INTEGER, INTENT (IN) :: ndm,ng3,jsdm,jspins REAL, INTENT (IN) :: dvyz(:,:),dvzx(:,:),dvxy(:,:)
REAL, INTENT (IN) :: vl(ndm,jsdm)
REAL, INTENT (IN) :: dvx(ndm,jsdm),dvy(ndm,jsdm),dvz(ndm,jsdm) TYPE(t_gradients),INTENT(INOUT)::grad
REAL, INTENT (IN) :: dvxx(ndm,jsdm),dvyy(ndm,jsdm),dvzz(ndm,jsdm)
REAL, INTENT (IN) :: dvyz(ndm,jsdm),dvzx(ndm,jsdm),dvxy(ndm,jsdm) REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,&
vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,&
TYPE(t_gradients),INTENT(INOUT)::grad vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,&
dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt, dagrzu,dzdx,dzdy,dzdz,sml
& vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu, INTEGER i,js,jspins,nsp
& vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,
& dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd, nsp=SIZE(dvx,1)
+ dagrzu,dzdx,dzdy,dzdz, jspins=SIZE(dvx,2)
+ sml sml = 1.e-14
INTEGER i
IF (ALLOCATED(grad%gr)) THEN
sml = 1.e-14 ! Gradients for libxc
DO js=1,jspins
IF(ALLOCATED(grad%sigma)) THEN DO i=1,nsp
!Use only contracted gradients for libxc grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
if (jspins==1) THEN ENDDO
DO i=1,ng3 END DO
grad%sigma(1,i)= IF(ALLOCATED(grad%sigma)) THEN
+ dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1) !Use only contracted gradients for libxc
ENDDO IF (jspins==1) THEN
ELSE DO i=1,nsp
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(1,i)= ENDDO
+ dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1) ELSE
grad%sigma(2,i)= DO i=1,nsp
+ dvx(i,1)*dvx(i,2)+dvy(i,1)*dvy(i,2)+dvz(i,1)*dvz(i,2) grad%sigma(1,i)= dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
grad%sigma(3,i)= grad%sigma(2,i)= dvx(i,1)*dvx(i,2)+dvy(i,1)*dvy(i,2)+dvz(i,1)*dvz(i,2)
+ dvx(i,2)*dvx(i,2)+dvy(i,2)*dvy(i,2)+dvz(i,2)*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 ENDDO
ENDIF ENDIF
RETURN END IF
ENDIF IF(ALLOCATED(grad%laplace)) THEN
DO js=1,jspins
DO i = 1,ndm DO i=1,nsp
grad%agrt(i) = 0.0 grad%laplace(i,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
grad%agru(i) = 0.0 ENDDO
grad%agrd(i) = 0.0 ENDDO
grad%gggrt(i) = 0.0 ENDIF
grad%gggru(i) = 0.0 RETURN
grad%gggrd(i) = 0.0 ENDIF
grad%gzgr(i) = 0.0
grad%g2rt(i) = 0.0 IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
grad%g2ru(i) = 0.0
grad%g2rd(i) = 0.0 DO i = 1,size(grad%agrt)
ENDDO grad%agrt(i) = 0.0
grad%agru(i) = 0.0
IF (jspins.eq.1) THEN grad%agrd(i) = 0.0
grad%gggrt(i) = 0.0
DO 10 i = 1,ng3 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) vlu=max(vl(i,1)/2,sml)
dvxu=dvx(i,1)/2 dvxu=dvx(i,1)/2
...@@ -108,7 +118,7 @@ c.....------------------------------------------------------------------ ...@@ -108,7 +118,7 @@ c.....------------------------------------------------------------------
dvzxt = dvzxu + dvzxd dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd 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%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**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. ...@@ -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%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd 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 dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2 dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/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 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%g2rt(i) = dvxxt + dvyyt + dvzzt
grad%g2ru(i) = dvxxu + dvyyu + dvzzu grad%g2ru(i) = dvxxu + dvyyu + dvzzu
grad%g2rd(i) = dvxxd + dvyyd + dvzzd 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) vlu = max(vl(i,1),sml)
dvxu=dvx(i,1) dvxu=dvx(i,1)
...@@ -187,7 +197,7 @@ c g2r: grad(grad(ro)) ...@@ -187,7 +197,7 @@ c g2r: grad(grad(ro))
dvzxt = dvzxu + dvzxd dvzxt = dvzxu + dvzxd
dvxyt = dvxyu + dvxyd 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%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml)
grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**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. ...@@ -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%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd