Greetings,
Sorry for the long delay in responding...
I tried setting the values manually for a small set of c's (see data.dat in
attachment) with
call xc_f90_mgga_x_tb09_set_par(xc_func, cc)
and ran the silicon structure with the tpss pseudopotential in
the QE site.
Indeed, for c=1 I got the same value as the default calculation...
(And a rather large gap although that had already been mentioned in
another thread http://qe-forge.org/pipermail/pw_forum/2015-April/106571.html)
I also tried computing the c parameter as per eq. 3 (see attach. cc_tb09.f90).
For a (default) calculation with the same pseudo for Si I get cc = 1.0135.
From the cc runs however, the 'optimal' value should be around 0.85 ~ 1.00,
so there is a certain discrepancy with the original article...
Just out of curiosity, is there a specific reason for the usage of TPSS in
correlation instead of the original PW92? The former seems to be rather
incompatible with non-tpss pseudos...
Quoting Paolo Giannozzi <[email protected]>:
Hard to say. Since you already had a look at the problem, you are in the
best position to find whether the code computes what it is expected to
compute or not. Try to set c to 1, then to a different parameter, see what
happens.
Paolo
On Wed, Feb 24, 2016 at 4:05 PM, Pedro Miguel Castro Borlido <
[email protected]> wrote:
Greetings,
I understand that the implementation of the TB09 meta gga is done
in Espresso
though Libxc, as described in the article by Germaneau and colleagues.
However, while looking through the code it seems that at no point
the c parameter
(equation 3 of PRL 102, 226401 (2009)) is set with
xc_f90_mgga_x_tb09_set_par thus
leaving c equal to 1 by default (i.e. the original BJ is used).
Could someone please check if this is actually a bug of if c is
actually set?
Best regards,
Pedro Borlido
Friedrich-Schiller Universität Jena
Institut für Festkörpertheorie und -optik
Max-Wien-Platz 1
07743 Jena
_______________________________________________
Pw_forum mailing list
[email protected]
http://pwscf.org/mailman/listinfo/pw_forum
--
Paolo Giannozzi, Dip. Scienze Matematiche Informatiche e Fisiche,
Univ. Udine, via delle Scienze 208, 33100 Udine, Italy
Phone +39-0432-558216, fax +39-0432-558222
subroutine cc_tb09(rho, rho_core, rhog_core, cc)
USE kinds, ONLY : DP
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : omega, alat
USE gvect, ONLY : g, nl,ngm
USE scf, ONLY : scf_type
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum
USE io_global, ONLY : stdout
USE lsda_mod, ONLY : nspin
!
type(scf_type), intent(in) :: rho
real(DP), intent(in) :: rho_core(dfftp%nnr)
complex(DP), intent(in) :: rhog_core(ngm)
real(DP), intent(out) :: cc
!
! Local variables
integer :: is, ir, idx0, idx, i, j, k
real(dp) :: grho_abs
real(dp), allocatable :: rhoout(:), grho(:,:)
complex(dp), allocatable :: rhogsum(:)
!
allocate( rhoout(dfftp%nnr) )
allocate( rhogsum(ngm) )
allocate( grho(3, dfftp%nnr) )
!
! Get the total charge, including spins
! (We just use this routine for nspin=1 so not really relevant)
rhoout(:) = rho_core(:)
rhogsum(:) = rhog_core(:)
!
do is=1, nspin
rhoout(:) = rhoout(:) + rho%of_r(:,is)
rhogsum(:) = rhogsum(:) + rho%of_g(:, is)
end do
!
! Compute gradient
call gradrho( dfftp%nnr, rhogsum, ngm, g, nl, grho )
!
! Sum over real space
cc = 0.0_dp
!
#if defined(__MPI)
idx0 = dfftp%nr1x * dfftp%nr2x * dfftp%ipp(me_bgrp + 1)
#else
idx0 = 0
#endif
!
do ir = 1, dfftp%nnr
!
! ... Three dimentional indices
!
idx = idx0 + ir - 1
k = idx / (dfftp%nr1x*dfftp%nr2x)
idx = idx - (dfftp%nr1x*dfftp%nr2x)*k
j = idx / dfftp%nr1x
idx = idx - dfftp%nr1x*j
i = idx
!
! ... Do not include points outside physical range
!
if ( i >= dfftp%nr1 .or. j >= dfftp%nr2 .or. k >= dfftp%nr3 ) cycle
!
grho_abs = sqrt( grho(1, ir)**2 + grho(2, ir)**2 + grho(3, ir)**2 )
!
cc = cc + grho_abs/rhoout(ir)
!
end do
!
cc = cc / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
call mp_sum( cc, intra_bgrp_comm)
!
! ... Paramters from TB article
!
cc = -0.012 + 1.023*sqrt(cc)
!
deallocate( rhoout, rhogsum, grho )
end subroutine cc_tb09
#
# cc
# Final Energy (ry)
# HOMO eV
# LUMO eV
# Gap eV
#
# Default run: -11.34597955 3.8591 5.5980 1.7389
0.50 -11.59246838 6.4581 6.2345 -0.2236
0.70 -11.53709797 5.4177 5.8598 0.4421
0.80 -11.49099695 4.8963 5.7233 0.8270
0.85 -11.46199529 4.6358 5.6724 1.0366
0.90 -11.42850210 4.3760 5.6341 1.2581
1.00 -11.34597955 3.8591 5.5980 1.7389
1.10 -11.23796484 3.3463 5.6239 2.2776
1.20 -11.09816747 2.8445 5.7249 2.8804
&control
calculation = 'scf'
prefix='silicon',
pseudo_dir = '/work3_ifto/nu67qih/work_tran_blaha/pseudos.ncv',
outdir='./tmp',
verbosity = 'high',
wf_collect = .true.,
/
&system
ibrav= 2,
celldm(1) = 10.262630155830101,
nat = 2,
ntyp= 1,
smearing = 'methfessel-paxton',
!input_dft = 'tb09',
degauss=0.01,
ecutwfc = 70.0,
nbnd = 10,
/
&electrons
/
ATOMIC_SPECIES
Si 28.086 Si.tpss-mt.UPF
ATOMIC_POSITIONS
Si 0.00 0.00 0.00
Si 0.25 0.25 0.25
K_POINTS automatic
20 20 20 0 0 0
_______________________________________________
Pw_forum mailing list
[email protected]
http://pwscf.org/mailman/listinfo/pw_forum
