Dear Fabio,
you are right, I forgot to switch to the f03 version the call to the mgga
functional routine.
I do not guarantee anything about the results, but now you can properly set
the c parameter of TB09 as I told in the previous email. At line 24 I have
put zero as default value for c, but with no particular reason.
I don't know how old is you version of QE, but this file has been taken and
modified from the latest version in the develop branch.
Hope it helps,
Fabrizio
On Fri, Nov 8, 2019 at 5:44 PM Fabio Costa <fabiocos...@hotmail.com> wrote:
> Dear Fabrizio
>
> Thank you for your quick reply
>
> I followed your advice, and now I'm obtaining strange results. I'm running
> some tests with silicon, and to begin with, the scf run with mBJ reports a
> negative gap. In the pw.x output, the highest occupied level is higher in
> energy than the lowest unoccupied level, resulting on a "gap" of -0.972 eV.
> Just to compare, the same calculation done with PBE results on an expected
> result of ~0.8 eV gap. Going further, I plotted the density of states for
> both cases, and for my surprise, the result obtained with mBJ has no gap at
> all, as can be seen from the attached file.
>
> Also, I'm not sure if I proceeded in the right way, but I tried to change
> the value of the parameter "c" by changing it at line 24 from
> xc_mgga_drivers.f90, and recompiling pw.x . After doing this, the result
> was exactly the same, so my guess is that maybe I misunderstood how to
> change it.
>
> Sorry for the trouble, and thank you for your patience
> Fábio Costa
> ------------------------------
> *De:* users <users-boun...@lists.quantum-espresso.org> em nome de
> Fabrizio Ferrari <ferrariruffino...@gmail.com>
> *Enviado:* sexta-feira, 8 de novembro de 2019 11:51
> *Para:* Quantum ESPRESSO users Forum <users@lists.quantum-espresso.org>
> *Assunto:* Re: [QE-users] On the use of the modified Becke-Johnson (TB09)
> functional
>
> Dear Fabio,
> I added some lines in the attached file 'Modules/xc_mgga_drivers.f90' so
> that, by replacing it with the original one, you should be able to set the
> c parameter in tb09. You can do it directly into that module by setting the
> initial value of cc_param or in any part of the program by calling the
> routine 'set_cc_for_tb09(xx)' (but in case you may need to adjust the
> module dependencies).
> Before the compilation you need to add the '-lxcf03' in LD_LIBS in the
> make.inc file
> Hope it helps,
> Fabrizio
>
> On Wed, Nov 6, 2019 at 8:36 PM Fabio Costa <fabiocos...@hotmail.com>
> wrote:
>
> Dear users and developers
>
> After some more attempts to use the TB09 functional, I ended up finding
> this previous post:
> https://lists.quantum-espresso.org/pipermail/users/2019-April/042619.html ,
> and after changing v_of_rho.f90, the code is now able to end without any
> error.
>
> In my attempts, I'm performing the calculations as an initial scf run with
> a PBE pseudopotencial, and then another scf run, now with restart_mode =
> 'restart' and input_dft = 'MGGA_X_TB09'. By following this procedure, as
> described in E. Germaneau's paper, the calculations do converge.
>
> The issue now is that the gap energies I'm finding are very large. For
> example, for silicon, the calculations result in a 2.8 eV gap, larger even
> than the results presented in the paper, where the values ranged between
> 1.9 to 2.3. Despite the large gaps, the band structure and DOS agree well
> with the ones obtained with PBE.
>
> Is there any way to improve these results. My first guess is to adjust the
> "c" parameter in the mBJ functional. After some searching through the
> routines and Libxc files, I'm still clueless on how to proceed on this
> matter.
>
> Any further suggestions will be much appreciated
>
> Thank you
> Fábio Costa
>
>
>
> ------------------------------
> *De:* users <users-boun...@lists.quantum-espresso.org> em nome de Fabio
> Costa <fabiocos...@hotmail.com>
> *Enviado:* segunda-feira, 14 de outubro de 2019 13:10
> *Para:* users@lists.quantum-espresso.org <users@lists.quantum-espresso.org
> >
> *Assunto:* [QE-users] On the use of the modified Becke-Johnson (TB09)
> functional
>
> Dear all
>
> I'm having some difficulties to work with the TB09 functional. After some
> reading, I found in some old posts that this can be done by compiling QE
> with Libxc, and by adding the line 'input_dft='tb09'' in the pw.x input.
>
> When I try to use this functional, the scf calculation always ends like
> this:
>
> total energy = NaN Ry Harris-Foulkes
> estimate = -11.24358313 Ry estimated scf accuracy <
> 0.01858366 Ry iteration # 18 ecut= 120.00 Ry beta= 0.10
> Davidson diagonalization with
> overlap
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> Error in routine cdiaghg (7): eigenvectors failed to
> converge
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> What draws my attention is that always in the last cycle before the error
> the total energy is printed as NaN. The parameters to the calculation seem
> to be fine, as it ends smoothly without the input_dft line.
>
> By searching in the QE user guide, there is a reference to the work of E.
> Germaneau, which lead me to his paper (DOI: 10.1016/j.cpc.2013.02.020
> <http://dx.doi.org/10.1016/j.cpc.2013.02.020>), where the results of his
> implementation are presented. Even though, it is unclear to me how to use
> the functional correctly.
>
> Thanks for any assistance
> Fábio Costa
> _______________________________________________
> Quantum ESPRESSO is supported by MaX (www.max-centre.eu/quantum-espresso)
> users mailing list users@lists.quantum-espresso.org
> https://lists.quantum-espresso.org/mailman/listinfo/users
>
>
MODULE xc_mgga
!
USE kinds, ONLY: DP
USE funct, ONLY: get_meta, get_metac, is_libxc, &
exx_is_active, scan_exx, get_exx_fraction
!
IMPLICIT NONE
!
PRIVATE
SAVE
!
! GGA exchange-correlation drivers
PUBLIC :: xc_metagcx
PUBLIC :: tau_xc, tau_xc_spin
PUBLIC :: change_threshold_mgga, set_cc_for_tb09
!
!
! input thresholds (default values)
REAL(DP) :: rho_threshold = 1.0E-8_DP
REAL(DP) :: grho2_threshold = 1.0E-12_DP
REAL(DP) :: tau_threshold = 1.0E-8_DP
!
!TB09 ONLY********************************************+++
REAL(DP) :: cc_param = 0.0_DP
!***************************************************
!
CONTAINS
!
!
!-------------------------------------------------------------------------------------
SUBROUTINE change_threshold_mgga( rho_thr_in, grho2_thr_in, tau_thr_in )
!------------------------------------------------------------------------------------
!! Change rho, grho and tau thresholds.
!
REAL(DP), INTENT(IN) :: rho_thr_in
REAL(DP), INTENT(IN), OPTIONAL :: grho2_thr_in
REAL(DP), INTENT(IN), OPTIONAL :: tau_thr_in
!
rho_threshold = rho_thr_in
IF ( PRESENT(grho2_thr_in) ) grho2_threshold = grho2_thr_in
IF ( PRESENT(tau_thr_in) ) tau_threshold = tau_thr_in
!
RETURN
!
END SUBROUTINE change_threshold_mgga
!
!
! ONLY FOR TB09 ********************************************************
SUBROUTINE set_cc_for_tb09( cc_in )
!--------------------------------------------------------------------
!! Change rho threshold.
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: cc_in
!
cc_param = cc_in
!
RETURN
!
END SUBROUTINE
!*********************************************************************
!
!
!----------------------------------------------------------------------------------------
SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!-------------------------------------------------------------------------------------
!! Wrapper routine. Calls metaGGA drivers from internal libraries
!! of q-e or from the external libxc, depending on the input choice.
!
#if defined(__LIBXC)
USE funct, ONLY : get_libxc_flags_exc
USE xc_f90_types_m
USE xc_f90_lib_m
USE xc_f03_lib_m
#endif
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!! length of the I/O arrays
INTEGER, INTENT(IN) :: ns
!! spin components
INTEGER, INTENT(IN) :: np
!! first dimension of v2c
REAL(DP), INTENT(IN) :: rho(length,ns)
!! the charge density
REAL(DP), INTENT(IN) :: grho(3,length,ns)
!! grho = \nabla rho
REAL(DP), INTENT(IN) :: tau(length,ns)
!! kinetic energy density
REAL(DP), INTENT(OUT) :: ex(length)
!! sx = E_x(rho,grho)
REAL(DP), INTENT(OUT) :: ec(length)
!! sc = E_c(rho,grho)
REAL(DP), INTENT(OUT) :: v1x(length,ns)
!! v1x = D(E_x)/D(rho)
REAL(DP), INTENT(OUT) :: v2x(length,ns)
!! v2x = D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
REAL(DP), INTENT(OUT) :: v3x(length,ns)
!! v3x = D(E_x)/D(tau)
REAL(DP), INTENT(OUT) :: v1c(length,ns)
!! v1c = D(E_c)/D(rho)
REAL(DP), INTENT(OUT) :: v2c(np,length,ns)
!! v2c = D(E_c)/D( D rho/D r_alpha ) / |\nabla rho|
REAL(DP), INTENT(OUT) :: v3c(length,ns)
!! v3c = D(E_c)/D(tau)
!
! ... local variables
!
INTEGER :: is, imeta, imetac
REAL(DP), ALLOCATABLE :: grho2(:,:)
!
#if defined(__LIBXC)
TYPE(xc_f90_pointer_t) :: xc_func
TYPE(xc_f90_pointer_t) :: xc_info1, xc_info2
!
!TB09 ONLY******************************************************************
TYPE(xc_f03_func_t) :: xc_func03
REAL(DP) :: ext_params(1) !number of external parameters could be bigger,
!depending on the chosen functionals
!*********************************************************************
REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:), tau_lxc(:)
REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:), vx_tau(:)
REAL(DP), ALLOCATABLE :: vc_rho(:), vc_sigma(:), vc_tau(:)
REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in TPSS
!
REAL(DP) :: exx_fraction
INTEGER :: k, ipol, pol_unpol, eflag
LOGICAL :: POLARIZED
!
imeta = get_meta()
imetac = get_metac()
!
ex = 0.0_DP ; v1x = 0.0_DP ; v2x = 0.0_DP ; v3x = 0.0_DP
ec = 0.0_DP ; v1c = 0.0_DP ; v2c = 0.0_DP ; v3c = 0.0_DP
!
POLARIZED = .FALSE.
IF (ns == 2) THEN
POLARIZED = .TRUE.
ENDIF
!
pol_unpol = ns
!
ALLOCATE( rho_lxc(length*ns), sigma(length*ns), tau_lxc(length*ns) )
ALLOCATE( lapl_rho(length*ns) )
!
ALLOCATE( ex_lxc(length) , ec_lxc(length) )
ALLOCATE( vx_rho(length*ns) , vx_sigma(length*ns), vx_tau(length*ns) )
ALLOCATE( vc_rho(length*ns) , vc_sigma(length*ns), vc_tau(length*ns) )
ALLOCATE( vlapl_rho(length*ns) )
!
!
IF ( ns == 1 ) THEN
!
DO k = 1, length
rho_lxc(k) = ABS( rho(k,1) )
IF ( rho_lxc(k) > rho_threshold ) &
sigma(k) = grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2
tau_lxc(k) = tau(k,1)
ENDDO
!
ELSE
!
DO k = 1, length
rho_lxc(2*k-1) = rho(k,1)
rho_lxc(2*k) = rho(k,2)
!
sigma(3*k-2) = grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2
sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) + &
grho(3,k,1) * grho(3,k,2)
sigma(3*k) = grho(1,k,2)**2 + grho(2,k,2)**2 + grho(3,k,2)**2
!
tau_lxc(2*k-1) = tau(k,1)
tau_lxc(2*k) = tau(k,2)
ENDDO
!
ENDIF
!
IF ( .NOT.is_libxc(5) .AND. imetac==0 ) THEN
!
ALLOCATE( grho2(length,ns) )
!
DO is = 1, ns
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
IF (ns == 1) THEN
CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
ELSEIF (ns == 2) THEN
CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, &
v2c, v3c )
ENDIF
!
DEALLOCATE( grho2 )
!
ENDIF
!
! META EXCHANGE
!
IF ( is_libxc(5) ) THEN
!
!ONLY FOR TB09***********************************************************
IF ( imeta==208 ) THEN
CALL xc_f03_func_init( xc_func03, imeta, 1 )
ext_params(1) = cc_param
CALL xc_f03_func_set_ext_params( xc_func03, ext_params )
CALL xc_f03_mgga_vxc( xc_func03, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
CALL xc_f03_func_end( xc_func03 )
ELSE
!********************************************************************
!
CALL xc_f90_func_init( xc_func, xc_info1, imeta, pol_unpol )
CALL get_libxc_flags_exc( xc_info1, eflag )
IF (eflag==1) THEN
CALL xc_f90_mgga_exc_vxc( xc_func, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ex_lxc(1), vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
ELSE
CALL xc_f90_mgga_vxc( xc_func, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
ENDIF
CALL xc_f90_func_end( xc_func )
!**********************************************************
ENDIF
!**********************************************************
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
ex(k) = ex_lxc(k) * rho_lxc(k)
v1x(k,1) = vx_rho(k)
v2x(k,1) = vx_sigma(k) * 2.0_DP
v3x(k,1) = vx_tau(k)
ENDDO
ELSE
DO k = 1, length
ex(k) = ex_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1x(k,1) = vx_rho(2*k-1)
v1x(k,2) = vx_rho(2*k)
v2x(k,1) = vx_sigma(3*k-2)*2.d0
v2x(k,2) = vx_sigma(3*k)*2.d0
v3x(k,1) = vx_tau(2*k-1)
v3x(k,2) = vx_tau(2*k)
ENDDO
ENDIF
!
! ... only for HK/MCA: SCAN0 (used in CPV)
IF ( scan_exx ) THEN
exx_fraction = get_exx_fraction()
IF (exx_is_active()) THEN
ex = (1.0_DP - exx_fraction) * ex
v1x = (1.0_DP - exx_fraction) * v1x
v2x = (1.0_DP - exx_fraction) * v2x
v3x = (1.0_DP - exx_fraction) * v3x
ENDIF
ENDIF
!
ENDIF
!
! META CORRELATION
!
IF ( is_libxc(6) ) THEN
!
CALL xc_f90_func_init( xc_func, xc_info1, imetac, pol_unpol )
CALL xc_f90_mgga_exc_vxc( xc_func, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ec_lxc(1), vc_rho(1), vc_sigma(1), vlapl_rho(1), vc_tau(1) )
CALL xc_f90_func_end( xc_func )
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
ec(k) = ec_lxc(k) * rho_lxc(k) !* SIGN(1.0_DP, rho(k,1))
v1c(k,1) = vc_rho(k)
v2c(1,k,1) = vc_sigma(k) * 2.0_DP
v3c(k,1) = vc_tau(k)
ENDDO
ELSE
DO k = 1, length
ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1c(k,1) = vc_rho(2*k-1)
v1c(k,2) = vc_rho(2*k)
DO ipol = 1, 3
v2c(ipol,k,1) = vc_sigma(3*k-2)*grho(ipol,k,1)*2.D0 + vc_sigma(3*k-1)*grho(ipol,k,2)
v2c(ipol,k,2) = vc_sigma(3*k) *grho(ipol,k,2)*2.D0 + vc_sigma(3*k-1)*grho(ipol,k,1)
ENDDO
v3c(k,1) = vc_tau(2*k-1)
v3c(k,2) = vc_tau(2*k)
ENDDO
ENDIF
!
ENDIF
!
DEALLOCATE( rho_lxc, sigma, tau_lxc, lapl_rho )
DEALLOCATE( ex_lxc , ec_lxc )
DEALLOCATE( vx_rho , vx_sigma, vx_tau )
DEALLOCATE( vc_rho , vc_sigma, vc_tau, vlapl_rho )
!
#else
!
ALLOCATE( grho2(length,ns) )
!
DO is = 1, ns
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
IF (ns == 1) THEN
!
CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
!
ELSEIF (ns == 2) THEN
!
CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, &
v2c, v3c )
!
ENDIF
!
DEALLOCATE( grho2 )
!
#endif
!
RETURN
!
END SUBROUTINE xc_metagcx
!
!
!---------------------------------------------------------------------------------
SUBROUTINE tau_xc( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!-------------------------------------------------------------------------------
! gradient corrections for exchange and correlation - Hartree a.u.
! See comments at the beginning of module for implemented cases
!
! input: rho, grho=|\nabla rho|^2
!
! definition: E_x = \int e_x(rho,grho) dr
!
! output: sx = e_x(rho,grho) = grad corr
! v1x= D(E_x)/D(rho)
! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
! v3x= D(E_x)/D(tau)
!
! sc, v1c, v2c as above for correlation
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!
INTEGER :: k, imeta
REAL(DP) :: arho
REAL(DP), DIMENSION(length) :: rho, grho2, tau, &
ex, ec, v1x, v2x, v3x, v1c, v2c, v3c
!
imeta = get_meta()
!
v1x=0.d0 ; v2x=0.d0 ; v3x=0.d0 ; ex=0.d0
v1c=0.d0 ; v2c=0.d0 ; v3c=0.d0 ; ec=0.d0
!
DO k = 1, length
!
arho = ABS(rho(k))
!
IF ( (arho<=rho_threshold).OR.(grho2(k)<=grho2_threshold).OR.(ABS(tau(k))<=rho_threshold) ) CYCLE
!
SELECT CASE( imeta )
CASE( 1 )
CALL tpsscxc( arho, grho2(k), tau(k), ex(k), ec(k), v1x(k), v2x(k), v3x(k), v1c(k), v2c(k), v3c(k) )
CASE( 2 )
CALL m06lxc( arho, grho2(k), tau(k), ex(k), ec(k), v1x(k), v2x(k), v3x(k), v1c(k), v2c(k), v3c(k) )
CASE DEFAULT
CALL errore( 'tau_xc', 'This case is not implemented', imeta )
END SELECT
!
ENDDO
!
RETURN
!
END SUBROUTINE tau_xc
!
!
!----------------------------------------------------------------------------------------
SUBROUTINE tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!------------------------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN) :: rho(length,2), tau(length,2)
REAL(DP), INTENT(IN) :: grho(3,length,2)
!
REAL(DP), INTENT(OUT) :: ex(length), ec(length), v1x(length,2), v2x(length,2), &
v3x(length,2), v1c(length,2), v3c(length,2)
REAL(DP), INTENT(OUT) :: v2c(3,length,2)
!
! ... local variables
!
INTEGER :: k, ipol, imeta
REAL(DP) :: rh, zeta, atau, grho2(2), ggrho2
REAL(DP) :: v2cup, v2cdw
!
imeta = get_meta()
!
ex=0.0_DP ; v1x=0.0_DP ; v2x=0.0_DP ; v3x=0.0_DP
ec=0.0_DP ; v1c=0.0_DP ; v2c=0.0_DP ; v3c=0.0_DP
!
! FIXME: for SCAN, this will be calculated later
!
DO k = 1, length
!
rh = rho(k,1) + rho(k,2)
atau = tau(k,1) + tau(k,2) ! KE-density in Hartree
grho2(1) = SUM( grho(:,k,1)**2 )
grho2(2) = SUM( grho(:,k,2)**2 )
ggrho2 = ( grho2(1) + grho2(2) ) * 4.0_DP
!
IF ((rh <= rho_threshold) .OR. (ggrho2 <= grho2_threshold) .OR. (ABS(atau) <= tau_threshold)) CYCLE
!
SELECT CASE( imeta )
CASE( 1 )
!
CALL tpsscx_spin( rho(k,1), rho(k,2), grho2(1), grho2(2), tau(k,1), &
tau(k,2), ex(k), v1x(k,1), v1x(k,2), v2x(k,1), v2x(k,2), v3x(k,1), v3x(k,2) )
!
zeta = (rho(k,1) - rho(k,2)) / rh
!
CALL tpsscc_spin( rh, zeta, grho(:,k,1), grho(:,k,2), atau, ec(k), &
v1c(k,1), v1c(k,2), v2c(:,k,1), v2c(:,k,2), v3c(k,1), v3c(k,2) )
!
CASE( 2 )
!
CALL m06lxc_spin( rho(k,1), rho(k,2), grho2(1), grho2(2), tau(k,1), tau(k,2), ex(k), ec(k), &
v1x(k,1), v1x(k,2), v2x(k,1), v2x(k,2), v3x(k,1), v3x(k,2), &
v1c(k,1), v1c(k,2), v2cup , v2cdw , v3c(k,1), v3c(k,2) )
!
v2c(:,k,1) = v2cup*grho(:,k,1)
v2c(:,k,2) = v2cdw*grho(:,k,2)
!
CASE DEFAULT
!
CALL errore( 'tau_xc_spin', 'This case not implemented', imeta )
!
END SELECT
!
ENDDO
!
RETURN
!
END SUBROUTINE tau_xc_spin
!
!
END MODULE xc_mgga
_______________________________________________
Quantum ESPRESSO is supported by MaX (www.max-centre.eu/quantum-espresso)
users mailing list users@lists.quantum-espresso.org
https://lists.quantum-espresso.org/mailman/listinfo/users
