Commit 869a13a2 authored by Gregor Michalicek's avatar Gregor Michalicek

Bugfix for hybrid functionals in vgen/vgen_xcpot.F90

   ...vXC%pw_w was not initialized.

   + some code beautification
parent e10a5740
...@@ -19,9 +19,8 @@ CONTAINS ...@@ -19,9 +19,8 @@ CONTAINS
!! TE_VEFF: charge density-effective potential integral !! TE_VEFF: charge density-effective potential integral
!! 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, & SUBROUTINE vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,&
vacuum, sym, obsolete, cell, oneD, sliceplot, mpi, results, noco, & obsolete,cell,oneD,sliceplot,mpi,results,noco,den,vTot,vXC,vCoul)
den, vTot, vx, vCoul )
USE m_rotate_int_den_to_local USE m_rotate_int_den_to_local
USE m_bfield USE m_bfield
...@@ -52,7 +51,7 @@ CONTAINS ...@@ -52,7 +51,7 @@ 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,vx,vCoul TYPE(t_potden), INTENT(INOUT) :: vTot,vXC,vCoul
TYPE(t_potden) :: workden,denRot TYPE(t_potden) :: workden,denRot
...@@ -61,42 +60,38 @@ CONTAINS ...@@ -61,42 +60,38 @@ CONTAINS
CALL vTot%resetPotDen() CALL vTot%resetPotDen()
CALL vCoul%resetPotDen() CALL vCoul%resetPotDen()
CALL vx%resetPotDen() CALL vXC%resetPotDen()
ALLOCATE( vx%pw_w, vTot%pw_w, mold=vTot%pw ) ALLOCATE(vXC%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, & CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,workden,vCoul,results)
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
CALL denRot%init( stars, atoms, sphhar, vacuum, input%jspins, noco%l_noco, 0 ) CALL denRot%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,0)
denRot=den denRot=den
CALL rotate_int_den_to_local( DIMENSION, sym, stars, atoms, sphhar, vacuum, cell, input, & CALL rotate_int_den_to_local(dimension,sym,stars,atoms,sphhar,vacuum,cell,input,noco,oneD,denRot)
noco, oneD, denRot )
ENDIF ENDIF
call vgen_xcpot( hybrid, input, xcpot, DIMENSION, atoms, sphhar, stars, & CALL vgen_xcpot(hybrid,input,xcpot,dimension,atoms,sphhar,stars,vacuum,sym,&
vacuum, sym, obsolete, cell, oneD, sliceplot, mpi, noco, den, denRot, vTot, vx, results ) obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vXC,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, vx%pw_w ) DEALLOCATE(vcoul%pw_w,vXC%pw_w)
CALL bfield(input,noco,atoms,field,vTot)
CALL bfield( input, noco, atoms, field, vTot )
! broadcast potentials
#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, vx ) CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,vXC)
#endif #endif
END SUBROUTINE vgen END SUBROUTINE vgen
......
...@@ -4,10 +4,14 @@ ...@@ -4,10 +4,14 @@
! 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_vgen_xcpot MODULE m_vgen_xcpot
USE m_juDFT USE m_juDFT
CONTAINS CONTAINS
SUBROUTINE vgen_xcpot(hybrid,input,xcpot,DIMENSION, atoms,sphhar,stars,&
vacuum,sym, obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vx,results) SUBROUTINE vgen_xcpot(hybrid,input,xcpot,dimension,atoms,sphhar,stars,vacuum,sym,&
obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,vTot,vXC,results)
! *********************************************************** ! ***********************************************************
! FLAPW potential generator * ! FLAPW potential generator *
! *********************************************************** ! ***********************************************************
...@@ -17,6 +21,8 @@ CONTAINS ...@@ -17,6 +21,8 @@ CONTAINS
! TE_VEFF: charge density-effective potential integral ! TE_VEFF: charge density-effective potential integral
! TE_EXC : charge density-ex-corr.energy density integral ! TE_EXC : charge density-ex-corr.energy density integral
! *********************************************************** ! ***********************************************************
USE m_types
USE m_constants USE m_constants
USE m_intnv USE m_intnv
USE m_vmtxcg USE m_vmtxcg
...@@ -26,37 +32,34 @@ CONTAINS ...@@ -26,37 +32,34 @@ CONTAINS
USE m_visxc USE m_visxc
USE m_visxcg USE m_visxcg
USE m_checkdopall USE m_checkdopall
USE m_cdn_io USE m_cdn_io
USE m_types USE m_convol
IMPLICIT NONE IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_hybrid),INTENT(IN) :: hybrid CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_obsolete),INTENT(IN) :: obsolete TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_sliceplot),INTENT(IN) :: sliceplot TYPE(t_obsolete), INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input TYPE(t_sliceplot), INTENT(IN) :: sliceplot
TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_sym),INTENT(IN) :: sym TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_potden), INTENT(IN) :: den, denRot TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_potden),INTENT(INOUT) :: vTot,vx TYPE(t_potden), INTENT(IN) :: den,denRot
TYPE(t_results),INTENT(INOUT),OPTIONAL :: results TYPE(t_potden), INTENT(INOUT) :: vTot,vXC
! .. TYPE(t_results), INTENT(INOUT), OPTIONAL :: results
! Local type instances
! .. Local type instances .. TYPE(t_potden) :: workDen,exc,veff
TYPE(t_potden) :: workDen,exc,veff ! Local Scalars
! .. Local Scalars .. INTEGER ifftd,ifftd2,ifftxc3d,ispin
INTEGER i,i3,irec2,irec3,ivac,j,js,k,k3,lh,n,nzst1
INTEGER ifftd,ifftd2, ifftxc3d
INTEGER jsp,l
#ifdef CPP_MPI #ifdef CPP_MPI
include 'mpif.h' include 'mpif.h'
integer:: ierr integer:: ierr
...@@ -69,13 +72,11 @@ CONTAINS ...@@ -69,13 +72,11 @@ CONTAINS
ALLOCATE(veff%pw_w,mold=veff%pw) ALLOCATE(veff%pw_w,mold=veff%pw)
ENDIF ENDIF
! ******** exchange correlation potential****************** ! exchange correlation potential
! vacuum region
! ---> vacuum region
IF (mpi%irank == 0) THEN IF (mpi%irank == 0) THEN
IF (input%film) THEN IF (input%film) THEN
CALL timestart("Vxc in vacuum") CALL timestart("Vxc in vacuum")
ifftd2 = 9*stars%mx1*stars%mx2 ifftd2 = 9*stars%mx1*stars%mx2
...@@ -84,163 +85,115 @@ CONTAINS ...@@ -84,163 +85,115 @@ CONTAINS
IF (.NOT.xcpot%is_gga()) THEN ! LDA IF (.NOT.xcpot%is_gga()) THEN ! LDA
IF (.NOT.oneD%odi%d1) THEN IF (.NOT.oneD%odi%d1) THEN
CALL vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,Den,vTot,exc)
CALL vvacxc(ifftd2,stars,vacuum,xcpot,input,noco,Den,&
vTot,exc)
ELSE ELSE
CALL judft_error("OneD broken") CALL judft_error("OneD broken")
! CALL vvacxc(& ! CALL vvacxc(stars,oneD%M,vacuum,odi%n2d,dimension,ifftd2,&
! & stars,oneD%M,vacuum,odi%n2d,dimension,ifftd2,& ! xcpot,input,odi%nq2,odi%nst2,den,noco,odi%kimax2%igf,&
! & xcpot,input,odi%nq2,& ! odl%pgf,vTot%vacxy,vTot%vacz,excxy,excz)
! & odi%nst2,den,noco,& END IF
! & odi%kimax2%igf,odl%pgf,&
! & vTot%vacxy,vTot%vacz,&
! & excxy,excz)
ENDIF
ELSE ! GGA ELSE ! GGA
IF (oneD%odi%d1) THEN IF (oneD%odi%d1) THEN
CALL judft_error("OneD broken") CALL judft_error("OneD broken")
! CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,&
!CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,& ! cell,xcpot,input,obsolete,workDen, ichsmrg,&
! cell,xcpot,input,obsolete,workDen, ichsmrg,& ! vTot%vacxy,vTot%vacz,rhmn, exc%vacxy,exc%vacz)
! vTot%vacxy,vTot%vacz,rhmn, exc%vacxy,exc%vacz)
ELSE ELSE
CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,& CALL vvacxcg(ifftd2,stars,vacuum,noco,oneD,cell,xcpot,input,obsolete,Den,vTot,exc)
cell,xcpot,input,obsolete,Den, &
vTot, exc)
END IF END IF
END IF END IF
CALL timestop("Vxc in vacuum") CALL timestop("Vxc in vacuum")
!+odim
END IF END IF
! ----------------------------------------
! ---> interstitial region ! interstitial region
CALL timestart("Vxc in interstitial") CALL timestart("Vxc in interstitial")
ifftd=27*stars%mx1*stars%mx2*stars%mx3 ifftd = 27*stars%mx1*stars%mx2*stars%mx3
IF ( (.NOT. obsolete%lwb) .OR. ( .not.xcpot%is_gga() ) ) THEN IF ((.NOT.obsolete%lwb).OR.(.NOT.xcpot%is_gga())) THEN
! no White-Bird-trick ! no White-Bird-trick
ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
IF ( .NOT.xcpot%is_gga() ) THEN IF (.NOT.xcpot%is_gga()) THEN ! LDA
! LDA CALL visxc(ifftd,stars,noco,xcpot,input,den,vTot,vXC,exc)
CALL visxc(ifftd,stars,noco,xcpot,input,den,vTot,vx,exc)
ELSE ! GGA ELSE ! GGA
CALL visxcg(ifftd,stars,sym,ifftxc3d,cell,den,xcpot,input,obsolete,noco,vTot,vXC,exc)
CALL visxcg(ifftd,stars,sym,ifftxc3d,cell,den,xcpot,input,&
obsolete,noco,vTot,vx,exc)
END IF END IF
ELSE ELSE
! White-Bird-trick ! White-Bird-trick
WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",& WRITE(6,'(a)') "W+B trick cancelled out. visxcwb uses at present common block cpgft3.",&
"visxcwb needs to be reprogrammed according to visxcg.f" "visxcwb needs to be reprogrammed according to visxcg.f"
CALL juDFT_error("visxcwb",calledby ="vgen") CALL juDFT_error("visxcwb",calledby ="vgen")
END IF END IF
!
CALL timestop("Vxc in interstitial") CALL timestop("Vxc in interstitial")
END IF !irank==0
ENDIF !irank==0 ! muffin tin spheres region
! IF (mpi%irank == 0) CALL timestart ("Vxc in MT")
! ------------------------------------------
! ----> muffin tin spheres region
IF (mpi%irank == 0) THEN
CALL timestart ("Vxc in MT")
END IF
#ifdef CPP_MPI #ifdef CPP_MPI
CALL MPI_BCAST(den%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr) CALL MPI_BCAST(den%mt,atoms%jmtd*(1+sphhar%nlhd)*atoms%ntype*dimension%jspd,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
#endif #endif
IF (xcpot%is_gga()) THEN IF (xcpot%is_gga()) THEN
CALL vmtxcg(dimension,mpi,sphhar,atoms, den,xcpot,input,sym,& CALL vmtxcg(dimension,mpi,sphhar,atoms,den,xcpot,input,sym,obsolete,vTot,vXC,exc)
obsolete, vTot,vx,exc)
ELSE ELSE
CALL vmtxc(DIMENSION,sphhar,atoms, den,xcpot,input,sym, vTot,exc,vx) CALL vmtxc(DIMENSION,sphhar,atoms,den,xcpot,input,sym,vTot,exc,vXC)
ENDIF ENDIF
!
! add MT EXX potential to vr ! add MT EXX potential to vr
!
IF (mpi%irank == 0) THEN IF (mpi%irank == 0) THEN
CALL timestop ("Vxc in MT") CALL timestop ("Vxc in MT")
! ------------------------------------------
! ---> check continuity of total potential
IF (input%vchk) THEN ! check continuity of total potential
CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,& IF (input%vchk) CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,vTot,1)
cell,vTot,1)
END IF
! ! TOTAL
!============TOTAL======================================
!
IF (PRESENT(results)) THEN IF (PRESENT(results)) THEN
! ! CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2
! CALCULATE THE INTEGRAL OF n1*Veff1 + n2*Veff2 ! Veff = Vcoulomb + Vxc
! Veff = Vcoulomb + Vxc
!
IF (noco%l_noco) THEN IF (noco%l_noco) THEN
workDen = denRot workDen = denRot
ELSE ELSE
workden=den workden = den
ENDIF END IF
veff=vTot veff = vTot
IF( xcpot%is_hybrid() ) THEN IF(xcpot%is_hybrid()) THEN
veff%pw = vTot%pw -xcpot%get_exchange_weight() * vx%pw DO ispin = 1, input%jspins
veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vx%pw_w CALL convol(stars,vXC%pw_w(:,ispin),vXC%pw(:,ispin),stars%ufft)
veff%mt = vTot%mt - xcpot%get_exchange_weight() * vx%mt END DO
veff%pw = vTot%pw - xcpot%get_exchange_weight() * vXC%pw
veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vXC%pw_w
veff%mt = vTot%mt - xcpot%get_exchange_weight() * vXC%mt
END IF END IF
results%te_veff = 0.0 results%te_veff = 0.0
DO js = 1,input%jspins DO ispin = 1, input%jspins
WRITE (6,FMT=8050) js WRITE (6,FMT=8050) ispin
8050 FORMAT (/,10x,'density-effective potential integrals for spin ',i2,/) 8050 FORMAT (/,10x,'density-effective potential integrals for spin ',i2,/)
CALL int_nv(js,stars,vacuum,atoms,sphhar, cell,sym,input,oneD,veff,workden,results%te_veff) CALL int_nv(ispin,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,veff,workden,results%te_veff)
ENDDO END DO
WRITE (6,FMT=8060) results%te_veff WRITE (6,FMT=8060) results%te_veff
8060 FORMAT (/,10x,'total density-effective potential integral :', t40,f20.10) 8060 FORMAT (/,10x,'total density-effective potential integral :', t40,f20.10)
!
! CALCULATE THE INTEGRAL OF n*exc ! CALCULATE THE INTEGRAL OF n*exc
!
! ---> perform spin summation of charge densities ! perform spin summation of charge densities for the calculation of Exc
! ---> for the calculation of Exc
CALL workden%sum_both_spin() CALL workden%sum_both_spin()
WRITE (6,FMT=8070) WRITE (6,FMT=8070)
8070 FORMAT (/,10x,'charge density-energy density integrals',/) 8070 FORMAT (/,10x,'charge density-energy density integrals',/)
results%te_exc = 0.0 results%te_exc = 0.0
CALL int_nv(1,stars,vacuum,atoms,sphhar, cell,sym,input,oneD,& CALL int_nv(1,stars,vacuum,atoms,sphhar,cell,sym,input,oneD,exc,workDen,results%te_exc)
exc,workDen, results%te_exc)
WRITE (6,FMT=8080) results%te_exc WRITE (6,FMT=8080) results%te_exc
8080 FORMAT (/,10x,'total charge density-energy density integral :', t40,f20.10) 8080 FORMAT (/,10x,'total charge density-energy density integral :', t40,f20.10)
END IF END IF
END IF ! mpi%irank == 0
ENDIF ! mpi%irank == 0
END SUBROUTINE vgen_xcpot END SUBROUTINE vgen_xcpot
END MODULE m_vgen_xcpot END MODULE m_vgen_xcpot
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment