From 2f933d101e6f36e02fde4911769db28e4f2ff03a Mon Sep 17 00:00:00 2001 From: fabrizio22 Date: Mon, 4 May 2020 13:51:13 +0200 Subject: [PATCH] XClib - lda - scratch --- CPV/src/Makefile | 2 +- EPW/src/Makefile | 2 +- GWW/bse/Makefile | 2 +- GWW/gww/Makefile | 2 +- GWW/head/Makefile | 2 +- GWW/pw4gww/Makefile | 2 +- GWW/simple/Makefile | 2 +- HP/src/Makefile | 2 +- Makefile | 7 +- Modules/correlation_gga.f90 | 3 +- Modules/correlation_lda_lsda.f90 | 469 +----------------------------- Modules/dmxc_drivers.f90 | 11 +- Modules/exchange_lda_lsda.f90 | 160 ---------- Modules/funct.f90 | 24 ++ Modules/make.depend | 7 + Modules/metagga.f90 | 5 +- Modules/xc_lda_lsda_drivers.f90 | 241 +--------------- Modules/xc_vdW_DF.f90 | 2 + NEB/src/Makefile | 3 +- PHonon/FD/Makefile | 2 +- PHonon/Gamma/Makefile | 2 +- PHonon/PH/Makefile | 2 +- PP/src/Makefile | 2 +- PP/src/make.depend | 19 ++ PP/src/ppacf.f90 | 3 +- PP/src/xc_vdW_scale_mod.f90 | 3 +- PW/src/Makefile | 2 +- PWCOND/src/Makefile | 2 +- TDDFPT/src/Makefile | 2 +- XClib/Makefile | 34 +++ XClib/constants_l.f90 | 155 ++++++++++ XClib/correlation_lda.f90 | 482 +++++++++++++++++++++++++++++++ XClib/dft_par_mod.f90 | 15 + XClib/exchange_lda.f90 | 176 +++++++++++ XClib/get_ldaxc.f90 | 59 ++++ XClib/kind_l.f90 | 19 ++ XClib/ldaxc_interfaces.f90 | 120 ++++++++ XClib/make.depend | 14 + XClib/xc_lda.f90 | 242 ++++++++++++++++ XSpectra/src/Makefile | 2 +- install/make.inc.in | 1 + install/makedeps.sh | 10 +- 42 files changed, 1419 insertions(+), 897 deletions(-) create mode 100644 XClib/Makefile create mode 100644 XClib/constants_l.f90 create mode 100644 XClib/correlation_lda.f90 create mode 100644 XClib/dft_par_mod.f90 create mode 100644 XClib/exchange_lda.f90 create mode 100644 XClib/get_ldaxc.f90 create mode 100644 XClib/kind_l.f90 create mode 100644 XClib/ldaxc_interfaces.f90 create mode 100644 XClib/make.depend create mode 100644 XClib/xc_lda.f90 diff --git a/CPV/src/Makefile b/CPV/src/Makefile index 861b8a266..9223892d0 100644 --- a/CPV/src/Makefile +++ b/CPV/src/Makefile @@ -114,7 +114,7 @@ makov_payne.o LOBJS = \ entropy.o -QEMODS=../../Modules/libqemod.a ../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a +QEMODS=../../Modules/libqemod.a ../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../XClib/libqeldaxc.a TLDEPS= bindir libs mods diff --git a/EPW/src/Makefile b/EPW/src/Makefile index c54db1ac6..e4ec57071 100644 --- a/EPW/src/Makefile +++ b/EPW/src/Makefile @@ -80,7 +80,7 @@ PWOBJS = ../../PW/src/libpw.a W90LIB = ../../wannier90-3.1.0/libwannier.a LRMODS = ../../LR_Modules/liblrmod.a PWOBJS = ../../PW/src/libpw.a -QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a \ +QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../dft-d3/libdftd3qe.a LIBOBJS =../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../clib/clib.a diff --git a/GWW/bse/Makefile b/GWW/bse/Makefile index b97221dc1..d5d2c2a3e 100644 --- a/GWW/bse/Makefile +++ b/GWW/bse/Makefile @@ -47,7 +47,7 @@ qpcorrections.o QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \ ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \ - ../../LAXlib/libqela.a ../../UtilXlib/libutil.a + ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../XClib/libqeldaxc.a PWOBJS = ../../PW/src/libpw.a PW4GWWOBJ = ../pw4gww/fft_custom.o ../pw4gww/stop_pp.o ../pw4gww/mp_wave_parallel.o GWWOBJ = ../gww/libgww.a ../minpack/minpacklib.a diff --git a/GWW/gww/Makefile b/GWW/gww/Makefile index 596364791..88f28f978 100644 --- a/GWW/gww/Makefile +++ b/GWW/gww/Makefile @@ -44,7 +44,7 @@ vcprim.o QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \ - ../../KS_Solvers/libks_solvers.a \ + ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../LAXlib/libqela.a ../../UtilXlib/libutil.a LIBMIN= ../minpack/minpacklib.a diff --git a/GWW/head/Makefile b/GWW/head/Makefile index b000d40c1..ce4b175dc 100644 --- a/GWW/head/Makefile +++ b/GWW/head/Makefile @@ -20,7 +20,7 @@ phq_readin.o \ solve_head.o QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a \ - ../../KS_Solvers/libks_solvers.a \ + ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \ ../../dft-d3/libdftd3qe.a LIBPWPH = ../../PHonon/PH/libph.a ../..//LR_Modules/liblrmod.a ../../PW/src/libpw.a diff --git a/GWW/pw4gww/Makefile b/GWW/pw4gww/Makefile index 97671c498..c4141bac1 100644 --- a/GWW/pw4gww/Makefile +++ b/GWW/pw4gww/Makefile @@ -50,7 +50,7 @@ operator_1_vp.o \ operator_debug.o QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a \ - ../../KS_Solvers/libks_solvers.a \ + ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a \ ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a # dft-d3 required by xlf for obscure reasons diff --git a/GWW/simple/Makefile b/GWW/simple/Makefile index 9579a8fb4..e412733ff 100644 --- a/GWW/simple/Makefile +++ b/GWW/simple/Makefile @@ -25,7 +25,7 @@ SIMPLEOBJS = \ QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \ ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \ ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \ - ../../dft-d3/libdftd3qe.a + ../../dft-d3/libdftd3qe.a ../../XClib/libqeldaxc.a # dft-d3 required by xlf for obscure reasons PWOBJS = ../../PW/src/libpw.a diff --git a/HP/src/Makefile b/HP/src/Makefile index 508c9dc5a..37b231c83 100644 --- a/HP/src/Makefile +++ b/HP/src/Makefile @@ -56,7 +56,7 @@ LRMODS = ../../LR_Modules/liblrmod.a PWOBJS = ../../PW/src/libpw.a QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a \ ../../upflib/libupf.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a\ - ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a + ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a ../../XClib/libqeldaxc.a TLDEPS= hplibs diff --git a/Makefile b/Makefile index 13987be69..8bb801774 100644 --- a/Makefile +++ b/Makefile @@ -175,7 +175,7 @@ pw4gwwlib : phlibs if test -d GWW ; then \ ( cd GWW ; $(MAKE) pw4gwwa || exit 1 ) ; fi -mods : libfox libutil libla libfft libupf libbeef libmbd +mods : libfox libutil libla libfft libupf libbeef libmbd librxc ( cd Modules ; $(MAKE) TLDEPS= all || exit 1 ) libks_solvers : libs libutil libla @@ -187,6 +187,9 @@ libla : liblapack libutil libcuda libfft : ( cd FFTXlib ; $(MAKE) TLDEPS= all || exit 1 ) +librxc : + ( cd XClib ; $(MAKE) TLDEPS= all || exit 1 ) + libutil : ( cd UtilXlib ; $(MAKE) TLDEPS= all || exit 1 ) @@ -277,7 +280,7 @@ install : clean : touch make.inc for dir in \ - CPV LAXlib FFTXlib UtilXlib upflib Modules PP PW EPW KS_Solvers \ + CPV LAXlib FFTXlib XClib UtilXlib upflib Modules PP PW EPW KS_Solvers \ NEB ACFDT COUPLE GWW XSpectra PWCOND dft-d3 \ atomic clib LR_Modules pwtools upflib \ dev-tools extlibs Environ TDDFPT PHonon HP GWW Doc GUI \ diff --git a/Modules/correlation_gga.f90 b/Modules/correlation_gga.f90 index 5d8508456..b78f125ee 100644 --- a/Modules/correlation_gga.f90 +++ b/Modules/correlation_gga.f90 @@ -1,7 +1,8 @@ ! MODULE corr_gga !corr_gga_gpu> ! - USE corr_lda, ONLY : pw, pw_spin !pw_spin_d,pw=>pw_d,corr_lda=>corr_lda_gpu> + USE ldaxc_interfaces, ONLY: pw + USE corr_lda, ONLY : pw_spin !pw_spin_d,pw=>pw_d,corr_lda=>corr_lda_gpu> ! CONTAINS ! diff --git a/Modules/correlation_lda_lsda.f90 b/Modules/correlation_lda_lsda.f90 index bb4cb438d..933996edd 100644 --- a/Modules/correlation_lda_lsda.f90 +++ b/Modules/correlation_lda_lsda.f90 @@ -8,475 +8,10 @@ MODULE corr_lda !corr_lda_gpu> +USE ldaxc_interfaces, ONLY: pz + CONTAINS -!------------------------------------------------------------------------- -SUBROUTINE pz( rs, iflag, ec, vc ) ! - !----------------------------------------------------------------------- - !! LDA parametrization from Monte Carlo DATA: - ! - !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981); - !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994). - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: iflag ! - !! see routine comments - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) - REAL(DP) :: lnrs, rs12, ox, dox - ! - DATA a / 0.0311d0, 0.031091d0 /, b / -0.048d0, -0.046644d0 /, & - c / 0.0020d0, 0.00419d0 /, d / -0.0116d0, -0.00983d0 / - ! - DATA gc / -0.1423d0, -0.103756d0 /, b1 / 1.0529d0, 0.56371d0 /, & - b2 / 0.3334d0, 0.27358d0 / - ! - IF ( rs < 1.0d0 ) THEN - ! - ! high density formula - lnrs = LOG(rs) - ! - ec = a(iflag)*lnrs + b(iflag) + c(iflag)*rs*lnrs + d(iflag)*rs - ! - vc = a(iflag)*lnrs + ( b(iflag) - a(iflag)/3.d0 ) + 2.d0/3.d0 * & - c(iflag)*rs*lnrs + ( 2.d0*d(iflag) - c(iflag) )/3.d0*rs - ELSE - ! - ! interpolation formula - rs12 = SQRT(rs) - ! - ox = 1.d0 + b1(iflag)*rs12 + b2(iflag)*rs - dox = 1.d0 + 7.d0/6.d0*b1(iflag)*rs12 + 4.d0/3.d0*b2(iflag)*rs - ! - ec = gc(iflag)/ox - vc = ec*dox/ox - ! - ENDIF - ! - RETURN - ! -END SUBROUTINE pz -! -! -!----------------------------------------------------------------------- -SUBROUTINE pzKZK( rs, ec, vc, vol ) ! - !----------------------------------------------------------------------- - !! LDA parametrization from Monte Carlo DATA: - ! - !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981) - !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - REAL(DP) :: vol ! - !! volume element - ! - ! ... local variables - ! - INTEGER :: iflag, kr - REAL(DP) :: ec0(2), vc0(2), ec0p - REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) - REAL(DP) :: lnrs, rs12, ox, dox, lnrsk, rsk - REAL(DP) :: a1, grs, g1, g2, g3, g4, dL, gh, gl, grsp - REAL(DP) :: f3, f2, f1, f0, pi - REAL(DP) :: D1, D2, D3, P1, P2, ry2h - ! - DATA a / 0.0311_DP, 0.031091_DP /, b / -0.048_DP, -0.046644_DP /, & - c / 0.0020_DP, 0.00419_DP /, d / -0.0116_DP, -0.00983_DP / - ! - DATA gc / -0.1423_DP, -0.103756_DP /, b1 / 1.0529_DP, 0.56371_DP /, & - b2 / 0.3334_DP, 0.27358_DP / - ! - DATA a1 / -2.2037_DP /, g1 / 0.1182_DP/, g2 / 1.1656_DP/, g3 / -5.2884_DP/, & - g4 / -1.1233_DP / - ! - DATA ry2h / 0.5_DP / - ! - iflag = 1 - pi = 4.d0 * ATAN(1.d0) - dL = vol**(1.d0/3.d0) - gh = 0.5d0 * dL / (2.d0 * pi)**(1.d0/3.d0) - gl = dL * (3.d0 / 2.d0 / pi)**(1.d0/3.d0) - ! - rsk = gh - ! - DO kr = 1, 2 - ! - lnrsk = LOG(rsk) - IF (rsk < 1.0d0) THEN - ! high density formula - ec0(kr) = a(iflag) *lnrsk + b(iflag) + c(iflag) * rsk * lnrsk + d(iflag) * rsk - vc0(kr) = a(iflag) * lnrsk + (b(iflag) - a(iflag) / 3.d0) + 2.d0 / & - 3.d0 * c (iflag) * rsk * lnrsk + (2.d0 * d (iflag) - c (iflag) ) & - / 3.d0 * rsk - ! - ELSE - ! interpolation formula - rs12 = SQRT(rsk) - ox = 1.d0 + b1 (iflag) * rs12 + b2 (iflag) * rsk - dox = 1.d0 + 7.d0 / 6.d0 * b1 (iflag) * rs12 + 4.d0 / 3.d0 * & - b2 (iflag) * rsk - ec0(kr) = gc (iflag) / ox - vc0(kr) = ec0(kr) * dox / ox - ! - ENDIF - ! - grs = g1 * rsk * lnrsk + g2 * rsk + g3 * rsk**1.5d0 + g4 * rsk**2.d0 - grsp = g1 * lnrsk + g1 + g2 + 1.5d0 * g3 * rsk**0.5d0 + 2.d0 * g4 * rsk - ! - ec0(kr) = ec0(kr) + (-a1 * rsk / dL**2.d0 + grs / dL**3.d0) * ry2h - vc0(kr) = vc0(kr) + (-2.d0 * a1 * rsk / dL**2.d0 / 3.d0 + & - grs / dL**3.d0 - grsp * rsk / 3.d0 / dL**3.d0) * ry2h - ! - rsk = rs - ! - ENDDO - ! - lnrs = LOG(rs) - ! - IF (rs <= gh) THEN - ec = ec0(2) - vc = vc0(2) - ELSE - IF ( rs <= gl ) THEN - ec0p = 3.d0 * (ec0(1) - vc0(1)) / gh - P1 = 3.d0 * ec0(1) - gh * ec0p - P2 = ec0p - D1 = gl - gh - D2 = gl**2.d0 - gh**2.d0 - D3 = gl**3.d0 - gh**3.d0 - f2 = 2.d0 * gl**2.d0 * P2 * D1 + D2 * P1 - f2 = f2 / (-(2.d0*gl*D1)**2.d0 + 4.d0*gl*D1*D2 - D2**2.d0 ) - f3 = - (P2 + 2.d0*D1*f2) / (3.d0 * D2) - f1 = - (P1 + D2 * f2) / (2.d0 * D1) - f0 = - gl * (gl * f2 + 2.d0 * f1) / 3.d0 - ! - ec = f3 * rs**3.d0 + f2 * rs**2.d0 + f1 * rs + f0 - vc = f2 * rs**2.d0 / 3.d0 + f1 * 2.d0 * rs / 3.d0 + f0 - ELSE - ec = 0.d0 - vc = 0.d0 - ENDIF - ENDIF - ! - ! - RETURN - ! -END SUBROUTINE pzKZK -! -! -!----------------------------------------------------------------------- -SUBROUTINE vwn( rs, ec, vc ) ! - !----------------------------------------------------------------------- - !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: a = 0.0310907d0, b = 3.72744d0, & - c = 12.9352d0, x0 = -0.10498d0 - REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt - ! - ! - q = SQRT( 4.d0*c - b*b ) - f1 = 2.d0*b/q - f2 = b*x0 / ( x0*x0 + b*x0 + c ) - f3 = 2.d0*( 2.d0 * x0 + b ) / q - ! - rs12 = SQRT(rs) - fx = rs + b*rs12 + c - qx = ATAN( q / (2.d0*rs12 + b) ) - ! - ec = a * ( LOG( rs/fx ) + f1*qx - f2*( LOG( (rs12 - x0)**2 / fx) & - + f3*qx) ) - ! - tx = 2.d0*rs12 + b - tt = tx*tx + q*q - ! - vc = ec - rs12*a / 6.d0*( 2.d0 / rs12 - tx/fx - 4.d0*b/tt - f2 * & - (2.d0 / (rs12 - x0) - tx/fx - 4.d0*(2.d0 * x0 + b)/tt) ) - ! - RETURN - ! -END SUBROUTINE vwn -! -! -!----------------------------------------------------------------------- -SUBROUTINE vwn1_rpa( rs, ec, vc ) ! - !----------------------------------------------------------------------- - !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: a = 0.0310907_DP, b = 13.0720_DP, & - c = 42.7198_DP, x0 = -0.409286_DP - REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt - ! - q = SQRT(4.d0*c - b*b) - f1 = 2.d0*b/q - f2 = b*x0 / (x0*x0 + b*x0 + c) - f3 = 2.d0 * (2.d0 * x0 + b) / q - ! - rs12 = SQRT(rs) - fx = rs + b * rs12 + c - qx = ATAN(q / (2.d0 * rs12 + b) ) - ! - ec = a*( LOG(rs/fx) + f1*qx - f2*(LOG((rs12 - x0)**2/fx) + f3*qx) ) - ! - tx = 2.d0*rs12 + b - tt = tx*tx + q*q - ! - vc = ec - rs12*a/6.d0*( 2.d0/rs12 - tx/fx - 4.d0*b/tt - f2 * & - (2.d0/(rs12 - x0) - tx/fx - 4.d0*(2.d0*x0 + b)/tt) ) - ! - ! - RETURN - ! -END SUBROUTINE vwn1_rpa -! -! -!----------------------------------------------------------------------- -SUBROUTINE lyp( rs, ec, vc ) ! - !----------------------------------------------------------------------- - !! C. Lee, W. Yang, and R.G. Parr, PRB 37, 785 (1988). - !! LDA part only. - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: a=0.04918d0, b=0.132d0*2.87123400018819108d0 - REAL(DP), PARAMETER :: pi43=1.61199195401647d0 - ! pi43=(4pi/3)^(1/3) - REAL(DP), PARAMETER :: c=0.2533d0*pi43, d=0.349d0*pi43 - REAL(DP) :: ecrs, ox - ! - ! - ecrs = b*EXP( -c*rs ) - ox = 1.d0 / (1.d0 + d*rs) - ! - ec = - a*ox*(1.d0 + ecrs) - vc = ec - rs/3.d0*a*ox*( d*ox + ecrs*(d*ox + c) ) - ! - RETURN - ! -END SUBROUTINE lyp -! -! -!----------------------------------------------------------------------- -SUBROUTINE pw( rs, iflag, ec, vc ) ! - !----------------------------------------------------------------------- - !! * iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) - !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - INTEGER, INTENT(IN) :: iflag ! - !! see routine comments - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: a=0.031091d0, b1=7.5957d0, b2=3.5876d0, c0=a, & - c1=0.046644d0, c2=0.00664d0, c3=0.01043d0, d0=0.4335d0, & - d1=1.4408d0 - REAL(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog - REAL(DP) :: a1(2), b3(2), b4(2) - DATA a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, & - b4 / 0.49294d0, 0.13354d0 / - ! - ! high- and low-density formulae implemented but not used in PW case - ! (reason: inconsistencies in PBE/PW91 functionals). - ! - IF ( rs < 1d0 .AND. iflag == 2 ) THEN - ! - ! high density formula - lnrs = LOG(rs) - ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs - vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * & - lnrs - (2.d0 * c3 + c2) / 3.d0 * rs - ! - ELSEIF ( rs > 100.d0 .AND. iflag == 2 ) THEN - ! - ! low density formula - ec = - d0 / rs + d1 / rs**1.5d0 - vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0 - ! - ELSE - ! - ! interpolation formula - rs12 = SQRT(rs) - rs32 = rs * rs12 - rs2 = rs**2 - om = 2.d0*a*( b1*rs12 + b2*rs + b3(iflag) * rs32 + b4(iflag)*rs2 ) - dom = 2.d0*a*( 0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3(iflag) * & - rs32 + 2.d0 * b4(iflag) * rs2 ) - olog = LOG( 1.d0 + 1.0d0 / om ) - ! - ec = - 2.d0 * a * (1.d0 + a1(iflag) * rs) * olog - vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1(iflag) * rs) & - * olog - 2.d0 / 3.d0 * a * (1.d0 + a1(iflag) * rs) * dom / & - (om * (om + 1.d0) ) - ! - ENDIF - ! - ! - RETURN - ! -END SUBROUTINE pw -! -! -!----------------------------------------------------------------------- -SUBROUTINE wignerc( rs, ec, vc ) ! - !----------------------------------------------------------------------- - !! Wigner correlation. - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP) :: rho13 !rho13=rho^(1/3) - REAL(DP), PARAMETER :: pi34 = 0.6203504908994d0 - ! pi34 = (3/4pi)^(1/3) - ! - ! - rho13 = pi34 / rs - vc = - rho13 * ( (0.943656d0 + 8.8963d0 * rho13) / (1.d0 + & - 12.57d0 * rho13)**2 ) - ec = - 0.738d0 * rho13 * ( 0.959d0 / (1.d0 + 12.57d0 * rho13) ) - ! - RETURN - ! -END SUBROUTINE wignerc -! -! -!----------------------------------------------------------------------- -SUBROUTINE hl( rs, ec, vc ) ! - !----------------------------------------------------------------------- - !! L. Hedin and B.I. Lundqvist, J. Phys. C 4, 2064 (1971). - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP) :: a, x - ! - ! - a = LOG(1.0d0 + 21.d0/rs) - x = rs / 21.0d0 - ec = a + (x**3 * a - x*x) + x/2.d0 - 1.0d0/3.0d0 - ec = - 0.0225d0 * ec - vc = - 0.0225d0 * a - ! - RETURN - ! -END SUBROUTINE hl -! -! -!----------------------------------------------------------------------- -SUBROUTINE gl( rs, ec, vc ) ! - !--------------------------------------------------------------------- - !! O. Gunnarsson and B. I. Lundqvist, PRB 13, 4274 (1976). - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ec - !! correlation energy - REAL(DP), INTENT(OUT) :: vc - !! correlation potential - ! - ! ... local variables - ! - REAL(DP) :: x - REAL(DP), PARAMETER :: c=0.0333d0, r=11.4d0 - ! c=0.0203, r=15.9 for the paramagnetic case - ! - x = rs / r - vc = - c * LOG(1.d0 + 1.d0 / x) - ec = - c * ( (1.d0 + x**3) * LOG(1.d0 + 1.d0 / x) - 1.0d0 / & - 3.0d0 + x * (0.5d0 - x) ) - ! - RETURN - ! -END SUBROUTINE gl ! ! ! ... LSDA diff --git a/Modules/dmxc_drivers.f90 b/Modules/dmxc_drivers.f90 index 87488d5b6..702e276f5 100644 --- a/Modules/dmxc_drivers.f90 +++ b/Modules/dmxc_drivers.f90 @@ -18,7 +18,8 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc ) ! USE kinds, ONLY: DP USE funct, ONLY: get_iexch, get_icorr, is_libxc - USE xc_lda_lsda, ONLY: xc_lda, xc_lsda + USE ldaxc_interfaces, ONLY: xc_lda + USE xc_lda_lsda, ONLY: xc_lsda !xc_lda #if defined(__LIBXC) #include "xc_version.h" USE xc_f03_lib_m @@ -182,8 +183,9 @@ SUBROUTINE dmxc_lda( length, rho_in, dmuxc ) !! Computes the derivative of the xc potential with respect to the !! local density. ! - USE xc_lda_lsda, ONLY: xc_lda - USE exch_lda, ONLY: slater + !USE xc_lda_lsda, ONLY: xc_lda + !USE exch_lda, ONLY: slater + USE ldaxc_interfaces, ONLY: xc_lda, slater, pz USE funct, ONLY: get_iexch, get_icorr USE kinds, ONLY: DP ! @@ -309,7 +311,8 @@ SUBROUTINE dmxc_lsda( length, rho_in, dmuxc ) USE kinds, ONLY: DP USE funct, ONLY: get_iexch, get_icorr USE xc_lda_lsda, ONLY: xc_lsda - USE exch_lda, ONLY: slater + !USE exch_lda, ONLY: slater + USE ldaxc_interfaces, ONLY: slater USE corr_lda, ONLY: pz, pz_polarized ! IMPLICIT NONE diff --git a/Modules/exchange_lda_lsda.f90 b/Modules/exchange_lda_lsda.f90 index 91524344a..c968af8a4 100644 --- a/Modules/exchange_lda_lsda.f90 +++ b/Modules/exchange_lda_lsda.f90 @@ -14,166 +14,6 @@ MODULE exch_lda !exch_lda_gpu> ! CONTAINS !----------------------------------------------------------------------- -SUBROUTINE slater( rs, ex, vx ) ! - !--------------------------------------------------------------------- - !! Slater exchange with alpha=2/3 - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - !! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ex - !! Exchange energy (per unit volume) - REAL(DP), INTENT(OUT) :: vx - !! Exchange potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: f = -0.687247939924714_DP, alpha = 2.0_DP/3.0_DP - ! f = -9/8*(3/2pi)^(2/3) - ex = f * alpha / rs - vx = 4._DP / 3._DP * f * alpha / rs - ! - RETURN - ! -END SUBROUTINE slater -! -! -!----------------------------------------------------------------------- -SUBROUTINE slater1( rs, ex, vx ) ! - !--------------------------------------------------------------------- - !! Slater exchange with alpha=1, corresponding to -1.374/r_s Ry. - !! Used to recover old results. - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ex - !! Exchange energy (per unit volume) - REAL(DP), INTENT(OUT) :: vx - !! Exchange potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: f = -0.687247939924714d0, alpha = 1.0d0 - ! - ex = f * alpha / rs - vx = 4.d0 / 3.d0 * f * alpha / rs - ! - RETURN - ! -END SUBROUTINE slater1 -! -! -!----------------------------------------------------------------------- -SUBROUTINE slater_rxc( rs, ex, vx ) ! - !--------------------------------------------------------------------- - !! Slater exchange with alpha=2/3 and Relativistic exchange. - ! - USE kinds, ONLY: DP - USE constants, ONLY: pi, c_au - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ex - !! Exchange energy (per unit volume) - REAL(DP), INTENT(OUT) :: vx - !! Exchange potential - ! - ! ... local variables - ! - REAL(DP), PARAMETER :: zero=0.d0, one=1.d0, pfive=0.5d0, & - opf=1.5d0 !, C014=0.014D0 - REAL(DP) :: trd, ftrd, tftm, a0, alp, z, fz, fzp, vxp, xp, & - beta, sb, alb, c014 - ! - trd = one/3.d0 - ftrd = 4.d0*trd - tftm = 2**ftrd-2.d0 - A0 = (4.d0/(9.d0*PI))**trd - C014 = 1.0_DP/a0/c_au - ! --- X-alpha PARAMETER: - alp = 2.d0 * trd - ! - z = zero - fz = zero - fzp= zero - ! - vxp = -3.d0*alp/( 2.d0*PI*A0*rs ) - xp = 3.d0*vxp/4.d0 - beta= C014 / rs - sb = SQRT(1.d0+beta*beta) - alb = LOG(beta+sb) - vxp = vxp * ( -pfive + opf * alb / (beta*sb) ) - xp = xp * ( one-opf*((beta*sb-alb)/beta**2)**2 ) - ! vxf = 2**trd*vxp - ! exf = 2**trd*xp - vx = vxp - ex = xp - ! - RETURN - ! -END SUBROUTINE slater_rxc -! -! -!----------------------------------------------------------------------- -SUBROUTINE slaterKZK( rs, ex, vx, vol ) ! - !--------------------------------------------------------------------- - !! Slater exchange with alpha=2/3, Kwee, Zhang and Krakauer KE - !! correction. - ! - USE kinds, ONLY: DP - ! - IMPLICIT NONE - ! - REAL(DP), INTENT(IN) :: rs - !! Wigner-Seitz radius - REAL(DP), INTENT(OUT) :: ex - !! Exchange energy (per unit volume) - REAL(DP), INTENT(OUT) :: vx - !! Exchange potential - REAL(DP) :: vol ! - !! Finite size volume element - ! - ! ... local variables - ! - REAL(DP) :: dL, ga, pi, a0 - REAL(DP), PARAMETER :: a1 = -2.2037d0, & - a2 = 0.4710d0, a3 = -0.015d0, ry2h = 0.5d0 - REAL(DP), PARAMETER :: f = -0.687247939924714d0, alpha = 2.0d0/3.0d0 - ! f = -9/8*(3/2pi)^(2/3) - ! - pi = 4.d0 * ATAN(1.d0) - a0 = f * alpha * 2.d0 - ! - dL = vol**(1.d0/3.d0) - ga = 0.5d0 * dL *(3.d0 /pi)**(1.d0/3.d0) - ! - IF ( rs < ga ) THEN - ex = a0 / rs + a1 * rs / dL**2.d0 + a2 * rs**2.d0 / dL**3.d0 - vx = (4.d0 * a0 / rs + 2.d0 * a1 * rs / dL**2.d0 + & - a2 * rs**2.d0 / dL**3.d0 ) / 3.d0 - ELSE - ex = a0 / ga + a1 * ga / dL**2.d0 + a2 * ga**2.d0 / dL**3.d0 ! solids - vx = ex - ! ex = a3 * dL**5.d0 / rs**6.d0 ! molecules - ! vx = 3.d0 * ex - ENDIF - ! - ex = ry2h * ex ! Ry to Hartree - vx = ry2h * vx - ! - RETURN - ! -END SUBROUTINE slaterKZK -! ! ! ... LSDA ! diff --git a/Modules/funct.f90 b/Modules/funct.f90 index 75758b5f4..472e57414 100644 --- a/Modules/funct.f90 +++ b/Modules/funct.f90 @@ -411,6 +411,11 @@ CONTAINS !----------------------------------------------------------------------- !! Translates a string containing the exchange-correlation name !! into internal indices iexch, icorr, igcx, igcc, inlc, imeta. + ! + USE ldaxc_interfaces, ONLY: get_ldaxc_indexes, get_ldaxc_param +#if defined(__LIBXC) + USE xc_f03_lib_m +#endif ! IMPLICIT NONE ! @@ -851,6 +856,25 @@ CONTAINS CALL errore( 'set_dft_from_name', ' conflicting values for inlc', 1 ) ENDIF ! + + + IF ( iexch==8 .OR. icorr==10 ) THEN ! 'sla+kzk' + ! + IF ( .NOT. finite_size_cell_volume_set) CALL errore('XC', & + 'finite size corrected exchange used w/o initialization',1) + ! + CALL get_ldaxc_indexes( iexch, icorr ) + CALL get_ldaxc_param( finite_size_cell_volume ) !...to set properly... + ! + ELSE + ! + CALL get_ldaxc_indexes( iexch, icorr ) + ! + ENDIF + + + + RETURN ! END SUBROUTINE set_dft_from_name diff --git a/Modules/make.depend b/Modules/make.depend index fef1a3f3c..40cb30954 100644 --- a/Modules/make.depend +++ b/Modules/make.depend @@ -60,8 +60,10 @@ constraints_module.o : kind.o control_flags.o : kind.o control_flags.o : parameters.o correlation_gga.o : correlation_lda_lsda.o +correlation_gga.o : ../XClib/ldaxc_interfaces.o correlation_gga.o : kind.o correlation_lda_lsda.o : kind.o +correlation_lda_lsda.o : ../XClib/ldaxc_interfaces.o cryst_to_car.o : kind.o deviatoric.o : io_global.o deviatoric.o : kind.o @@ -79,6 +81,7 @@ dmxc_drivers.o : correlation_lda_lsda.o dmxc_drivers.o : exchange_lda_lsda.o dmxc_drivers.o : funct.o dmxc_drivers.o : kind.o +dmxc_drivers.o : ../XClib/ldaxc_interfaces.o dmxc_drivers.o : xc_lda_lsda_drivers.o dylmr2.o : kind.o electrons_base.o : constants.o @@ -122,6 +125,7 @@ fox_init_module.o : mp_images.o funct.o : beef_interface.o funct.o : io_global.o funct.o : kind.o +funct.o : ../XClib/ldaxc_interfaces.o funct.o : xc_rVV10.o funct.o : xc_vdW_DF.o generate_function.o : ../FFTXlib/fft_types.o @@ -184,6 +188,7 @@ mbdlib.o : tsvdw.o mdiis.o : ../UtilXlib/mp.o mdiis.o : kind.o metagga.o : constants.o +metagga.o : ../XClib/ldaxc_interfaces.o metagga.o : correlation_gga.o metagga.o : correlation_lda_lsda.o metagga.o : exchange_lda_lsda.o @@ -416,6 +421,7 @@ xc_gga_drivers.o : funct.o xc_gga_drivers.o : kind.o xc_lda_lsda_drivers.o : correlation_lda_lsda.o xc_lda_lsda_drivers.o : exchange_lda_lsda.o +xc_lda_lsda_drivers.o : ../XClib/ldaxc_interfaces.o xc_lda_lsda_drivers.o : funct.o xc_lda_lsda_drivers.o : kind.o xc_mgga_drivers.o : funct.o @@ -433,6 +439,7 @@ xc_rVV10.o : mp_bands.o xc_rVV10.o : recvec.o xc_vdW_DF.o : ../FFTXlib/fft_interfaces.o xc_vdW_DF.o : ../UtilXlib/mp.o +xc_vdW_DF.o : ../XClib/ldaxc_interfaces.o xc_vdW_DF.o : cell_base.o xc_vdW_DF.o : constants.o xc_vdW_DF.o : control_flags.o diff --git a/Modules/metagga.f90 b/Modules/metagga.f90 index 50035c8ee..7fe27e5b7 100644 --- a/Modules/metagga.f90 +++ b/Modules/metagga.f90 @@ -20,8 +20,9 @@ ! MODULE metagga !metagga_gpu> ! -USE exch_lda, ONLY : slater !slater_d,exch_lda=>exch_lda_gpu> -USE corr_lda, ONLY : pw, pw_spin !pw_d,pw_spin=>pw_spin_d,corr_lda=>corr_lda_gpu> +USE ldaxc_interfaces, ONLY: slater, pw +!USE exch_lda, ONLY : slater !slater_d,exch_lda=>exch_lda_gpu> +USE corr_lda, ONLY : pw_spin !pw_d,pw_spin=>pw_spin_d,corr_lda=>corr_lda_gpu> USE corr_gga, ONLY : pbec, pbec_spin !pbec_d,pbec_spin=>pbec_spin_d,corr_gga=>corr_gga_gpu> ! CONTAINS diff --git a/Modules/xc_lda_lsda_drivers.f90 b/Modules/xc_lda_lsda_drivers.f90 index 86618045b..2672e8fac 100644 --- a/Modules/xc_lda_lsda_drivers.f90 +++ b/Modules/xc_lda_lsda_drivers.f90 @@ -12,6 +12,7 @@ USE kinds, ONLY: DP USE funct, ONLY: get_iexch, get_icorr, is_libxc, & exx_is_active, get_exx_fraction, & get_finite_size_cell_volume +USE ldaxc_interfaces, ONLY: xc_lda ! IMPLICIT NONE ! @@ -19,7 +20,7 @@ PRIVATE SAVE ! ! LDA and LSDA exchange-correlation drivers -PUBLIC :: xc, xc_lda, xc_lsda +PUBLIC :: xc, xc_lsda !xc_lda PUBLIC :: change_threshold_lda ! ! density threshold (set to default value) @@ -267,244 +268,6 @@ END SUBROUTINE xc !------- LDA-LSDA DRIVERS -------------------------------------------------- !---------------------------------------------------------------------------- ! -!---------------------------------------------------------------------------- -SUBROUTINE xc_lda( length, rho_in, ex_out, ec_out, vx_out, vc_out ) - !-------------------------------------------------------------------------- - !! LDA exchange and correlation functionals - Hartree a.u. - ! - !! * Exchange: - !! * Slater; - !! * relativistic Slater. - !! * Correlation: - !! * Ceperley-Alder (Perdew-Zunger parameters); - !! * Vosko-Wilk-Nusair; - !! * Lee-Yang-Parr; - !! * Perdew-Wang; - !! * Wigner; - !! * Hedin-Lundqvist; - !! * Ortiz-Ballone (Perdew-Zunger formula); - !! * Ortiz-Ballone (Perdew-Wang formula); - !! * Gunnarsson-Lundqvist. - ! - !! NOTE: - !! $$ E_x = \int E_x(\text{rho}) dr, E_x(\text{rho}) = - !! \text{rho}\epsilon_c(\text{rho})\ . $$ - !! Same for correlation. - ! - USE exch_lda - USE corr_lda - ! - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: length - !! length of the I/O arrays - REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in - !! Charge density - REAL(DP), INTENT(OUT), DIMENSION(length) :: ex_out - !! \(\epsilon_x(rho)\) ( NOT \(E_x(\text{rho})\) ) - REAL(DP), INTENT(OUT), DIMENSION(length) :: vx_out - !! \(dE_x(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_x(\text{rho})/d\text{rho}\) ) - REAL(DP), INTENT(OUT), DIMENSION(length) :: ec_out - !! \(\epsilon_c(rho)\) ( NOT \(E_c(\text{rho})\) ) - REAL(DP), INTENT(OUT), DIMENSION(length) :: vc_out - !! \(dE_c(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_c(\text{rho})/d\text{rho}\) ) - ! - ! ... local variables - ! - INTEGER :: ir, iexch, icorr - REAL(DP) :: rho, rs - REAL(DP) :: ex, ec, ec_ - REAL(DP) :: vx, vc, vc_ - REAL(DP) :: exx_fraction - REAL(DP) :: finite_size_cell_volume - LOGICAL :: exx_started, is_there_finite_size_corr - REAL(DP), PARAMETER :: third = 1.0_DP/3.0_DP, & - pi34 = 0.6203504908994_DP, e2 = 2.0_DP - ! pi34 = (3/4pi)^(1/3) - ! -#if defined(_OPENMP) - INTEGER :: ntids - INTEGER, EXTERNAL :: omp_get_num_threads - ! - ntids = omp_get_num_threads() -#endif - ! - iexch = get_iexch() - icorr = get_icorr() - exx_started = exx_is_active() - exx_fraction = get_exx_fraction() - IF (iexch==8 .OR. icorr==10) THEN - CALL get_finite_size_cell_volume( is_there_finite_size_corr, & - finite_size_cell_volume ) - ! - IF (.NOT. is_there_finite_size_corr) CALL errore( 'XC',& - 'finite size corrected exchange used w/o initialization', 1 ) - ENDIF - ! -!$omp parallel if(ntids==1) -!$omp do private( rho, rs, ex, ec, ec_, vx, vc, vc_ ) - DO ir = 1, length - ! - rho = ABS(rho_in(ir)) - ! - ! ... RHO THRESHOLD - ! - IF ( rho > rho_threshold ) THEN - rs = pi34 / rho**third - ELSE - ex_out(ir) = 0.0_DP ; ec_out(ir) = 0.0_DP - vx_out(ir) = 0.0_DP ; vc_out(ir) = 0.0_DP - CYCLE - ENDIF - ! - ! ... EXCHANGE - ! - SELECT CASE( iexch ) - CASE( 1 ) ! 'sla' - ! - CALL slater( rs, ex, vx ) - ! - CASE( 2 ) ! 'sl1' - ! - CALL slater1( rs, ex, vx ) - ! - CASE( 3 ) ! 'rxc' - ! - CALL slater_rxc( rs, ex, vx ) - ! - CASE( 4, 5 ) ! 'oep','hf' - ! - IF ( exx_started ) THEN - ex = 0.0_DP - vx = 0.0_DP - ELSE - CALL slater( rs, ex, vx ) - ENDIF - ! - CASE( 6, 7 ) ! 'pb0x' or 'DF-cx-0', or 'DF2-0', - ! ! 'B3LYP' - CALL slater( rs, ex, vx ) - IF ( exx_started ) THEN - ex = (1.0_DP - exx_fraction) * ex - vx = (1.0_DP - exx_fraction) * vx - ENDIF - ! - CASE( 8 ) ! 'sla+kzk' - ! - CALL slaterKZK( rs, ex, vx, finite_size_cell_volume ) - ! - CASE( 9 ) ! 'X3LYP' - ! - CALL slater( rs, ex, vx ) - IF ( exx_started ) THEN - ex = (1.0_DP - exx_fraction) * ex - vx = (1.0_DP - exx_fraction) * vx - ENDIF - ! - CASE DEFAULT - ! - ex = 0.0_DP - vx = 0.0_DP - ! - END SELECT - ! - ! - ! ... CORRELATION - ! - SELECT CASE( icorr ) - CASE( 1 ) - ! - CALL pz( rs, 1, ec, vc ) - ! - CASE( 2 ) - ! - CALL vwn( rs, ec, vc ) - ! - CASE( 3 ) - ! - CALL lyp( rs, ec, vc ) - ! - CASE( 4 ) - ! - CALL pw( rs, 1, ec, vc ) - ! - CASE( 5 ) - ! - CALL wignerc( rs, ec, vc ) - ! - CASE( 6 ) - ! - CALL hl( rs, ec, vc ) - ! - CASE( 7 ) - ! - CALL pz( rs, 2, ec, vc ) - ! - CASE( 8 ) - ! - CALL pw( rs, 2, ec, vc ) - ! - CASE( 9 ) - ! - CALL gl( rs, ec, vc ) - ! - CASE( 10 ) - ! - CALL pzKZK( rs, ec, vc, finite_size_cell_volume ) - ! - CASE( 11 ) - ! - CALL vwn1_rpa( rs, ec, vc ) - ! - CASE( 12 ) ! 'B3LYP' - ! - CALL vwn( rs, ec, vc ) - ec = 0.19_DP * ec - vc = 0.19_DP * vc - ! - CALL lyp( rs, ec_, vc_ ) - ec = ec + 0.81_DP * ec_ - vc = vc + 0.81_DP * vc_ - ! - CASE( 13 ) ! 'B3LYP-V1R' - ! - CALL vwn1_rpa ( rs, ec, vc ) - ec = 0.19_DP * ec - vc = 0.19_DP * vc - ! - CALL lyp( rs, ec_, vc_ ) - ec = ec + 0.81_DP * ec_ - vc = vc + 0.81_DP * vc_ - ! - CASE( 14 ) ! 'X3LYP' - ! - CALL vwn1_rpa( rs, ec, vc ) - ec = 0.129_DP * ec - vc = 0.129_DP * vc - ! - CALL lyp( rs, ec_, vc_ ) - ec = ec + 0.871_DP * ec_ - vc = vc + 0.871_DP * vc_ - ! - CASE DEFAULT - ! - ec = 0.0_DP - vc = 0.0_DP - ! - END SELECT - ! - ex_out(ir) = ex ; ec_out(ir) = ec - vx_out(ir) = vx ; vc_out(ir) = vc - ! - ENDDO -!$omp end do -!$omp end parallel - ! - ! - RETURN - ! -END SUBROUTINE xc_lda -! ! !----------------------------------------------------------------------------- SUBROUTINE xc_lsda( length, rho_in, zeta_in, ex_out, ec_out, vx_out, vc_out ) diff --git a/Modules/xc_vdW_DF.f90 b/Modules/xc_vdW_DF.f90 index 412444bd0..9c263caa1 100644 --- a/Modules/xc_vdW_DF.f90 +++ b/Modules/xc_vdW_DF.f90 @@ -723,6 +723,8 @@ CONTAINS SUBROUTINE get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas) + USE ldaxc_interfaces, ONLY: pw + IMPLICIT NONE REAL(DP), INTENT(IN) :: total_rho(:), grad_rho(:,:) ! Input variables needed. diff --git a/NEB/src/Makefile b/NEB/src/Makefile index 9b7a2d198..463fbb1be 100644 --- a/NEB/src/Makefile +++ b/NEB/src/Makefile @@ -39,7 +39,8 @@ stop_run_path.o PWOBJS= ../../PW/src/libpw.a QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \ - ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a + ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a \ + ../../XClib/libqeldaxc.a TLDEPS=pwlibs all : tldeps neb.x path_interpolation.x diff --git a/PHonon/FD/Makefile b/PHonon/FD/Makefile index c2a2efa6c..c5af3e882 100644 --- a/PHonon/FD/Makefile +++ b/PHonon/FD/Makefile @@ -12,7 +12,7 @@ FDOBJS = \ QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../upflib/libupf.a \ - ../../KS_Solvers/libks_solvers.a \ + ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a PWOBJS = ../../PW/src/libpw.a diff --git a/PHonon/Gamma/Makefile b/PHonon/Gamma/Makefile index 401569f5f..42e1e4ce4 100644 --- a/PHonon/Gamma/Makefile +++ b/PHonon/Gamma/Makefile @@ -35,7 +35,7 @@ solve_ph.o \ writedyn.o PWOBJS = ../../PW/src/libpw.a -QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \ +QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a ../../XClib/libqeldaxc.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a ../../LR_Modules/liblrmod.a TLDEPS= pwlibs diff --git a/PHonon/PH/Makefile b/PHonon/PH/Makefile index bbf812df1..1f0b33110 100644 --- a/PHonon/PH/Makefile +++ b/PHonon/PH/Makefile @@ -200,7 +200,7 @@ write_eigenvectors.o LRMODS = ../../LR_Modules/liblrmod.a PWOBJS = ../../PW/src/libpw.a -QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \ +QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a ../../XClib/libqeldaxc.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a TLDEPS= phlibs diff --git a/PP/src/Makefile b/PP/src/Makefile index 0e5a3fc79..21329d1ac 100644 --- a/PP/src/Makefile +++ b/PP/src/Makefile @@ -56,7 +56,7 @@ xc_vdW_scale_mod.o # added by Yang Jiao PWOBJS = ../../PW/src/libpw.a QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \ - ../../dft-d3/libdftd3qe.a + ../../dft-d3/libdftd3qe.a ../../XClib/libqeldaxc.a MODULES = $(PWOBJS) $(QEMODS) diff --git a/PP/src/make.depend b/PP/src/make.depend index 9ac6116b5..c1ba0165a 100644 --- a/PP/src/make.depend +++ b/PP/src/make.depend @@ -90,6 +90,25 @@ bands.o : ../../Modules/wavefunctions.o bands.o : ../../PW/src/pwcom.o bands.o : ../../UtilXlib/mp.o bands.o : ../../upflib/uspp.o +<<<<<<< HEAD +======= +<<<<<<< HEAD +xctest_qe_libxc.o : ../../Modules/correlation_lda_lsda.o +xctest_qe_libxc.o : ../../Modules/funct.o +xctest_qe_libxc.o : ../../Modules/xc_gga_drivers.o +xctest_qe_libxc.o : ../../Modules/xc_lda_lsda_drivers.o +xctest_qe_libxc.o : ../../Modules/xc_mgga_drivers.o +xctest_qe_libxc.o : ../../Modules/kind.o +======= +benchmark_libxc.o : ../../Modules/correlation_lda_lsda.o +benchmark_libxc.o : ../../Modules/funct.o +benchmark_libxc.o : ../../Modules/libxc.o +benchmark_libxc.o : ../../Modules/xc_gga_drivers.o +benchmark_libxc.o : ../../Modules/xc_lda_lsda_drivers.o +benchmark_libxc.o : ../../XClib/ldaxc_interfaces.o +benchmark_libxc.o : ../../Modules/xc_mgga_drivers.o +>>>>>>> XClib - lda - scratch +>>>>>>> XClib - lda - scratch chdens_bspline.o : ../../Modules/bspline.o chdens_bspline.o : ../../Modules/cell_base.o chdens_bspline.o : ../../Modules/fft_base.o diff --git a/PP/src/ppacf.f90 b/PP/src/ppacf.f90 index e52ddbb5f..7bb7ff640 100644 --- a/PP/src/ppacf.f90 +++ b/PP/src/ppacf.f90 @@ -56,7 +56,8 @@ PROGRAM do_ppacf USE funct, ONLY : get_iexch, get_icorr, get_igcx, get_igcc USE funct, ONLY : set_exx_fraction, set_auxiliary_flags, & enforce_input_dft - USE exch_lda, ONLY : slater, slater_spin + USE exch_lda, ONLY : slater_spin + USE ldaxc_interfaces, ONLY : slater USE xc_gga, ONLY : gcxc, gcx_spin, gcc_spin USE xc_lda_lsda, ONLY : xc USE wvfct, ONLY : npw, npwx diff --git a/PP/src/xc_vdW_scale_mod.f90 b/PP/src/xc_vdW_scale_mod.f90 index 63a9c9168..51db04f05 100644 --- a/PP/src/xc_vdW_scale_mod.f90 +++ b/PP/src/xc_vdW_scale_mod.f90 @@ -367,7 +367,8 @@ CONTAINS SUBROUTINE get_q0cc_on_grid (cc,lecnl_qx,total_rho, grad_rho, q0, thetas) USE vdW_DF, ONLY : spline_interpolation - USE corr_lda, ONLY : pw + !USE corr_lda, ONLY : pw + use ldaxc_interfaces, only : pw implicit none diff --git a/PW/src/Makefile b/PW/src/Makefile index 3e322c183..86b19af71 100644 --- a/PW/src/Makefile +++ b/PW/src/Makefile @@ -266,7 +266,7 @@ wannier_check.o \ wannier_clean.o \ wannier_occ.o -QEMODS=../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \ +QEMODS=../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a ../../XClib/libqeldaxc.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a TLDEPS= bindir libs mods libks_solvers dftd3 diff --git a/PWCOND/src/Makefile b/PWCOND/src/Makefile index 4b24ca6cc..8d741d80f 100644 --- a/PWCOND/src/Makefile +++ b/PWCOND/src/Makefile @@ -49,7 +49,7 @@ transmit.o PWOBJS = ../../PW/src/libpw.a QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/libks_solvers.a ../../upflib/libupf.a \ ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \ - ../../dft-d3/libdftd3qe.a + ../../dft-d3/libdftd3qe.a ../../XClib/libqeldaxc.a TLDEPS= pwlibs diff --git a/TDDFPT/src/Makefile b/TDDFPT/src/Makefile index 84284abc1..2840eb0e2 100644 --- a/TDDFPT/src/Makefile +++ b/TDDFPT/src/Makefile @@ -9,7 +9,7 @@ MODFLAGS= $(BASEMOD_FLAGS) \ $(MOD_FLAG)../../LR_Modules QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \ - ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a \ + ../../upflib/libupf.a ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a PWOBJS = ../../PW/src/libpw.a LRMODS = ../../LR_Modules/liblrmod.a diff --git a/XClib/Makefile b/XClib/Makefile new file mode 100644 index 000000000..dffe5cc4a --- /dev/null +++ b/XClib/Makefile @@ -0,0 +1,34 @@ +# Makefile for XClib + +include ../make.inc + +MODFLAGS= $(BASEMOD_FLAGS) + +LDAL = \ +kind_l.o \ +constants_l.o \ +dft_par_mod.o \ +ldaxc_interfaces.o \ +correlation_lda.o \ +get_ldaxc.o \ +exchange_lda.o \ +xc_lda.o + +all: libqeldaxc.a + +libqeldaxc.a: $(LDAL) + $(AR) $(ARFLAGS) $@ $? + $(RANLIB) $@ + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L + +# .PHONY forces execution of a rule irrespective of the presence of an +# updated file with the same name of the rule. In this way, the script +# that generates version.f90 always runs, updating the version if you +# execute "svn update". The update_version script takes care of not +# changing the file if the svn version did not change + +.PHONY: all clean + +include make.depend diff --git a/XClib/constants_l.f90 b/XClib/constants_l.f90 new file mode 100644 index 000000000..751f28932 --- /dev/null +++ b/XClib/constants_l.f90 @@ -0,0 +1,155 @@ +! +! Copyright (C) 2002-2006 Quantum ESPRESSO group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +!---------------------------------------------------------------------------- +MODULE constants_l + !---------------------------------------------------------------------------- + ! + USE kind_l, ONLY : DP + ! + ! ... The constants needed everywhere + ! + IMPLICIT NONE + ! + SAVE + ! + ! ... Mathematical constants + ! + REAL(DP), PARAMETER :: pi = 3.14159265358979323846_DP + REAL(DP), PARAMETER :: tpi = 2.0_DP * pi + REAL(DP), PARAMETER :: fpi = 4.0_DP * pi + REAL(DP), PARAMETER :: sqrtpi = 1.77245385090551602729_DP + REAL(DP), PARAMETER :: sqrtpm1= 1.0_DP / sqrtpi + REAL(DP), PARAMETER :: sqrt2 = 1.41421356237309504880_DP + ! + ! ... Physical constants, SI (NIST 2018) + ! http://physics.nist.gov/constants + ! + REAL(DP), PARAMETER :: H_PLANCK_SI = 6.62607015E-34_DP ! J s + REAL(DP), PARAMETER :: K_BOLTZMANN_SI = 1.380649E-23_DP ! J K^-1 + REAL(DP), PARAMETER :: ELECTRON_SI = 1.602176634E-19_DP ! C + REAL(DP), PARAMETER :: ELECTRONVOLT_SI = 1.602176634E-19_DP ! J + REAL(DP), PARAMETER :: ELECTRONMASS_SI = 9.1093837015E-31_DP ! Kg + REAL(DP), PARAMETER :: HARTREE_SI = 4.3597447222071E-18_DP ! J + REAL(DP), PARAMETER :: RYDBERG_SI = HARTREE_SI/2.0_DP ! J + REAL(DP), PARAMETER :: BOHR_RADIUS_SI = 0.529177210903E-10_DP ! m + REAL(DP), PARAMETER :: AMU_SI = 1.66053906660E-27_DP ! Kg + REAL(DP), PARAMETER :: C_SI = 2.99792458E+8_DP ! m sec^-1 + REAL(DP), PARAMETER :: MUNOUGHT_SI = fpi*1.0E-7_DP ! N A^-2 + REAL(DP), PARAMETER :: EPSNOUGHT_SI = 1.0_DP / (MUNOUGHT_SI * & + C_SI**2) ! F m^-1 + ! + ! ... Physical constants, atomic units: + ! ... AU for "Hartree" atomic units (e = m = hbar = 1) + ! ... RY for "Rydberg" atomic units (e^2=2, m=1/2, hbar=1) + ! + REAL(DP), PARAMETER :: K_BOLTZMANN_AU = K_BOLTZMANN_SI / HARTREE_SI + REAL(DP), PARAMETER :: K_BOLTZMANN_RY = K_BOLTZMANN_SI / RYDBERG_SI + ! + ! ... Unit conversion factors: energy and masses + ! + REAL(DP), PARAMETER :: AUTOEV = HARTREE_SI / ELECTRONVOLT_SI + REAL(DP), PARAMETER :: RYTOEV = AUTOEV / 2.0_DP + REAL(DP), PARAMETER :: AMU_AU = AMU_SI / ELECTRONMASS_SI + REAL(DP), PARAMETER :: AMU_RY = AMU_AU / 2.0_DP + ! + ! ... Unit conversion factors: atomic unit of time, in s and ps + ! + REAL(DP), PARAMETER :: AU_SEC = H_PLANCK_SI/tpi/HARTREE_SI + REAL(DP), PARAMETER :: AU_PS = AU_SEC * 1.0E+12_DP + ! + ! ... Unit conversion factors: pressure (1 Pa = 1 J/m^3, 1GPa = 10 Kbar ) + ! + REAL(DP), PARAMETER :: AU_GPA = HARTREE_SI / BOHR_RADIUS_SI ** 3 & + / 1.0E+9_DP + REAL(DP), PARAMETER :: RY_KBAR = 10.0_DP * AU_GPA / 2.0_DP + ! + ! ... Unit conversion factors: 1 debye = 10^-18 esu*cm + ! ... = 3.3356409519*10^-30 C*m + ! ... = 0.208194346 e*A + ! ... ( 1 esu = (0.1/c) Am, c=299792458 m/s) + ! + REAL(DP), PARAMETER :: DEBYE_SI = 3.3356409519_DP * 1.0E-30_DP ! C*m + REAL(DP), PARAMETER :: AU_DEBYE = ELECTRON_SI * BOHR_RADIUS_SI / & + DEBYE_SI + ! + REAL(DP), PARAMETER :: eV_to_kelvin = ELECTRONVOLT_SI / K_BOLTZMANN_SI + REAL(DP), PARAMETER :: ry_to_kelvin = RYDBERG_SI / K_BOLTZMANN_SI + ! + ! .. Unit conversion factors: Energy to wavelength + ! + REAL(DP), PARAMETER :: EVTONM = 1E+9_DP * H_PLANCK_SI * C_SI / & + &ELECTRONVOLT_SI + REAL(DP), PARAMETER :: RYTONM = 1E+9_DP * H_PLANCK_SI * C_SI / RYDBERG_SI + ! + ! Speed of light in atomic units + ! + REAL(DP), PARAMETER :: C_AU = C_SI / BOHR_RADIUS_SI * AU_SEC + ! + ! ... zero up to a given accuracy + ! + REAL(DP), PARAMETER :: eps4 = 1.0E-4_DP + REAL(DP), PARAMETER :: eps6 = 1.0E-6_DP + REAL(DP), PARAMETER :: eps8 = 1.0E-8_DP + REAL(DP), PARAMETER :: eps12 = 1.0E-12_DP + REAL(DP), PARAMETER :: eps14 = 1.0E-14_DP + REAL(DP), PARAMETER :: eps16 = 1.0E-16_DP + REAL(DP), PARAMETER :: eps24 = 1.0E-24_DP + REAL(DP), PARAMETER :: eps32 = 1.0E-32_DP + ! + REAL(DP), PARAMETER :: gsmall = 1.0E-12_DP + ! + REAL(DP), PARAMETER :: e2 = 2.0_DP ! the square of the electron charge + REAL(DP), PARAMETER :: degspin = 2.0_DP ! the number of spins per level + ! + !!!!!! COMPATIBIILITY + ! + REAL(DP), PARAMETER :: BOHR_RADIUS_CM = BOHR_RADIUS_SI * 100.0_DP + REAL(DP), PARAMETER :: BOHR_RADIUS_ANGS = BOHR_RADIUS_CM * 1.0E8_DP + REAL(DP), PARAMETER :: ANGSTROM_AU = 1.0_DP/BOHR_RADIUS_ANGS + REAL(DP), PARAMETER :: DIP_DEBYE = AU_DEBYE + REAL(DP), PARAMETER :: AU_TERAHERTZ = AU_PS + REAL(DP), PARAMETER :: AU_TO_OHMCMM1 = 46000.0_DP ! (ohm cm)^-1 + REAL(DP), PARAMETER :: RY_TO_THZ = 1.0_DP / AU_TERAHERTZ / FPI + REAL(DP), PARAMETER :: RY_TO_GHZ = RY_TO_THZ*1000.0_DP + REAL(DP), PARAMETER :: RY_TO_CMM1 = 1.E+10_DP * RY_TO_THZ / C_SI + ! + REAL(DP), PARAMETER :: AVOGADRO = 6.02214076D+23 + +END MODULE constants_l + +! perl script to create a program to list the available constants: +! extract with: grep '^!XX!' constants.f90 | sed 's,!XX!,,' > mkconstlist.pl +! then run: perl mkconstlist.pl constants.f90 > testme.f90 +! and compile and run: testme.f90 +!XX!#!/usr/bin/perl -w +!XX! +!XX!use strict; +!XX! +!XX!print <) { +!XX! if ( /REAL\s*\(DP\)\s*,\s*PARAMETER\s*::\s*([a-zA-Z_0-9]+)\s*=.*$/ ) { +!XX! print " WRITE (*,'(A18,G24.17)') '$1:',$1\n"; +!XX! } +!XX!} +!XX! +!XX!print <corr_lda_gpu> + +! CONTAINS + +!------------------------------------------------------------------------- +SUBROUTINE pz_l( rs, iflag, ec, vc ) ! + !----------------------------------------------------------------------- + !! LDA parametrization from Monte Carlo DATA: + ! + !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981); + !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994). + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + INTEGER, INTENT(IN) :: iflag ! + !! see routine comments + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) + REAL(DP) :: lnrs, rs12, ox, dox + ! + DATA a / 0.0311d0, 0.031091d0 /, b / -0.048d0, -0.046644d0 /, & + c / 0.0020d0, 0.00419d0 /, d / -0.0116d0, -0.00983d0 / + ! + DATA gc / -0.1423d0, -0.103756d0 /, b1 / 1.0529d0, 0.56371d0 /, & + b2 / 0.3334d0, 0.27358d0 / + ! + IF ( rs < 1.0d0 ) THEN + ! + ! high density formula + lnrs = LOG(rs) + ! + ec = a(iflag)*lnrs + b(iflag) + c(iflag)*rs*lnrs + d(iflag)*rs + ! + vc = a(iflag)*lnrs + ( b(iflag) - a(iflag)/3.d0 ) + 2.d0/3.d0 * & + c(iflag)*rs*lnrs + ( 2.d0*d(iflag) - c(iflag) )/3.d0*rs + ELSE + ! + ! interpolation formula + rs12 = SQRT(rs) + ! + ox = 1.d0 + b1(iflag)*rs12 + b2(iflag)*rs + dox = 1.d0 + 7.d0/6.d0*b1(iflag)*rs12 + 4.d0/3.d0*b2(iflag)*rs + ! + ec = gc(iflag)/ox + vc = ec*dox/ox + ! + ENDIF + ! + RETURN + ! +END SUBROUTINE pz_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE pzKZK_l( rs, ec, vc, vol ) ! + !----------------------------------------------------------------------- + !! LDA parametrization from Monte Carlo DATA: + ! + !! * iflag=1: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981) + !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + REAL(DP) :: vol ! + !! volume element + ! + ! ... local variables + ! + INTEGER :: iflag, kr + REAL(DP) :: ec0(2), vc0(2), ec0p + REAL(DP) :: a(2), b(2), c(2), d(2), gc(2), b1(2), b2(2) + REAL(DP) :: lnrs, rs12, ox, dox, lnrsk, rsk + REAL(DP) :: a1, grs, g1, g2, g3, g4, dL, gh, gl, grsp + REAL(DP) :: f3, f2, f1, f0, pi + REAL(DP) :: D1, D2, D3, P1, P2, ry2h + ! + DATA a / 0.0311_DP, 0.031091_DP /, b / -0.048_DP, -0.046644_DP /, & + c / 0.0020_DP, 0.00419_DP /, d / -0.0116_DP, -0.00983_DP / + ! + DATA gc / -0.1423_DP, -0.103756_DP /, b1 / 1.0529_DP, 0.56371_DP /, & + b2 / 0.3334_DP, 0.27358_DP / + ! + DATA a1 / -2.2037_DP /, g1 / 0.1182_DP/, g2 / 1.1656_DP/, g3 / -5.2884_DP/, & + g4 / -1.1233_DP / + ! + DATA ry2h / 0.5_DP / + ! + iflag = 1 + pi = 4.d0 * ATAN(1.d0) + dL = vol**(1.d0/3.d0) + gh = 0.5d0 * dL / (2.d0 * pi)**(1.d0/3.d0) + gl = dL * (3.d0 / 2.d0 / pi)**(1.d0/3.d0) + ! + rsk = gh + ! + DO kr = 1, 2 + ! + lnrsk = LOG(rsk) + IF (rsk < 1.0d0) THEN + ! high density formula + ec0(kr) = a(iflag) *lnrsk + b(iflag) + c(iflag) * rsk * lnrsk + d(iflag) * rsk + vc0(kr) = a(iflag) * lnrsk + (b(iflag) - a(iflag) / 3.d0) + 2.d0 / & + 3.d0 * c (iflag) * rsk * lnrsk + (2.d0 * d (iflag) - c (iflag) ) & + / 3.d0 * rsk + ! + ELSE + ! interpolation formula + rs12 = SQRT(rsk) + ox = 1.d0 + b1 (iflag) * rs12 + b2 (iflag) * rsk + dox = 1.d0 + 7.d0 / 6.d0 * b1 (iflag) * rs12 + 4.d0 / 3.d0 * & + b2 (iflag) * rsk + ec0(kr) = gc (iflag) / ox + vc0(kr) = ec0(kr) * dox / ox + ! + ENDIF + ! + grs = g1 * rsk * lnrsk + g2 * rsk + g3 * rsk**1.5d0 + g4 * rsk**2.d0 + grsp = g1 * lnrsk + g1 + g2 + 1.5d0 * g3 * rsk**0.5d0 + 2.d0 * g4 * rsk + ! + ec0(kr) = ec0(kr) + (-a1 * rsk / dL**2.d0 + grs / dL**3.d0) * ry2h + vc0(kr) = vc0(kr) + (-2.d0 * a1 * rsk / dL**2.d0 / 3.d0 + & + grs / dL**3.d0 - grsp * rsk / 3.d0 / dL**3.d0) * ry2h + ! + rsk = rs + ! + ENDDO + ! + lnrs = LOG(rs) + ! + IF (rs <= gh) THEN + ec = ec0(2) + vc = vc0(2) + ELSE + IF ( rs <= gl ) THEN + ec0p = 3.d0 * (ec0(1) - vc0(1)) / gh + P1 = 3.d0 * ec0(1) - gh * ec0p + P2 = ec0p + D1 = gl - gh + D2 = gl**2.d0 - gh**2.d0 + D3 = gl**3.d0 - gh**3.d0 + f2 = 2.d0 * gl**2.d0 * P2 * D1 + D2 * P1 + f2 = f2 / (-(2.d0*gl*D1)**2.d0 + 4.d0*gl*D1*D2 - D2**2.d0 ) + f3 = - (P2 + 2.d0*D1*f2) / (3.d0 * D2) + f1 = - (P1 + D2 * f2) / (2.d0 * D1) + f0 = - gl * (gl * f2 + 2.d0 * f1) / 3.d0 + ! + ec = f3 * rs**3.d0 + f2 * rs**2.d0 + f1 * rs + f0 + vc = f2 * rs**2.d0 / 3.d0 + f1 * 2.d0 * rs / 3.d0 + f0 + ELSE + ec = 0.d0 + vc = 0.d0 + ENDIF + ENDIF + ! + ! + RETURN + ! +END SUBROUTINE pzKZK_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE vwn_l( rs, ec, vc ) ! + !----------------------------------------------------------------------- + !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: a = 0.0310907d0, b = 3.72744d0, & + c = 12.9352d0, x0 = -0.10498d0 + REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt + ! + ! + q = SQRT( 4.d0*c - b*b ) + f1 = 2.d0*b/q + f2 = b*x0 / ( x0*x0 + b*x0 + c ) + f3 = 2.d0*( 2.d0 * x0 + b ) / q + ! + rs12 = SQRT(rs) + fx = rs + b*rs12 + c + qx = ATAN( q / (2.d0*rs12 + b) ) + ! + ec = a * ( LOG( rs/fx ) + f1*qx - f2*( LOG( (rs12 - x0)**2 / fx) & + + f3*qx) ) + ! + tx = 2.d0*rs12 + b + tt = tx*tx + q*q + ! + vc = ec - rs12*a / 6.d0*( 2.d0 / rs12 - tx/fx - 4.d0*b/tt - f2 * & + (2.d0 / (rs12 - x0) - tx/fx - 4.d0*(2.d0 * x0 + b)/tt) ) + ! + RETURN + ! +END SUBROUTINE vwn_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE vwn1_rpa_l( rs, ec, vc ) ! + !----------------------------------------------------------------------- + !! S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: a = 0.0310907_DP, b = 13.0720_DP, & + c = 42.7198_DP, x0 = -0.409286_DP + REAL(DP) :: q, f1, f2, f3, rs12, fx, qx, tx, tt + ! + q = SQRT(4.d0*c - b*b) + f1 = 2.d0*b/q + f2 = b*x0 / (x0*x0 + b*x0 + c) + f3 = 2.d0 * (2.d0 * x0 + b) / q + ! + rs12 = SQRT(rs) + fx = rs + b * rs12 + c + qx = ATAN(q / (2.d0 * rs12 + b) ) + ! + ec = a*( LOG(rs/fx) + f1*qx - f2*(LOG((rs12 - x0)**2/fx) + f3*qx) ) + ! + tx = 2.d0*rs12 + b + tt = tx*tx + q*q + ! + vc = ec - rs12*a/6.d0*( 2.d0/rs12 - tx/fx - 4.d0*b/tt - f2 * & + (2.d0/(rs12 - x0) - tx/fx - 4.d0*(2.d0*x0 + b)/tt) ) + ! + ! + RETURN + ! +END SUBROUTINE vwn1_rpa_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE lyp_l( rs, ec, vc ) ! + !----------------------------------------------------------------------- + !! C. Lee, W. Yang, and R.G. Parr, PRB 37, 785 (1988). + !! LDA part only. + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: a=0.04918d0, b=0.132d0*2.87123400018819108d0 + REAL(DP), PARAMETER :: pi43=1.61199195401647d0 + ! pi43=(4pi/3)^(1/3) + REAL(DP), PARAMETER :: c=0.2533d0*pi43, d=0.349d0*pi43 + REAL(DP) :: ecrs, ox + ! + ! + ecrs = b*EXP( -c*rs ) + ox = 1.d0 / (1.d0 + d*rs) + ! + ec = - a*ox*(1.d0 + ecrs) + vc = ec - rs/3.d0*a*ox*( d*ox + ecrs*(d*ox + c) ) + ! + RETURN + ! +END SUBROUTINE lyp_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE pw_l( rs, iflag, ec, vc ) ! + !----------------------------------------------------------------------- + !! * iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) + !! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + INTEGER, INTENT(IN) :: iflag ! + !! see routine comments + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: a=0.031091d0, b1=7.5957d0, b2=3.5876d0, c0=a, & + c1=0.046644d0, c2=0.00664d0, c3=0.01043d0, d0=0.4335d0, & + d1=1.4408d0 + REAL(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog + REAL(DP) :: a1(2), b3(2), b4(2) + DATA a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, & + b4 / 0.49294d0, 0.13354d0 / + ! + ! high- and low-density formulae implemented but not used in PW case + ! (reason: inconsistencies in PBE/PW91 functionals). + ! + IF ( rs < 1d0 .AND. iflag == 2 ) THEN + ! + ! high density formula + lnrs = LOG(rs) + ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs + vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * & + lnrs - (2.d0 * c3 + c2) / 3.d0 * rs + ! + ELSEIF ( rs > 100.d0 .AND. iflag == 2 ) THEN + ! + ! low density formula + ec = - d0 / rs + d1 / rs**1.5d0 + vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0 + ! + ELSE + ! + ! interpolation formula + rs12 = SQRT(rs) + rs32 = rs * rs12 + rs2 = rs**2 + om = 2.d0*a*( b1*rs12 + b2*rs + b3(iflag) * rs32 + b4(iflag)*rs2 ) + dom = 2.d0*a*( 0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3(iflag) * & + rs32 + 2.d0 * b4(iflag) * rs2 ) + olog = LOG( 1.d0 + 1.0d0 / om ) + ! + ec = - 2.d0 * a * (1.d0 + a1(iflag) * rs) * olog + vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1(iflag) * rs) & + * olog - 2.d0 / 3.d0 * a * (1.d0 + a1(iflag) * rs) * dom / & + (om * (om + 1.d0) ) + ! + ENDIF + ! + ! + RETURN + ! +END SUBROUTINE pw_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE wignerc_l( rs, ec, vc ) ! + !----------------------------------------------------------------------- + !! Wigner correlation. + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP) :: rho13 !rho13=rho^(1/3) + REAL(DP), PARAMETER :: pi34 = 0.6203504908994d0 + ! pi34 = (3/4pi)^(1/3) + ! + ! + rho13 = pi34 / rs + vc = - rho13 * ( (0.943656d0 + 8.8963d0 * rho13) / (1.d0 + & + 12.57d0 * rho13)**2 ) + ec = - 0.738d0 * rho13 * ( 0.959d0 / (1.d0 + 12.57d0 * rho13) ) + ! + RETURN + ! +END SUBROUTINE wignerc_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE hl_l( rs, ec, vc ) ! + !----------------------------------------------------------------------- + !! L. Hedin and B.I. Lundqvist, J. Phys. C 4, 2064 (1971). + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP) :: a, x + ! + ! + a = LOG(1.0d0 + 21.d0/rs) + x = rs / 21.0d0 + ec = a + (x**3 * a - x*x) + x/2.d0 - 1.0d0/3.0d0 + ec = - 0.0225d0 * ec + vc = - 0.0225d0 * a + ! + RETURN + ! +END SUBROUTINE hl_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE gl_l( rs, ec, vc ) ! + !--------------------------------------------------------------------- + !! O. Gunnarsson and B. I. Lundqvist, PRB 13, 4274 (1976). + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ec + !! correlation energy + REAL(DP), INTENT(OUT) :: vc + !! correlation potential + ! + ! ... local variables + ! + REAL(DP) :: x + REAL(DP), PARAMETER :: c=0.0333d0, r=11.4d0 + ! c=0.0203, r=15.9 for the paramagnetic case + ! + x = rs / r + vc = - c * LOG(1.d0 + 1.d0 / x) + ec = - c * ( (1.d0 + x**3) * LOG(1.d0 + 1.d0 / x) - 1.0d0 / & + 3.0d0 + x * (0.5d0 - x) ) + ! + RETURN + ! +END SUBROUTINE gl_l +! +! +! END MODULE diff --git a/XClib/dft_par_mod.f90 b/XClib/dft_par_mod.f90 new file mode 100644 index 000000000..71a796156 --- /dev/null +++ b/XClib/dft_par_mod.f90 @@ -0,0 +1,15 @@ +! +! ... LDAlib +! +MODULE dft_par_mod + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + INTEGER :: iexch, icorr + LOGICAL :: exx_started !, is_there_finite_size_corr + REAL(DP) :: finite_size_cell_volume, exx_fraction + REAL(DP) :: rho_threshold + ! +END MODULE diff --git a/XClib/exchange_lda.f90 b/XClib/exchange_lda.f90 new file mode 100644 index 000000000..31be95424 --- /dev/null +++ b/XClib/exchange_lda.f90 @@ -0,0 +1,176 @@ +! +! Copyright (C) 2001-2014 Quantum ESPRESSO group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +! +! MODULE exch_lda !exch_lda_gpu> +! +! ! Module containing LDA functionals. +! +! CONTAINS +!----------------------------------------------------------------------- +SUBROUTINE slater_l( rs, ex, vx ) ! + !--------------------------------------------------------------------- + !! Slater exchange with alpha=2/3 + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + !! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ex + !! Exchange energy (per unit volume) + REAL(DP), INTENT(OUT) :: vx + !! Exchange potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: f = -0.687247939924714_DP, alpha = 2.0_DP/3.0_DP + ! f = -9/8*(3/2pi)^(2/3) + ex = f * alpha / rs + vx = 4._DP / 3._DP * f * alpha / rs + ! + RETURN + ! +END SUBROUTINE slater_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE slater1_l( rs, ex, vx ) ! + !--------------------------------------------------------------------- + !! Slater exchange with alpha=1, corresponding to -1.374/r_s Ry. + !! Used to recover old results. + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ex + !! Exchange energy (per unit volume) + REAL(DP), INTENT(OUT) :: vx + !! Exchange potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: f = -0.687247939924714d0, alpha = 1.0d0 + ! + ex = f * alpha / rs + vx = 4.d0 / 3.d0 * f * alpha / rs + ! + RETURN + ! +END SUBROUTINE slater1_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE slater_rxc_l( rs, ex, vx ) ! + !--------------------------------------------------------------------- + !! Slater exchange with alpha=2/3 and Relativistic exchange. + ! + USE kind_l, ONLY: DP + USE constants_l, ONLY: pi, c_au + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ex + !! Exchange energy (per unit volume) + REAL(DP), INTENT(OUT) :: vx + !! Exchange potential + ! + ! ... local variables + ! + REAL(DP), PARAMETER :: zero=0.d0, one=1.d0, pfive=0.5d0, & + opf=1.5d0 !, C014=0.014D0 + REAL(DP) :: trd, ftrd, tftm, a0, alp, z, fz, fzp, vxp, xp, & + beta, sb, alb, c014 + ! + trd = one/3.d0 + ftrd = 4.d0*trd + tftm = 2**ftrd-2.d0 + A0 = (4.d0/(9.d0*PI))**trd + C014 = 1.0_DP/a0/c_au + ! --- X-alpha PARAMETER: + alp = 2.d0 * trd + ! + z = zero + fz = zero + fzp= zero + ! + vxp = -3.d0*alp/( 2.d0*PI*A0*rs ) + xp = 3.d0*vxp/4.d0 + beta= C014 / rs + sb = SQRT(1.d0+beta*beta) + alb = LOG(beta+sb) + vxp = vxp * ( -pfive + opf * alb / (beta*sb) ) + xp = xp * ( one-opf*((beta*sb-alb)/beta**2)**2 ) + ! vxf = 2**trd*vxp + ! exf = 2**trd*xp + vx = vxp + ex = xp + ! + RETURN + ! +END SUBROUTINE slater_rxc_l +! +! +!----------------------------------------------------------------------- +SUBROUTINE slaterKZK_l( rs, ex, vx, vol ) ! + !--------------------------------------------------------------------- + !! Slater exchange with alpha=2/3, Kwee, Zhang and Krakauer KE + !! correction. + ! + USE kind_l, ONLY: DP + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rs + !! Wigner-Seitz radius + REAL(DP), INTENT(OUT) :: ex + !! Exchange energy (per unit volume) + REAL(DP), INTENT(OUT) :: vx + !! Exchange potential + REAL(DP) :: vol ! + !! Finite size volume element + ! + ! ... local variables + ! + REAL(DP) :: dL, ga, pi, a0 + REAL(DP), PARAMETER :: a1 = -2.2037d0, & + a2 = 0.4710d0, a3 = -0.015d0, ry2h = 0.5d0 + REAL(DP), PARAMETER :: f = -0.687247939924714d0, alpha = 2.0d0/3.0d0 + ! f = -9/8*(3/2pi)^(2/3) + ! + pi = 4.d0 * ATAN(1.d0) + a0 = f * alpha * 2.d0 + ! + dL = vol**(1.d0/3.d0) + ga = 0.5d0 * dL *(3.d0 /pi)**(1.d0/3.d0) + ! + IF ( rs < ga ) THEN + ex = a0 / rs + a1 * rs / dL**2.d0 + a2 * rs**2.d0 / dL**3.d0 + vx = (4.d0 * a0 / rs + 2.d0 * a1 * rs / dL**2.d0 + & + a2 * rs**2.d0 / dL**3.d0 ) / 3.d0 + ELSE + ex = a0 / ga + a1 * ga / dL**2.d0 + a2 * ga**2.d0 / dL**3.d0 ! solids + vx = ex + ! ex = a3 * dL**5.d0 / rs**6.d0 ! molecules + ! vx = 3.d0 * ex + ENDIF + ! + ex = ry2h * ex ! Ry to Hartree + vx = ry2h * vx + ! + RETURN + ! +END SUBROUTINE slaterKZK_l +! +! +! END MODULE diff --git a/XClib/get_ldaxc.f90 b/XClib/get_ldaxc.f90 new file mode 100644 index 000000000..7c7298046 --- /dev/null +++ b/XClib/get_ldaxc.f90 @@ -0,0 +1,59 @@ +! +! ... LDAlib +! +SUBROUTINE get_ldaxclib( iexch_, icorr_ ) + ! + USE dft_par_mod + ! + IMPLICIT NONE + ! + INTEGER, INTENT(IN) :: iexch_, icorr_ + ! + iexch = iexch_ + icorr = icorr_ + ! + RETURN + ! +END SUBROUTINE + + +SUBROUTINE get_ldaxcparlib( finite_size_cell_volume_, exx_started_, exx_fraction_ ) + ! + USE kind_l, ONLY: DP + USE dft_par_mod + ! + IMPLICIT NONE + ! + REAL(DP), OPTIONAL, INTENT(IN) :: finite_size_cell_volume_ + LOGICAL, OPTIONAL, INTENT(IN) :: exx_started_ + REAL(DP), OPTIONAL, INTENT(IN) :: exx_fraction_ + ! + finite_size_cell_volume = -1.d0 + IF ( present(finite_size_cell_volume_) ) & + finite_size_cell_volume = finite_size_cell_volume_ + exx_started = .FALSE. + IF ( present(exx_started_) ) & + exx_started = exx_started_ + exx_fraction = 0.d0 + IF ( present(exx_fraction_) ) & + exx_fraction = exx_fraction_ + ! + RETURN + ! +END SUBROUTINE + + +SUBROUTINE get_threshold( rho_threshold_ ) + ! + USE kind_l, ONLY: DP + USE dft_par_mod + ! + IMPLICIT NONE + ! + REAL(DP), INTENT(IN) :: rho_threshold_ + ! + rho_threshold = rho_threshold_ + ! + RETURN + ! +END SUBROUTINE diff --git a/XClib/kind_l.f90 b/XClib/kind_l.f90 new file mode 100644 index 000000000..a0c966bbe --- /dev/null +++ b/XClib/kind_l.f90 @@ -0,0 +1,19 @@ +! +! ... LDAlib +! +MODULE kind_l +!------------------------------------------------------------------------------! +! #if defined(__MPI) +! #if defined(__MPI_MODULE) +! USE mpi +! #else +! INCLUDE 'mpif.h' + ! #endif + ! #endif + ! + IMPLICIT NONE + SAVE + INTEGER, PARAMETER :: DP = selected_real_kind(14,200) + PUBLIC :: DP + ! +END MODULE diff --git a/XClib/ldaxc_interfaces.f90 b/XClib/ldaxc_interfaces.f90 new file mode 100644 index 000000000..591a1ce3e --- /dev/null +++ b/XClib/ldaxc_interfaces.f90 @@ -0,0 +1,120 @@ +! +! ... LDAlib +! +!=---------------------------------------------------------------------------=! +MODULE ldaxc_interfaces + ! + IMPLICIT NONE + PRIVATE + ! + PUBLIC :: XC_LDA, SLATER, PZ, PW + PUBLIC :: GET_LDAXC_INDEXES, GET_LDAXC_PARAM, GET_LDA_THRESHOLD + ! + ! + INTERFACE GET_LDAXC_INDEXES + ! + SUBROUTINE get_ldaxclib( iexch_, icorr_ ) + ! + !USE dft_par_mod + USE kind_l, ONLY: DP + IMPLICIT NONE + INTEGER, INTENT(IN) :: iexch_, icorr_ + ! + END SUBROUTINE + ! + END INTERFACE + ! + INTERFACE GET_LDAXC_PARAM + ! + SUBROUTINE get_ldaxcparlib( finite_size_cell_volume_, exx_started_, exx_fraction_ ) + ! + !USE dft_par_mod + USE kind_l, ONLY: DP + IMPLICIT NONE + REAL(DP), OPTIONAL, INTENT(IN) :: finite_size_cell_volume_ + LOGICAL , OPTIONAL, INTENT(IN) :: exx_started_ + REAL(DP), OPTIONAL, INTENT(IN) :: exx_fraction_ + ! + END SUBROUTINE + ! + END INTERFACE + ! + INTERFACE GET_LDA_THRESHOLD + ! + SUBROUTINE get_threshold( rho_threshold_ ) + ! + !USE dft_par_mod + USE kind_l, ONLY: DP + IMPLICIT NONE + REAL(DP), INTENT(IN) :: rho_threshold_ + ! + END SUBROUTINE + ! + END INTERFACE + !-------------------------------------------------------- + ! + INTERFACE XC_LDA + ! + SUBROUTINE xc_lda_l( length, rho_in, ex_out, ec_out, vx_out, vc_out ) + ! + USE dft_par_mod + USE kind_l, ONLY: DP + IMPLICIT NONE + INTEGER, INTENT(IN) :: length + REAL(DP), INTENT(IN) :: rho_in(length) + REAL(DP), INTENT(OUT) :: ec_out(length), vc_out(length), & + ex_out(length), vx_out(length) + ! + END SUBROUTINE xc_lda_l + ! + END INTERFACE + ! + !-------------------------------------------------------- + ! + INTERFACE SLATER + ! + SUBROUTINE slater_l( rs, ex, vx ) + ! + USE kind_l, ONLY: DP + IMPLICIT NONE + REAL(DP), INTENT(IN) :: rs + REAL(DP), INTENT(OUT) :: ex + REAL(DP), INTENT(OUT) :: vx + ! + END SUBROUTINE slater_l + ! + END INTERFACE + ! + ! + INTERFACE PZ + ! + SUBROUTINE pz_l( rs, iflag, ec, vc ) + ! + USE kind_l, ONLY: DP + IMPLICIT NONE + REAL(DP), INTENT(IN) :: rs + REAL(DP), INTENT(OUT) :: ec, vc + INTEGER, INTENT(IN) :: iflag + ! + END SUBROUTINE pz_l + ! + END INTERFACE + ! + ! + INTERFACE PW + ! + SUBROUTINE pw_l( rs, iflag, ec, vc ) + ! + USE kind_l, ONLY: DP + IMPLICIT NONE + REAL(DP), INTENT(IN) :: rs + REAL(DP), INTENT(OUT) :: ec, vc + INTEGER, INTENT(IN) :: iflag + ! + END SUBROUTINE pw_l + ! + END INTERFACE + ! + ! +END MODULE ldaxc_interfaces +!=---------------------------------------------------------------------------=! diff --git a/XClib/make.depend b/XClib/make.depend new file mode 100644 index 000000000..d880251af --- /dev/null +++ b/XClib/make.depend @@ -0,0 +1,14 @@ +xc_lda.o : kind_l.o +xc_lda.o : dft_par_mod.o +xc_lda.o : correlation_lda.o +xc_lda.o : exchange_lda.o +ldaxc_interfaces.o : dft_par_mod.o +ldaxc_interfaces.o : kind_l.o +correlation_lda.o : kind_l.o +exchange_lda.o : kind_l.o +exchange_lda.o : constants_l.o +get_ldaxc.o : kind_l.o +get_ldaxc.o : dft_par_mod.o +constants_l.o : kind_l.o +dft_par_mod.o : kind_l.o +kind_l.o : diff --git a/XClib/xc_lda.f90 b/XClib/xc_lda.f90 new file mode 100644 index 000000000..05bc0150e --- /dev/null +++ b/XClib/xc_lda.f90 @@ -0,0 +1,242 @@ +!---------------------------------------------------------------------------- +SUBROUTINE xc_lda_l( length, rho_in, ex_out, ec_out, vx_out, vc_out ) + !-------------------------------------------------------------------------- + !! LDA exchange and correlation functionals - Hartree a.u. + ! + !! * Exchange: + !! * Slater; + !! * relativistic Slater. + !! * Correlation: + !! * Ceperley-Alder (Perdew-Zunger parameters); + !! * Vosko-Wilk-Nusair; + !! * Lee-Yang-Parr; + !! * Perdew-Wang; + !! * Wigner; + !! * Hedin-Lundqvist; + !! * Ortiz-Ballone (Perdew-Zunger formula); + !! * Ortiz-Ballone (Perdew-Wang formula); + !! * Gunnarsson-Lundqvist. + ! + !! NOTE: + !! $$ E_x = \int E_x(\text{rho}) dr, E_x(\text{rho}) = + !! \text{rho}\epsilon_c(\text{rho})\ . $$ + !! Same for correlation. + ! + USE kind_l, ONLY: DP + USE dft_par_mod + + ! USE exch_lda !-gpu + ! USE corr_lda !-gpu + ! + IMPLICIT NONE + ! + INTEGER, INTENT(IN) :: length + !! length of the I/O arrays + REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in + !! Charge density + REAL(DP), INTENT(OUT), DIMENSION(length) :: ex_out + !! \(\epsilon_x(rho)\) ( NOT \(E_x(\text{rho})\) ) + REAL(DP), INTENT(OUT), DIMENSION(length) :: vx_out + !! \(dE_x(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_x(\text{rho})/d\text{rho}\) ) + REAL(DP), INTENT(OUT), DIMENSION(length) :: ec_out + !! \(\epsilon_c(rho)\) ( NOT \(E_c(\text{rho})\) ) + REAL(DP), INTENT(OUT), DIMENSION(length) :: vc_out + !! \(dE_c(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_c(\text{rho})/d\text{rho}\) ) + ! + ! ... local variables + ! + INTEGER :: ir + REAL(DP) :: rho, rs + REAL(DP) :: ex, ec, ec_ + REAL(DP) :: vx, vc, vc_ + !REAL(DP) :: exx_fraction + !REAL(DP) :: finite_size_cell_volume + !LOGICAL :: exx_started, is_there_finite_size_corr + REAL(DP), PARAMETER :: third = 1.0_DP/3.0_DP, & + pi34 = 0.6203504908994_DP, e2 = 2.0_DP + ! pi34 = (3/4pi)^(1/3) + ! +#if defined(_OPENMP) + INTEGER :: ntids + INTEGER, EXTERNAL :: omp_get_num_threads + ! + ntids = omp_get_num_threads() +#endif + ! + !exx_started = exx_is_active() !....to adjust... + !exx_fraction = get_exx_fraction() + + rho_threshold = 1.E-10_DP !....out + + !IF (iexch==8 .OR. icorr==10) THEN + ! CALL get_finite_size_cell_volume( is_there_finite_size_corr, & + ! finite_size_cell_volume ) + ! + ! IF (.NOT. is_there_finite_size_corr) CALL errore( 'XC',& + ! 'finite size corrected exchange used w/o initialization', 1 ) + !ENDIF + ! +!$omp parallel if(ntids==1) +!$omp do private( rho, rs, ex, ec, ec_, vx, vc, vc_ ) + DO ir = 1, length + ! + rho = ABS(rho_in(ir)) + ! + ! ... RHO THRESHOLD + ! + IF ( rho > rho_threshold ) THEN + rs = pi34 / rho**third + ELSE + ex_out(ir) = 0.0_DP ; ec_out(ir) = 0.0_DP + vx_out(ir) = 0.0_DP ; vc_out(ir) = 0.0_DP + CYCLE + ENDIF + ! + ! ... EXCHANGE + ! + SELECT CASE( iexch ) + CASE( 1 ) ! 'sla' + ! + CALL slater_l( rs, ex, vx ) + ! + CASE( 2 ) ! 'sl1' + ! + CALL slater1_l( rs, ex, vx ) + ! + CASE( 3 ) ! 'rxc' + ! + CALL slater_rxc_l( rs, ex, vx ) + ! + CASE( 4, 5 ) ! 'oep','hf' + ! + IF ( exx_started ) THEN + ex = 0.0_DP + vx = 0.0_DP + ELSE + CALL slater_l( rs, ex, vx ) + ENDIF + ! + CASE( 6, 7 ) ! 'pb0x' or 'DF-cx-0', or 'DF2-0', + ! ! 'B3LYP' + CALL slater_l( rs, ex, vx ) + IF ( exx_started ) THEN + ex = (1.0_DP - exx_fraction) * ex + vx = (1.0_DP - exx_fraction) * vx + ENDIF + ! + CASE( 8 ) ! 'sla+kzk' + ! + CALL slaterKZK_l( rs, ex, vx, finite_size_cell_volume ) + ! + CASE( 9 ) ! 'X3LYP' + ! + CALL slater_l( rs, ex, vx ) + IF ( exx_started ) THEN + ex = (1.0_DP - exx_fraction) * ex + vx = (1.0_DP - exx_fraction) * vx + ENDIF + ! + CASE DEFAULT + ! + ex = 0.0_DP + vx = 0.0_DP + ! + END SELECT + ! + ! + ! ... CORRELATION + ! + SELECT CASE( icorr ) + CASE( 1 ) + ! + CALL pz_l( rs, 1, ec, vc ) + ! + CASE( 2 ) + ! + CALL vwn_l( rs, ec, vc ) + ! + CASE( 3 ) + ! + CALL lyp_l( rs, ec, vc ) + ! + CASE( 4 ) + ! + CALL pw_l( rs, 1, ec, vc ) + ! + CASE( 5 ) + ! + CALL wignerc_l( rs, ec, vc ) + ! + CASE( 6 ) + ! + CALL hl_l( rs, ec, vc ) + ! + CASE( 7 ) + ! + CALL pz_l( rs, 2, ec, vc ) + ! + CASE( 8 ) + ! + CALL pw_l( rs, 2, ec, vc ) + ! + CASE( 9 ) + ! + CALL gl_l( rs, ec, vc ) + ! + CASE( 10 ) + ! + CALL pzKZK_l( rs, ec, vc, finite_size_cell_volume ) + ! + CASE( 11 ) + ! + CALL vwn1_rpa_l( rs, ec, vc ) + ! + CASE( 12 ) ! 'B3LYP' + ! + CALL vwn_l( rs, ec, vc ) + ec = 0.19_DP * ec + vc = 0.19_DP * vc + ! + CALL lyp_l( rs, ec_, vc_ ) + ec = ec + 0.81_DP * ec_ + vc = vc + 0.81_DP * vc_ + ! + CASE( 13 ) ! 'B3LYP-V1R' + ! + CALL vwn1_rpa_l( rs, ec, vc ) + ec = 0.19_DP * ec + vc = 0.19_DP * vc + ! + CALL lyp_l( rs, ec_, vc_ ) + ec = ec + 0.81_DP * ec_ + vc = vc + 0.81_DP * vc_ + ! + CASE( 14 ) ! 'X3LYP' + ! + CALL vwn1_rpa_l( rs, ec, vc ) + ec = 0.129_DP * ec + vc = 0.129_DP * vc + ! + CALL lyp_l( rs, ec_, vc_ ) + ec = ec + 0.871_DP * ec_ + vc = vc + 0.871_DP * vc_ + ! + CASE DEFAULT + ! + ec = 0.0_DP + vc = 0.0_DP + ! + END SELECT + ! + ex_out(ir) = ex ; ec_out(ir) = ec + vx_out(ir) = vx ; vc_out(ir) = vc + ! + ENDDO +!$omp end do +!$omp end parallel + ! + ! + RETURN + ! +END SUBROUTINE xc_lda_l + diff --git a/XSpectra/src/Makefile b/XSpectra/src/Makefile index 320d07684..e10c9b121 100644 --- a/XSpectra/src/Makefile +++ b/XSpectra/src/Makefile @@ -39,7 +39,7 @@ GIPAWOBJS=./paw_gipaw.o \ MANIP_XS_OBJ=./gaunt_mod.o QEMODS = ../../Modules/libqemod.a ../../upflib/libupf.a \ - ../../FFTXlib/libqefft.a ../../KS_Solvers/libks_solvers.a \ + ../../FFTXlib/libqefft.a ../../KS_Solvers/libks_solvers.a ../../XClib/libqeldaxc.a \ ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a PWOBJS = ../../PW/src/libpw.a diff --git a/install/make.inc.in b/install/make.inc.in index fcf09fa14..2bef8180a 100644 --- a/install/make.inc.in +++ b/install/make.inc.in @@ -59,6 +59,7 @@ MOD_FLAG = @imod@ # Each Makefile can add directories to MODFLAGS and libraries to QEMODS BASEMOD_FLAGS= $(MOD_FLAG)$(TOPDIR)/upflib \ + $(MOD_FLAG)$(TOPDIR)/XClib \ $(MOD_FLAG)$(TOPDIR)/Modules \ $(MOD_FLAG)$(TOPDIR)/FFTXlib \ $(MOD_FLAG)$(TOPDIR)/LAXlib \ diff --git a/install/makedeps.sh b/install/makedeps.sh index 7e9b3686b..a22c2e015 100755 --- a/install/makedeps.sh +++ b/install/makedeps.sh @@ -18,7 +18,7 @@ then dirs=" LAXlib FFTXlib UtilXlib clib \ KS_Solvers/Davidson KS_Solvers/Davidson_RCI KS_Solvers/CG \ KS_Solvers/PPCG KS_Solvers/ParO KS_Solvers/DENSE \ - upflib Modules LR_Modules PW/src CPV/src PW/tools PP/src PWCOND/src \ + upflib XClib Modules LR_Modules PW/src CPV/src PW/tools PP/src PWCOND/src \ PHonon/Gamma PHonon/PH PHonon/FD HP/src atomic/src \ EPW/src XSpectra/src ACFDT/src NEB/src TDDFPT/src \ GWW/pw4gww GWW/gww GWW/head GWW/bse GWW/simple \ @@ -57,10 +57,10 @@ for dir in $dirs; do # default DEPENDS="$LEVEL1/include" # for convenience, used later - DEPEND1="$LEVEL1/include $LEVEL1/FFTXlib $LEVEL1/LAXlib $LEVEL1/UtilXlib \ + DEPEND1="$LEVEL1/include $LEVEL1/FFTXlib $LEVEL1/XClib $LEVEL1/LAXlib $LEVEL1/UtilXlib \ $LEVEL1/upflib" DEPEND3="$LEVEL2/include $LEVEL2/FFTXlib $LEVEL2/LAXlib $LEVEL2/UtilXlib" - DEPEND2="$DEPEND3 $LEVEL2/upflib $LEVEL2/Modules" + DEPEND2="$DEPEND3 $LEVEL2/upflib $LEVEL2/XClib $LEVEL2/Modules" case $DIR in Modules ) DEPENDS="$DEPEND1" ;; @@ -140,6 +140,10 @@ for dir in $dirs; do if test "$DIR" = "Modules" then sed '/@mbd@/d' make.depend > tmp; mv tmp make.depend + + if test "$DIR" = "XClib" + then + sed '/@elpa1@/d' make.depend > tmp; mv tmp make.depend fi if test "$DIR" = "PW/src" || test "$DIR" = "TDDFPT/src"