From 779b36c4d27eacae6fc8ccfb134457daec9e921a Mon Sep 17 00:00:00 2001 From: Frederick Stein Date: Thu, 13 May 2021 20:58:48 +0200 Subject: [PATCH] Remove SRLDA functional - LibXC has the same functional with higher derivatives available --- src/common/bibliography.F | 75 +- src/input_cp2k_xc.F | 44 +- src/xc/xc_derivatives.F | 11 - src/xc/xc_sr_lda.F | 917 ------------------ .../QS/regtest-hfx-wfn-fitting/CH4-rsLDA.inp | 17 +- tests/QS/regtest-hfx-wfn-fitting/TEST_FILES | 2 +- tests/QS/regtest-rs-dhft/CH3-rsLDAlrMP2.inp | 17 +- ...{H2O-srLDAlrMP2.inp => H2O-rsLDAlrMP2.inp} | 17 +- tests/QS/regtest-rs-dhft/TEST_FILES | 4 +- tests/TEST_DIRS | 4 +- 10 files changed, 49 insertions(+), 1059 deletions(-) delete mode 100644 src/xc/xc_sr_lda.F rename tests/QS/regtest-rs-dhft/{H2O-srLDAlrMP2.inp => H2O-rsLDAlrMP2.inp} (88%) diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 599ee648df..d9375e52e0 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -84,8 +84,8 @@ MODULE bibliography Brieuc2016, Barca2018, Scheiber2018, Huang2011, Heaton_Burgess2007, & Schuett2018, Holmberg2018, Togo2018, Staub2019, Grimme2013, Grimme2016, & Grimme2017, Kondov2007, Clabaut2020, & - Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, Paziani2006, & - Toulouse2004, Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, & + Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, & + Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, & Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Lass2018, cp2kqs2020, & Behler2007, Behler2011, Schran2020a, Schran2020b, & Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Shigeta2001, Heinecke2016, & @@ -4306,77 +4306,6 @@ CONTAINS "ER"), & DOI="10.1103/PhysRevB.61.8433") - CALL add_reference(key=Paziani2006, ISI_record=s2a( & - "AU Paziani, S", & - " Moroni, S", & - " Gori-Giorgi, P", & - " Bachelet, GB", & - "AF Paziani, S", & - " Moroni, S", & - " Gori-Giorgi, P", & - " Bachelet, GB", & - "TI Local-spin-density functional for multideterminant density functional", & - " theory", & - "SO PHYSICAL REVIEW B", & - "NR 62", & - "TC 64", & - "Z9 64", & - "PU AMER PHYSICAL SOC", & - "PI COLLEGE PK", & - "PA ONE PHYSICS ELLIPSE, COLLEGE PK, MD 20740-3844 USA", & - "SN 2469-9950", & - "EI 2469-9969", & - "J9 PHYS REV B", & - "JI Phys. Rev. B", & - "PD APR", & - "PY 2006", & - "VL 73", & - "IS 15", & - "AR 155111", & - "DI 10.1103/PhysRevB.73.155111", & - "PG 9", & - "SC Materials Science; Physics", & - "GA 037OA", & - "UT WOS:000237155100035", & - "ER"), & - DOI="10.1103/PhysRevB.73.155111") - - CALL add_reference(key=Toulouse2004, ISI_record=s2a( & - "AU Toulouse, J", & - " Savin, A", & - " Flad, HJ", & - "AF Toulouse, J", & - " Savin, A", & - " Flad, HJ", & - "TI Short-range exchange-correlation energy of a uniform electron gas with", & - " modified electron-electron interaction", & - "SO INTERNATIONAL JOURNAL OF QUANTUM CHEMISTRY", & - "CT 43rd International Symposium on Theory and Computations in Molecular and", & - " Materials Sciences, Biology, and Pharmacology", & - "RP Savin, A (reprint author), CNRS, Chim Theor Lab, 4 Pl Jussieu, F-75252 Paris, France.", & - "NR 28", & - "TC 90", & - "Z9 90", & - "PU WILEY-BLACKWELL", & - "PI HOBOKEN", & - "PA 111 RIVER ST, HOBOKEN 07030-5774, NJ USA", & - "SN 0020-7608", & - "EI 1097-461X", & - "J9 INT J QUANTUM CHEM", & - "JI Int. J. Quantum Chem.", & - "PD DEC 20", & - "PY 2004", & - "VL 100", & - "IS 6", & - "BP 1047", & - "EP 1056", & - "DI 10.1002/qua.20259", & - "PG 10", & - "GA 866QS", & - "UT WOS:000224788600025", & - "ER"), & - DOI="10.1002/qua.20259") - CALL add_reference(key=Limpanuparb2011, ISI_record=s2a( & "AU Limpanuparb, Taweetham", & " Gill, Peter M. W.", & diff --git a/src/input_cp2k_xc.F b/src/input_cp2k_xc.F index fc4b84c024..ffab3a684d 100644 --- a/src/input_cp2k_xc.F +++ b/src/input_cp2k_xc.F @@ -14,9 +14,8 @@ MODULE input_cp2k_xc USE bibliography, ONLY: & Becke1988, Becke1997, BeckeRoussel1989, Goedecker1996, Grimme2006, Grimme2010, Grimme2011, & - Heyd2004, Kruse2012, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Paziani2006, & - Perdew1981, Perdew1996, Perdew2008, Proynov2007, Tao2003, Toulouse2004, Tran2013, & - Vosko1980, Wellendorff2012, Zhang1998 + Heyd2004, Kruse2012, Lee1988, Lehtola2018, Marques2012, Ortiz1994, Perdew1981, Perdew1996, & + Perdew2008, Proynov2007, Tao2003, Tran2013, Vosko1980, Wellendorff2012, Zhang1998 USE cp_output_handling, ONLY: add_last_numeric,& cp_print_key_section_create,& high_print_level @@ -731,45 +730,6 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) - CALL section_create(subsection, __LOCATION__, name="SRLDA", & - description="Uses the short-range LDA functional", & - n_keywords=0, n_subsections=0, repeats=.FALSE., & - citations=(/Paziani2006, Toulouse2004/)) - CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & - description="activates the functional", & - lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="scale_x", & - description="scales the exchange part of the functional", & - default_r_val=1._dp) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="scale_c", & - description="scales the correlation part of the functional", & - default_r_val=1._dp) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="omega", & - description="provide the range-separation parameter of the functional", & - default_r_val=1._dp) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="PARAMETRIZATION", & - description="Which one parametrizations of the underlying PW92 functional should be used", & - usage="PARAMETRIZATION DMC", & - enum_c_vals=(/ & - "ORIGINAL", & - "DMC ", & - "VMC "/), & - enum_i_vals=(/c_pw92, c_pw92dmc, c_pw92vmc/), & - default_i_val=c_pw92) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - - CALL section_add_subsection(section, subsection) - CALL section_release(subsection) - END SUBROUTINE create_xc_fun_section ! ************************************************************************************************** diff --git a/src/xc/xc_derivatives.F b/src/xc/xc_derivatives.F index 98e207d3cc..72cf3eb381 100644 --- a/src/xc/xc_derivatives.F +++ b/src/xc/xc_derivatives.F @@ -67,9 +67,6 @@ MODULE xc_derivatives xc_rho_cflags_type USE xc_rho_set_types, ONLY: xc_rho_set_get,& xc_rho_set_type - USE xc_sr_lda, ONLY: sr_lda_eval,& - sr_lda_info,& - sr_lsd_eval USE xc_tfw, ONLY: tfw_lda_eval,& tfw_lda_info,& tfw_lsd_eval,& @@ -320,8 +317,6 @@ CONTAINS ELSE CALL xbr_pbe_lda_hole_tc_lr_lda_info(reference, shortform, needs, max_deriv) END IF - CASE ("SRLDA") - CALL sr_lda_info(reference, shortform, lsd, needs, max_deriv) CASE default ! If the functional has not been implemented internally, it's from LibXC IF (lsd) THEN @@ -549,12 +544,6 @@ CONTAINS CALL xbr_pbe_lda_hole_tc_lr_lda_eval(rho_set, deriv_set, deriv_order, & functional) END IF - CASE ("SRLDA") - IF (lsd) THEN - CALL sr_lsd_eval(rho_set, deriv_set, deriv_order, functional) - ELSE - CALL sr_lda_eval(rho_set, deriv_set, deriv_order, functional) - END IF CASE default ! If functional not natively supported, ask LibXC IF (lsd) THEN diff --git a/src/xc/xc_sr_lda.F b/src/xc/xc_sr_lda.F deleted file mode 100644 index f6bd5bdfd8..0000000000 --- a/src/xc/xc_sr_lda.F +++ /dev/null @@ -1,917 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright 2000-2021 CP2K developers group ! -! ! -! SPDX-License-Identifier: GPL-2.0-or-later ! -!--------------------------------------------------------------------------------------------------! - -! ************************************************************************************************** -!> \brief Calculates the short range correlation LDA energy (Improved version of Paola Gori-Giorgi's Code -!> \par History -!> 18-MAR-2002, TCH, working version -!> fawzi (04.2004) : adapted to the new xc interface -!> \see functionals_utilities -! ************************************************************************************************** -MODULE xc_sr_lda -#:include "xc_perdew_wang.fypp" - - USE bibliography, ONLY: Paziani2006, & - Toulouse2004, & - cite_reference - USE input_section_types, ONLY: section_vals_type, & - section_vals_val_get - USE kinds, ONLY: dp - USE mathconstants, ONLY: pi - USE xc_input_constants, ONLY: pw_dmc, & - pw_orig, & - pw_vmc - USE xc_derivative_set_types, ONLY: xc_derivative_set_type, & - xc_dset_get_derivative - USE xc_derivative_types, ONLY: xc_derivative_get, & - xc_derivative_type - USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type - USE xc_rho_set_types, ONLY: xc_rho_set_get, & - xc_rho_set_type - USE xc_functionals_utilities, ONLY: set_util -#include "../base/base_uses.f90" - - IMPLICIT NONE - PRIVATE - - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_sr_lda' - -@:global_var_pw92() - REAL(KIND=dp), PARAMETER, PRIVATE :: & - Acoul = 2._dp*(LOG(2._dp) - 1._dp)/pi**2, & - aQ2 = 5.84605_dp, & - cQ2 = 3.91744_dp, & - dQ2 = 3.44851_dp, & - bQ2 = dQ2 - 3._dp/(2._dp*pi*Acoul)*(4._dp/(9._dp*pi))**(1._dp/3._dp), & - f02 = 4._dp/(9._dp*(2._dp**(1._dp/3._dp) - 1._dp)), & - alpha = (4._dp/9._dp/pi)**(1._dp/3._dp), & - cf = (9._dp*pi/4._dp)**(1._dp/3._dp), & - p2p = 0.04_dp, & - p3p = 0.4319_dp, & - Cg0 = 0.0819306_dp, & - Fg0 = 0.752411_dp, & - Dg0 = -0.0127713_dp, & - Eg0 = 0.00185898_dp, & - Bg0 = 0.7317_dp - Fg0, & - adib = 0.784949_dp, & - q1a = -0.388_dp, & - q2a = 0.676_dp, & - q3a = 0.547_dp, & - t1a = -4.95_dp, & - t2a = 1._dp, & - t3a = 0.31_dp - - PUBLIC :: sr_lda_info, sr_lda_eval, sr_lsd_eval - -CONTAINS - -! ************************************************************************************************** -!> \brief Return some info on the functionals. -!> \param reference full reference -!> \param shortform short reference -!> \param lsd ... -!> \param needs ... -!> \param max_deriv ... -! ************************************************************************************************** - SUBROUTINE sr_lda_info(reference, shortform, lsd, needs, max_deriv) - CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform - LOGICAL, INTENT(IN), OPTIONAL :: lsd - TYPE(xc_rho_cflags_type), INTENT(INOUT), OPTIONAL :: needs - INTEGER, INTENT(OUT), OPTIONAL :: max_deriv - - CALL cite_reference(Toulouse2004) - CALL cite_reference(Paziani2006) - IF (PRESENT(reference)) THEN - reference = "J. Toulouse, A. Savin, and H.-J. Flad," & - //" Int. J. Quantum Chem. 100, 1074-1056 (2004)" - END IF - IF (PRESENT(shortform)) THEN - shortform = "J. Toulouse et al., IJQC 100, 1074-1056 (2004)" - END IF - IF (PRESENT(needs)) THEN - IF (lsd) THEN - needs%rho_spin = .TRUE. - ELSE - needs%rho = .TRUE. - END IF - END IF - IF (PRESENT(max_deriv)) max_deriv = 1 - - END SUBROUTINE sr_lda_info - -@:init_pw92() - -! ************************************************************************************************** -!> \brief Calculate the correlation energy and its derivatives -!> wrt to rho (the electron density) up to 3rd order. This -!> is the short-range LDA version of the Perdew-Wang correlation energy -!> If no order argument is given, then the routine calculates -!> just the energy. -!> \param rho_set ... -!> \param deriv_set ... -!> \param order order of derivatives to calculate -!> order must lie between -2 and 2. If it is negative then only -!> that order will be calculated, otherwise all derivatives up to -!> that order will be calculated. -!> \param sr_section ... -! ************************************************************************************************** - SUBROUTINE sr_lda_eval(rho_set, deriv_set, order, sr_section) - - TYPE(xc_rho_set_type), POINTER :: rho_set - TYPE(xc_derivative_set_type), POINTER :: deriv_set - INTEGER, INTENT(in) :: order - TYPE(section_vals_type), POINTER :: sr_section - - CHARACTER(len=*), PARAMETER :: routineN = 'sr_lda_eval' - - INTEGER :: npoints, handle, method - INTEGER, DIMENSION(:, :), POINTER :: bo - REAL(KIND=dp) :: omega, rho_cutoff, sc, sx - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dummy, e_0, e_rho, rho - TYPE(xc_derivative_type), POINTER :: deriv - - CALL timeset(routineN, handle) - - CALL section_vals_val_get(sr_section, 'SCALE_X', r_val=sx) - CALL section_vals_val_get(sr_section, 'SCALE_C', r_val=sc) - CALL section_vals_val_get(sr_section, 'OMEGA', r_val=omega) - CALL section_vals_val_get(sr_section, 'PARAMETRIZATION', i_val=method) - - NULLIFY (bo, rho, e_0, e_rho, dummy) - CPASSERT(ASSOCIATED(rho_set)) - CPASSERT(rho_set%ref_count > 0) - CPASSERT(ASSOCIATED(deriv_set)) - CPASSERT(deriv_set%ref_count > 0) - CALL xc_rho_set_get(rho_set, rho=rho, & - local_bounds=bo, rho_cutoff=rho_cutoff) - - CALL perdew_wang_init(method, rho_cutoff) - - npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) - - dummy => rho - - e_0 => dummy - e_rho => dummy - - IF (order >= 0) THEN - deriv => xc_dset_get_derivative(deriv_set, "", & - allocate_deriv=.TRUE.) - CALL xc_derivative_get(deriv, deriv_data=e_0) - END IF - IF (order >= 1 .OR. order == -1) THEN - deriv => xc_dset_get_derivative(deriv_set, "(rho)", & - allocate_deriv=.TRUE.) - CALL xc_derivative_get(deriv, deriv_data=e_rho) - END IF - IF (order > 1 .OR. order < -1) THEN - CPABORT("derivatives bigger than 1 not implemented") - END IF - - CALL sr_lda_calc(rho, omega, rho_cutoff, e_0, e_rho, & - npoints, order, sx, sc) - - CALL timestop(handle) - - END SUBROUTINE sr_lda_eval - -! ************************************************************************************************** -!> \brief ... -!> \param rho ... -!> \param omega ... -!> \param rho_cutoff ... -!> \param e_0 ... -!> \param e_rho ... -!> \param npoints ... -!> \param order ... -!> \param sx ... -!> \param sc ... -! ************************************************************************************************** - SUBROUTINE sr_lda_calc(rho, omega, rho_cutoff, e_0, e_rho, npoints, order, sx, sc) - !FM low level calc routine - REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho - REAL(KIND=dp), INTENT(IN) :: omega, rho_cutoff - REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho - INTEGER, INTENT(IN) :: npoints, order - REAL(KIND=dp), INTENT(IN) :: sx, sc - - CHARACTER(len=*), PARAMETER :: routineN = 'sr_lda_calc' - - INTEGER :: handle, k - REAL(KIND=dp) :: my_rho, rs - REAL(KIND=dp), DIMENSION(0:1) :: ed - - CALL timeset(routineN, handle) - - IF (sc /= 0.0_dp .OR. sx /= 0.0_dp) THEN -!$OMP PARALLEL DO PRIVATE (k, ed, my_rho, rs) DEFAULT(NONE)& -!$OMP SHARED(npoints,rho,rho_cutoff,omega,e_0,e_rho,order,sc,sx) - DO k = 1, npoints - - my_rho = rho(k) - IF (rho(k) > rho_cutoff) THEN - rs = (3.0_dp/4.0_dp/pi/my_rho)**(1.0_dp/3.0_dp) - - CALL ldasr(rs, omega, ed(0), ed(1), sx, sc) - - IF (order >= 0) THEN - e_0(k) = e_0(k) + ed(0)*my_rho - END IF - IF (order >= 1 .OR. order == -1) THEN - e_rho(k) = e_rho(k) + ed(1) - END IF - END IF - - END DO -!$OMP END PARALLEL DO - END IF - - CALL timestop(handle) - - END SUBROUTINE sr_lda_calc - -! ************************************************************************************************** -!> \brief Calculate the correlation energy and its derivatives -!> wrt to rho (the electron density) up to 3rd order. This -!> is the short-range LSD version of the Perdew-Wang correlation energy -!> If no order argument is given, then the routine calculates -!> just the energy. -!> \param rho_set ... -!> \param deriv_set ... -!> \param order order of derivatives to calculate -!> order must lie between -3 and 3. If it is negative then only -!> that order will be calculated, otherwise all derivatives up to -!> that order will be calculated. -!> \param sr_section ... -! ************************************************************************************************** - SUBROUTINE sr_lsd_eval(rho_set, deriv_set, order, sr_section) - - TYPE(xc_rho_set_type), POINTER :: rho_set - TYPE(xc_derivative_set_type), POINTER :: deriv_set - INTEGER, INTENT(IN), OPTIONAL :: order - TYPE(section_vals_type), POINTER :: sr_section - - CHARACTER(len=*), PARAMETER :: routineN = 'sr_lsd_eval' - - INTEGER :: npoints, handle, method - INTEGER, DIMENSION(:, :), POINTER :: bo - REAL(KIND=dp) :: omega, rho_cutoff, sc, sx - REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: a, b, dummy, e_0, ea, eb - TYPE(xc_derivative_type), POINTER :: deriv - - CALL timeset(routineN, handle) - - CALL section_vals_val_get(sr_section, 'SCALE_X', r_val=sx) - CALL section_vals_val_get(sr_section, 'SCALE_C', r_val=sc) - CALL section_vals_val_get(sr_section, 'OMEGA', r_val=omega) - CALL section_vals_val_get(sr_section, 'PARAMETRIZATION', i_val=method) - - NULLIFY (bo, a, b, e_0, ea, eb) - CPASSERT(ASSOCIATED(rho_set)) - CPASSERT(rho_set%ref_count > 0) - CPASSERT(ASSOCIATED(deriv_set)) - CPASSERT(deriv_set%ref_count > 0) - CALL xc_rho_set_get(rho_set, rhoa=a, rhob=b, & - local_bounds=bo, rho_cutoff=rho_cutoff) - - CALL perdew_wang_init(method, rho_cutoff) - - npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1) - - ! meaningful default for the arrays we don't need: let us make compiler - ! and debugger happy... - dummy => a - - e_0 => dummy - ea => dummy; eb => dummy - - IF (order >= 0) THEN - deriv => xc_dset_get_derivative(deriv_set, "", & - allocate_deriv=.TRUE.) - CALL xc_derivative_get(deriv, deriv_data=e_0) - END IF - IF (order >= 1 .OR. order == -1) THEN - deriv => xc_dset_get_derivative(deriv_set, "(rhoa)", & - allocate_deriv=.TRUE.) - CALL xc_derivative_get(deriv, deriv_data=ea) - deriv => xc_dset_get_derivative(deriv_set, "(rhob)", & - allocate_deriv=.TRUE.) - CALL xc_derivative_get(deriv, deriv_data=eb) - END IF - IF (order > 1 .OR. order < -1) THEN - CPABORT("derivatives bigger than 1 not implemented") - END IF - - CALL sr_lsd_calc(a, b, omega, rho_cutoff, e_0, ea, eb, npoints, order, sx, sc) - - CALL timestop(handle) - - END SUBROUTINE sr_lsd_eval - -! ************************************************************************************************** -!> \brief ... -!> \param rhoa ... -!> \param rhob ... -!> \param omega ... -!> \param rho_cutoff ... -!> \param e_0 ... -!> \param ea ... -!> \param eb ... -!> \param npoints ... -!> \param order ... -!> \param sx ... -!> \param sc ... -! ************************************************************************************************** - SUBROUTINE sr_lsd_calc(rhoa, rhob, omega, rho_cutoff, e_0, ea, eb, npoints, order, sx, sc) - !FM low-level computation routine - REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob - REAL(KIND=dp), INTENT(IN) :: omega, rho_cutoff - REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, ea, eb - INTEGER, INTENT(IN) :: npoints, order - REAL(KIND=dp), INTENT(IN) :: sx, sc - - CHARACTER(len=*), PARAMETER :: routineN = 'sr_lsd_calc' - - INTEGER :: handle, k - REAL(KIND=dp) :: my_rhoa, my_rhob, rho, rs, zeta - REAL(KIND=dp), DIMENSION(0:5) :: ed - - CALL timeset(routineN, handle) - - IF (sc /= 0.0_dp .OR. sx /= 0.0_dp) THEN -!$OMP PARALLEL DO PRIVATE (k, rho, ed, my_rhoa, my_rhob, rs, zeta) DEFAULT(NONE)& -!$OMP SHARED(npoints,rhoa,rhob,rho_cutoff,omega,order,e_0,ea,eb,sx,sc) - DO k = 1, npoints - - my_rhoa = rhoa(k) - my_rhob = rhob(k) - rho = my_rhoa + my_rhob - IF (rho > rho_cutoff) THEN - rs = (3.0_dp/4.0_dp/pi/rho)**(1.0_dp/3.0_dp) - zeta = (my_rhoa - my_rhob)/rho - - CALL lsdsr(rs, zeta, omega, ed(0), ed(1), ed(2), sx, sc) - IF (order >= 0) THEN - e_0(k) = e_0(k) + ed(0)*rho - END IF - IF (order >= 1 .OR. order == -1) THEN - ea(k) = ea(k) + ed(1) - eb(k) = eb(k) + ed(2) - END IF - END IF - - END DO -!$OMP END PARALLEL DO - - CALL timestop(handle) - - END IF - - END SUBROUTINE sr_lsd_calc - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param mu ... -!> \param excsr ... -!> \param vxcsr ... -!> \param sx ... -!> \param sc ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE ldasr(rs, mu, excsr, vxcsr, sx, sc) - REAL(KIND=dp), INTENT(IN) :: rs, mu - REAL(KIND=dp), INTENT(OUT) :: excsr, vxcsr - REAL(KIND=dp), INTENT(IN) :: sx, sc - - REAL(KIND=dp) :: ec, ecd, eclr, ex, exlr, vc, vclr, vx, & - vxlr - - IF (sx /= 0.0_dp) THEN - ex = -3._dp*cf/rs/4._dp/pi - vx = -(3._dp/2._dp/pi)**(2._dp/3._dp)/rs - - CALL exchangelr_lda(rs, mu, exlr, vxlr) - ELSE - ex = 0.0_dp - vx = 0.0_dp - - exlr = 0.0_dp - vxlr = 0.0_dp - END IF - - IF (sc /= 0.0_dp) THEN - CALL ecPW_lda(rs, ec, ecd) - vc = ec - rs/3._dp*ecd - - CALL ecorrlr_lda(rs, mu, eclr, vclr, ec, ecd) - ELSE - ec = 0.0_dp - vc = 0.0_dp - - eclr = 0.0_dp - vclr = 0.0_dp - END IF - - excsr = sx*ex + sc*ec - (sx*exlr + sc*eclr) - vxcsr = sx*vx + sc*vc - (sx*vxlr + sc*vclr) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param mu ... -!> \param eclr ... -!> \param vclr ... -!> \param ec ... -!> \param ecd ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE ecorrlr_lda(rs, mu, eclr, vclr, ec, ecd) - REAL(KIND=dp), INTENT(IN) :: rs, mu - REAL(KIND=dp), INTENT(OUT) :: eclr, vclr - REAL(KIND=dp), INTENT(IN) :: ec, ecd - - REAL(KIND=dp) :: a1, a1rs, a2, a2rs, a3, a3rs, a4, a4rs, a5, a5rs, b0, coe2, coe2rs, coe3, & - coe3rs, coe4, coe4rs, coe5, coe5rs, d2anti, d2antid, d3anti, d3antid, eclrrs, x, z - - b0 = adib*rs - z = 0._dp - - d2anti = (q1a + q2a*rs)*EXP(-q3a*rs)/rs - d3anti = (t1a + t2a*rs)*EXP(-t3a*rs)/rs**2 - - d2antid = -((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*EXP(-q3a*rs) - d3antid = -((rs*t2a*(1._dp + rs*t3a) + t1a*(2._dp + rs*t3a))/rs**3)*EXP(-rs*t3a) - - coe2 = -3._dp/8._dp/rs**3*(g0(rs) - 0.5_dp) - coe2rs = -3._dp/8._dp/rs**3*g0d(rs) + 9._dp/8._dp/rs**4*(g0(rs) - 0.5_dp) - - coe3 = -g0(rs)/SQRT(2._dp*pi)/rs**3 - coe3rs = -g0d(rs)/SQRT(2._dp*pi)/rs**3 + 3._dp*g0(rs)/SQRT(2._dp*pi)/rs**4 - - coe4 = -9._dp/64._dp/rs**3*(.5_dp*dpol(rs*2._dp**(1._dp/3._dp)) + d2anti - cf**2/5._dp/rs**2) - coe4rs = -3._dp/rs*coe4 - 9._dp/64._dp/rs**3*(((1._dp + z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp + z))** & - (1._dp/3._dp)) + ((1._dp - z)/2._dp)**(5._dp/3._dp)* & - dpold(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) + (1._dp - z**2)*d2antid & - + cf**2/5._dp*((1._dp + z)**(8._dp/3._dp) + (1._dp - z)**(8._dp/3._dp))/rs**3) - - coe5 = -9._dp/40._dp/SQRT(2._dp*pi)/rs**3*(((1._dp + z)/2._dp)**2*dpol(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) & - + ((1._dp - z)/2._dp)**2*dpol(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) & - + (1._dp - z**2)*d3anti) - coe5rs = -3._dp/rs*coe5 - 9._dp/(40._dp*SQRT(2._dp*pi)*rs**3)*( & - ((1._dp + z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) & - + ((1._dp - z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) + (1._dp - z**2)*d3antid) - - a1 = 4._dp*b0**6*coe3 + b0**8*coe5 - a1rs = 24._dp*adib*b0**5*coe3 + 4._dp*b0**6*coe3rs + 8._dp*adib*b0**7*coe5 + b0**8*coe5rs - - a2 = 4._dp*b0**6*coe2 + b0**8*coe4 + 6._dp*b0**4*ec - a2rs = 24._dp*adib*b0**5*coe2 + 4._dp*b0**6*coe2rs + 8._dp*adib*b0**7*coe4 + b0**8*coe4rs & - + 24._dp*adib*b0**3*ec + 6._dp*b0**4*ecd - - a3 = b0**8*coe3 - a3rs = 8._dp*adib*b0**7*coe3 + b0**8*coe3rs - - a4 = b0**6*(b0**2*coe2 + 4._dp*ec) - a4rs = 8._dp*adib*b0**7*coe2 + b0**8*coe2rs + 24._dp*adib*b0**5*ec + 4._dp*b0**6*ecd - - a5 = b0**8*ec - a5rs = 8._dp*adib*b0**7*ec + b0**8*ecd - - x = mu*SQRT(rs) - - eclr = (Qrpa(x) + mu**3*(a1 + mu*(a2 + mu*(a3 + mu*(a4 + a5*mu**2)))))/((1._dp + b0**2*mu**2)**4) - - eclrrs = -8._dp*adib/(1._dp + b0**2*mu**2)*b0*mu**2*eclr + & - 1._dp/((1._dp + b0**2*mu**2)**4)*(mu/(2._dp*SQRT(rs))*Qrpad(x) + & - mu**3*(a1rs + mu*(a2rs + mu*(a3rs + mu*(a4rs + a5rs*mu**2))))) - - vclr = eclr - rs/3._dp*eclrrs - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param z ... -!> \param mu ... -!> \param vxlrup ... -!> \param vxlrdown ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE exchangelr_lda(rs, mu, exlr, vxlr) - REAL(KIND=dp), INTENT(IN) :: rs, mu - REAL(KIND=dp), INTENT(OUT) :: exlr, vxlr - - REAL(KIND=dp) :: derrs, fx, fx1, y - - y = alpha/2._dp*mu*rs - fx = -((y*(-3._dp + 4._dp*y**2 + (2._dp - 4._dp*y**2)*EXP(-.25_dp/y**2)) + SQRT(pi)*ERF(.5_dp/y))/pi) - exlr = mu*fx - fx1 = (3._dp*(1._dp + (-4._dp + 4._dp*EXP(-.25_dp/y**2))*y**2))/pi - derrs = alpha/4._dp*mu**2*fx1 - vxlr = 2._dp/3._dp*rs*derrs - - vxlr = exlr - vxlr - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief PW92 energy functional -!> \param rs ... -!> \param ec ... -!> \param ecd ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE ecPW_lda(rs, ec, ecd) - REAL(KIND=dp), INTENT(IN) :: rs - REAL(KIND=dp), INTENT(OUT) :: ec, ecd - - REAL(KIND=dp) :: G(0:1) - - CALL calc_g(rs, 0, G, 1) - ec = G(0) - ecd = G(1) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param z ... -!> \param mu ... -!> \param excsr ... -!> \param vxcsrup ... -!> \param vxcsrdown ... -!> \param sx ... -!> \param sc ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE lsdsr(rs, z, mu, excsr, vxcsrup, vxcsrdown, sx, sc) - REAL(KIND=dp), INTENT(IN) :: rs, z, mu - REAL(KIND=dp), INTENT(OUT) :: excsr, vxcsrup, vxcsrdown - REAL(KIND=dp), INTENT(IN) :: sx, sc - - REAL(KIND=dp) :: ec, ecd, eclr, ecz, ex, exlr, vcdown, & - vclrdown, vclrup, vcup, vxdown, & - vxlrdown, vxlrup, vxup - - IF (sx /= 0.0_dp) THEN - ex = -3._dp*cf/rs/8._dp/pi*((1._dp + z)**(4._dp/3._dp) + & - (1._dp - z)**(4._dp/3._dp)) - - vxup = -(1._dp + z)**(1._dp/3._dp)*(3._dp/2._dp/pi)**(2._dp/3._dp)/rs - vxdown = -(1._dp - z)**(1._dp/3._dp)*(3._dp/2._dp/pi)**(2._dp/3._dp)/rs - - CALL exchangelr_lsd(rs, z, mu, exlr, vxlrup, vxlrdown) - ELSE - ex = 0.0_dp - vxup = 0.0_dp - vxdown = 0.0_dp - - exlr = 0.0_dp - vxlrup = 0.0_dp - vxlrdown = 0.0_dp - END IF - - IF (sc /= 0.0_dp) THEN - CALL ecPW_lsd(rs, z, ec, ecd, ecz) - vcup = ec - rs/3._dp*ecd - (z - 1._dp)*ecz - vcdown = ec - rs/3._dp*ecd - (z + 1._dp)*ecz - - CALL ecorrlr_lsd(rs, z, mu, eclr, vclrup, vclrdown, ec, ecd, ecz) - ELSE - ec = 0.0_dp - vcup = 0.0_dp - vcdown = 0.0_dp - - eclr = 0.0_dp - vclrup = 0.0_dp - vclrdown = 0.0_dp - END IF - - excsr = sx*ex + sc*ec - (sx*exlr + sc*eclr) - vxcsrup = sx*vxup + sc*vcup - (sx*vxlrup + sc*vclrup) - vxcsrdown = sx*vxdown + sc*vcdown - (sx*vxlrdown + sc*vclrdown) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param z ... -!> \param mu ... -!> \param eclr ... -!> \param ec ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE ecorrlr_lsd(rs, z, mu, eclr, vclrup, vclrdown, ec, ecd, ecz) - REAL(KIND=dp), INTENT(IN) :: rs, z, mu - REAL(KIND=dp), INTENT(OUT) :: eclr, vclrup, vclrdown - REAL(KIND=dp), INTENT(IN) :: ec, ecd, ecz - - REAL(KIND=dp) :: a1, a1rs, a1z, a2, a2rs, a2z, a3, a3rs, a3z, a4, a4rs, a4z, a5, a5rs, a5z, & - b0, coe2, coe2rs, coe2z, coe3, coe3rs, coe3z, coe4, coe4rs, coe4z, coe5, coe5rs, coe5z, & - d2anti, d2antid, d3anti, d3antid, eclrrs, eclrz, phi, x - - phi = ((1._dp + z)**(2._dp/3._dp) + (1._dp - z)**(2._dp/3._dp))/2._dp - - b0 = adib*rs - - d2anti = (q1a + q2a*rs)*EXP(-q3a*rs)/rs - d3anti = (t1a + t2a*rs)*EXP(-t3a*rs)/rs**2 - - d2antid = -((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*EXP(-q3a*rs) - d3antid = -((rs*t2a*(1._dp + rs*t3a) + t1a*(2._dp + rs*t3a))/rs**3)*EXP(-rs*t3a) - - coe2 = -3._dp/8._dp/rs**3*(1._dp - z**2)*(g0(rs) - 0.5_dp) - coe2rs = -3._dp/8._dp/rs**3*(1._dp - z**2)*g0d(rs) + 9._dp/8._dp/rs**4*(1._dp - z**2)*(g0(rs) - 0.5_dp) - coe2z = -3._dp/8._dp/rs**3*(-2._dp*z)*(g0(rs) - 0.5_dp) - - coe3 = -(1._dp - z**2)*g0(rs)/SQRT(2._dp*pi)/rs**3 - coe3rs = -(1._dp - z**2)*g0d(rs)/SQRT(2._dp*pi)/rs**3 + 3._dp*(1._dp - z**2)*g0(rs)/SQRT(2._dp*pi)/rs**4 - coe3z = 2._dp*z*g0(rs)/(SQRT(2._dp*pi)*rs**3) - - IF (ABS(z) >= 1._dp) THEN - - coe4 = -9._dp/64._dp/rs**3*(dpol(rs) - cf**2*2**(5._dp/3._dp)/5._dp/rs**2) - coe4rs = -3._dp/rs*coe4 - 9._dp/64._dp/rs**3*(dpold(rs) + 2._dp*cf**2*2**(5._dp/3._dp)/5._dp/rs**3) - coe4z = -9._dp/64._dp/rs**3*(dpol(rs) - rs/6._dp*dpold(rs) - 2._dp*d2anti & - - 4._dp/15._dp*cf**2*2._dp**(5._dp/3._dp)/rs**2)*z - - coe5 = -9._dp/40._dp/SQRT(2._dp*pi)/rs**3*dpol(rs) - coe5rs = -3._dp/rs*coe5 - 9._dp/40._dp/SQRT(2._dp*pi)/rs**3*dpold(rs) - coe5z = -9._dp/40._dp/SQRT(2._dp*pi)/rs**3*(dpol(rs) - rs/6._dp*dpold(rs) - 2._dp*d3anti)*z - - ELSE - - coe4 = -9._dp/64._dp/rs**3*(((1._dp + z)/2._dp)**2* & - dpol(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) + ((1._dp - z)/2._dp)**2 & - *dpol(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) + & - (1._dp - z**2)*d2anti - cf**2/10._dp*((1._dp + z)**(8._dp/3._dp) & - + (1._dp - z)**(8._dp/3._dp))/rs**2) - coe4rs = -3._dp/rs*coe4 - 9._dp/64._dp/rs**3*(((1._dp + z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp + z))** & - (1._dp/3._dp)) + ((1._dp - z)/2._dp)**(5._dp/3._dp)* & - dpold(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) + (1._dp - z**2)*d2antid & - + cf**2/5._dp*((1._dp + z)**(8._dp/3._dp) + (1._dp - z)**(8._dp/3._dp))/rs**3) - coe4z = -9._dp/64._dp/rs**3*(1._dp/2._dp*(1._dp + z)*dpol(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) & - - 1._dp/2._dp*(1._dp - z)*dpol(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) & - - rs/6._dp*((1._dp + z)/2._dp)**(2._dp/3._dp)*dpold(rs*(2/(1._dp + z))**(1._dp/3._dp)) & - + rs/6._dp*((1._dp - z)/2._dp)**(2._dp/3._dp)*dpold(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) & - - 2._dp*z*d2anti - 4._dp/15._dp*cf**2/rs**2*((1._dp + z)**(5._dp/3._dp) & - - (1._dp - z)**(5._dp/3._dp))) - - coe5 = -9._dp/40._dp/SQRT(2._dp*pi)/rs**3*(((1._dp + z)/2._dp)**2*dpol(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) & - + ((1._dp - z)/2._dp)**2*dpol(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) & - + (1._dp - z**2)*d3anti) - coe5rs = -3._dp/rs*coe5 - 9._dp/(40._dp*SQRT(2._dp*pi)*rs**3)*( & - ((1._dp + z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp + z))**(1._dp/3._dp)) & - + ((1._dp - z)/2._dp)**(5._dp/3._dp)*dpold(rs*(2._dp/(1._dp - z))**(1._dp/3._dp)) + (1._dp - z**2)*d3antid) - coe5z = -9._dp/40._dp/SQRT(2._dp*pi)/rs**3*(1._dp/2._dp*(1._dp + z)*dpol(rs*(2/(1._dp + z))**(1._dp/3._dp)) & - - 1._dp/2._dp*(1._dp - z)*dpol(rs*(2/(1._dp - z))**(1._dp/3._dp)) & - - rs/6._dp*((1._dp + z)/2._dp)**(2._dp/3._dp)*dpold(rs*(2/(1._dp + z)) & - **(1._dp/3._dp)) + rs/6._dp*((1._dp - z)/2._dp)**(2._dp/3._dp) & - *dpold(rs*(2/(1._dp - z))**(1._dp/3._dp)) - 2._dp*z*d3anti) - - END IF - - a1 = 4._dp*b0**6*coe3 + b0**8*coe5 - a1rs = 24._dp*adib*b0**5*coe3 + 4._dp*b0**6*coe3rs + 8._dp*adib*b0**7*coe5 + b0**8*coe5rs - a1z = 4._dp*b0**6*coe3z + b0**8*coe5z - - a2 = 4._dp*b0**6*coe2 + b0**8*coe4 + 6._dp*b0**4*ec - a2rs = 24._dp*adib*b0**5*coe2 + 4._dp*b0**6*coe2rs + 8._dp*adib*b0**7*coe4 + b0**8*coe4rs & - + 24._dp*adib*b0**3*ec + 6._dp*b0**4*ecd - a2z = 4._dp*b0**6*coe2z + b0**8*coe4z + 6._dp*b0**4*ecz - - a3 = b0**8*coe3 - a3rs = 8._dp*adib*b0**7*coe3 + b0**8*coe3rs - a3z = b0**8*coe3z - - a4 = b0**6*(b0**2*coe2 + 4._dp*ec) - a4rs = 8._dp*adib*b0**7*coe2 + b0**8*coe2rs + 24._dp*adib*b0**5*ec + 4._dp*b0**6*ecd - a4z = b0**6*(b0**2*coe2z + 4._dp*ecz) - - a5 = b0**8*ec - a5rs = 8._dp*adib*b0**7*ec + b0**8*ecd - a5z = b0**8*ecz - - x = mu*SQRT(rs)/phi - - eclr = (phi**3*Qrpa(x) + mu**3*(a1 + mu*(a2 + mu*(a3 + mu*(a4 + a5*mu**2)))))/((1._dp + b0**2*mu**2)**4) - - eclrrs = -8._dp*adib/(1._dp + b0**2*mu**2)*b0*mu**2*eclr + & - 1._dp/((1._dp + b0**2*mu**2)**4)*(phi**2*mu/(2._dp*SQRT(rs))*Qrpad(x) + & - mu**3*(a1rs + mu*(a2rs + mu*(a3rs + mu*(a4rs + a5rs*mu**2))))) - - IF (z >= 1._dp) THEN - vclrup = eclr - rs/3._dp*eclrrs - vclrdown = 0._dp - ELSE IF (z <= -1._dp) THEN - vclrup = 0._dp - vclrdown = eclr - rs/3._dp*eclrrs - ELSE - - eclrz = (phi**2*((1._dp + z)**(-1._dp/3._dp) - (1._dp - z)**(-1._dp/3._dp)) & - *Qrpa(x) - phi*Qrpad(x)*mu*SQRT(rs)*((1._dp + z)**(-1._dp/3._dp) & - - (1._dp - z)**(-1._dp/3._dp))/3._dp + & - mu**3*(a1z + mu*(a2z + mu*(a3z + mu*(a4z + a5z*mu**2)))))/((1._dp + b0**2*mu**2)**4) - - vclrup = eclr - rs/3._dp*eclrrs - (z - 1._dp)*eclrz - vclrdown = eclr - rs/3._dp*eclrrs - (z + 1._dp)*eclrz - END IF - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION g0(rs) - REAL(KIND=dp), INTENT(IN) :: rs - REAL(KIND=dp) :: g0 - - g0 = (1._dp - (0.7317_dp - Fg0)*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4)*EXP(-ABS(Fg0)*rs)/2._dp - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION g0d(rs) - REAL(KIND=dp), INTENT(IN) :: rs - REAL(KIND=dp) :: g0d - - g0d = (-Bg0 + 2.0_dp*Cg0*rs + 3.0_dp*Dg0*rs**2 + 4.0_dp*Eg0*rs**3)/2._dp*EXP(-Fg0*rs) & - - (Fg0*(1.0_dp - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/ & - 2._dp*EXP(-Fg0*rs) - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION dpol(rs) - REAL(KIND=dp), INTENT(IN) :: rs - REAL(KIND=dp) :: dpol - - dpol = 2._dp**(5._dp/3._dp)/5._dp*cf**2/rs**2*(1._dp + (p3p - 0.454555_dp)*rs) & - /(1._dp + p3p*rs + p2p*rs**2) - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION dpold(rs) - REAL(KIND=dp), INTENT(IN) :: rs - REAL(KIND=dp) :: dpold - - dpold = 2._dp**(5._dp/3._dp)/5._dp*cf**2* & - (-2._dp + (0.454555 - 4._dp*p3p)*rs + & - (-4._dp*p2p + (0.90911 - 2.*p3p)*p3p)*rs**2 & - + p2p*(1.363665 - 3._dp*p3p)*rs**3)/ & - (rs**3*(1._dp + p3p*rs + p2p*rs**2)**2) - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param x ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION Qrpa(x) - REAL(KIND=dp), INTENT(IN) :: x - REAL(KIND=dp) :: Qrpa - - Qrpa = Acoul*LOG((1._dp + aQ2*x + bQ2*x**2 + cQ2*x**3)/(1._dp + aQ2*x + dQ2*x**2)) - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param x ... -!> \return ... -! ************************************************************************************************** - ELEMENTAL FUNCTION Qrpad(x) - REAL(KIND=dp), INTENT(IN) :: x - REAL(KIND=dp) :: Qrpad - - Qrpad = Acoul*((x*(bQ2*(2._dp + aQ2*x) + cQ2*x*(3._dp + 2._dp*aQ2*x) + dQ2*(-2._dp - aQ2*x + cQ2*x**3)))/ & - ((1._dp + aQ2*x + dQ2*x**2)*(1._dp + aQ2*x + bQ2*x**2 + cQ2*x**3))) - - END FUNCTION - -! ************************************************************************************************** -!> \brief ... -!> \param rs ... -!> \param z ... -!> \param mu ... -!> \param vxlrup ... -!> \param vxlrdown ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE exchangelr_lsd(rs, z, mu, exlr, vxlrup, vxlrdown) - REAL(KIND=dp), INTENT(IN) :: rs, z, mu - REAL(KIND=dp), INTENT(OUT) :: exlr, vxlrup, vxlrdown - - REAL(KIND=dp) :: derrs, derz, fx, fx1, x, y - - IF (z >= 1._dp) THEN - x = rs*alpha*mu - y = .5_dp**(4._dp/3._dp)*x - fx = -((y*(-3._dp + 4._dp*y**2 + (2._dp - 4._dp*y**2)*EXP(-.25/y**2)) + SQRT(pi)*ERF(.5_dp/y))/pi) - exlr = mu*fx - vxlrup = mu*(x/(2._dp**(1._dp/3._dp)*pi) - x/(2._dp**(1._dp/3._dp)*pi)* & - EXP(-2._dp**(2._dp/3._dp)/x**2) - & - ERF(2._dp**(1._dp/3._dp)/x)/SQRT(pi)) - vxlrdown = 0._dp - ELSE IF (z <= -1._dp) THEN - x = rs*alpha*mu - y = .5_dp**(4._dp/3._dp)*x - fx = -((y*(-3._dp + 4._dp*y**2 + (2._dp - 4._dp*y**2)*EXP(-.25/y**2)) + SQRT(pi)*ERF(.5_dp/y))/pi) - exlr = mu*fx - vxlrdown = mu*(x/(2._dp**(1._dp/3._dp)*pi) - x/(2._dp**(1._dp/3._dp)*pi)* & - EXP(-2._dp**(2._dp/3._dp)/x**2) - & - ERF(2._dp**(1._dp/3._dp)/x)/SQRT(pi)) - vxlrup = 0._dp - ELSE - y = alpha/2._dp/(1.+z)**(1._dp/3._dp)*mu*rs - fx = -((y*(-3._dp + 4._dp*y**2 + (2._dp - 4._dp*y**2)*EXP(-.25_dp/y**2)) + & - SQRT(pi)*ERF(.5_dp/y))/pi) - exlr = (1._dp + z)*mu*fx/2._dp - fx1 = (3._dp*(1._dp + (-4._dp + 4._dp*EXP(-.25_dp/y**2))*y**2))/pi - derrs = alpha/4._dp*(1._dp + z)**(2._dp/3._dp)*mu**2*fx1 - derz = 1._dp/2._dp*mu*fx - 1._dp/6._dp*fx1*mu*y - vxlrup = rs/3._dp*derrs + (z - 1._dp)*derz - vxlrdown = rs/3._dp*derrs + (z + 1._dp)*derz - - y = alpha/2._dp/(1.-z)**(1._dp/3._dp)*mu*rs - fx = -((y*(-3._dp + 4._dp*y**2 + (2._dp - 4._dp*y**2)*EXP(-.25_dp/y**2)) + & - SQRT(pi)*ERF(.5_dp/y))/pi) - exlr = exlr + (1._dp - z)*mu*fx/2._dp - fx1 = (3._dp*(1._dp + (-4._dp + 4._dp*EXP(-.25_dp/y**2))*y**2))/pi - derrs = alpha/4._dp*(1._dp - z)**(2._dp/3._dp)*mu**2*fx1 - derz = -1._dp/2._dp*mu*fx + 1._dp/6._dp*fx1*mu*y - vxlrup = vxlrup + rs/3._dp*derrs + (z - 1._dp)*derz - vxlrdown = vxlrdown + rs/3._dp*derrs + (z + 1._dp)*derz - - vxlrup = exlr - vxlrup - vxlrdown = exlr - vxlrdown - ENDIF - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief PW92 energy functional -!> \param rs ... -!> \param z ... -!> \param ec ... -!> \param ecd ... -!> \param ecz ... -! ************************************************************************************************** - ELEMENTAL SUBROUTINE ecPW_lsd(rs, z, ec, ecd, ecz) - REAL(KIND=dp), INTENT(IN) :: rs, z - REAL(KIND=dp), INTENT(OUT) :: ec, ecd, ecz - - REAL(KIND=dp) :: alfac(0:1), ec0(0:1), ec1(0:1), ff - - IF (ABS(z) >= 1._dp) THEN - CALL calc_g(rs, 0, ec0, 0) - CALL calc_g(rs, 1, ec1, 1) - CALL calc_g(rs, -1, alfac, 0) - alfac = -alfac - - ec = ec1(0) - ecd = ec1(1) - ecz = SIGN(-4._dp/f02*alfac(0) + (ec1(0) - ec0(0))*(4._dp + 2._dp**(4._dp/3._dp)/3._dp/ & - (2._dp**(1._dp/3._dp) - 1._dp)), z) - ELSE - ff = ((1._dp + z)**(4._dp/3._dp) + (1._dp - z)**(4._dp/3._dp) - & - 2._dp)/(2._dp**(4._dp/3._dp) - 2._dp) - - CALL calc_g(rs, 0, ec0, 1) - CALL calc_g(rs, 1, ec1, 1) - CALL calc_g(rs, -1, alfac, 1) - alfac = -alfac - - ec = ec0(0) + alfac(0)*ff/f02*(1._dp - z**4) + (ec1(0) - ec0(0))*ff*z**4 - ecd = ec0(1) + alfac(1)*ff/f02*(1._dp - z**4) + (ec1(1) - ec0(1))*ff*z**4 - ecz = alfac(0)*(-4._dp*z**3)*ff/f02 + alfac(0)*(1._dp - z**4)/f02* & - 4._dp/3._dp*((1._dp + z)**(1._dp/3._dp) - (1._dp - z)**(1._dp/3._dp))/ & - (2._dp**(4._dp/3._dp) - 2._dp) + (ec1(0) - ec0(0))*(4._dp*z**3*ff + & - 4._dp/3._dp*((1._dp + z)**(1._dp/3._dp) - (1._dp - z)**(1._dp/3._dp))/ & - (2._dp**(4._dp/3._dp) - 2._dp)*z**4) - END IF - - END SUBROUTINE - -@:calc_g() - -END MODULE xc_sr_lda diff --git a/tests/QS/regtest-hfx-wfn-fitting/CH4-rsLDA.inp b/tests/QS/regtest-hfx-wfn-fitting/CH4-rsLDA.inp index 570caf41f8..04effe5206 100644 --- a/tests/QS/regtest-hfx-wfn-fitting/CH4-rsLDA.inp +++ b/tests/QS/regtest-hfx-wfn-fitting/CH4-rsLDA.inp @@ -1,3 +1,5 @@ +@SET MY_OMEGA 0.5 + &FORCE_EVAL METHOD Quickstep &DFT @@ -28,8 +30,17 @@ &END SCF &XC &XC_FUNCTIONAL - &SRLDA - OMEGA 0.5 + &LDA_X + &END LDA_X + &LDA_X_ERF + _OMEGA ${MY_OMEGA} + SCALE -1.0 + &END + &LDA_C_PMGB06 + _OMEGA ${MY_OMEGA} + SCALE -1.0 + &END + &LDA_C_PW &END &END XC_FUNCTIONAL &HF @@ -43,7 +54,7 @@ &END &INTERACTION_POTENTIAL POTENTIAL_TYPE LONGRANGE - OMEGA 0.5 + OMEGA ${MY_OMEGA} &END FRACTION 1.0 &END diff --git a/tests/QS/regtest-hfx-wfn-fitting/TEST_FILES b/tests/QS/regtest-hfx-wfn-fitting/TEST_FILES index 53dbd2a156..b59622424c 100644 --- a/tests/QS/regtest-hfx-wfn-fitting/TEST_FILES +++ b/tests/QS/regtest-hfx-wfn-fitting/TEST_FILES @@ -7,7 +7,7 @@ CH3-PBE0_TC.inp 1 1e-13 CH4-HSE06.inp 1 2e-13 -8.07752172778785 CH4-HSE06_2.inp 1 2e-13 -8.07752172778785 CH4-HSE06_TC_2.inp 1 2e-13 -8.07752172778785 -CH4-rsLDA.inp 1 6e-13 -8.07876568953425 +CH4-rsLDA.inp 1 6e-13 -8.48091146595489 CH4-PBE0.inp 1 2e-13 -8.07859057522753 CH4-PBE0_TC.inp 1 2e-13 -8.06493647354302 #EOF diff --git a/tests/QS/regtest-rs-dhft/CH3-rsLDAlrMP2.inp b/tests/QS/regtest-rs-dhft/CH3-rsLDAlrMP2.inp index 56aaaa8107..01edcbfaed 100644 --- a/tests/QS/regtest-rs-dhft/CH3-rsLDAlrMP2.inp +++ b/tests/QS/regtest-rs-dhft/CH3-rsLDAlrMP2.inp @@ -31,10 +31,19 @@ MAX_SCF 100 &END SCF &XC - &XC_FUNCTIONAL NONE - &SRLDA - OMEGA ${MY_OMEGA} - &END SRLDA + &XC_FUNCTIONAL + &LDA_X + &END LDA_X + &LDA_X_ERF + _OMEGA ${MY_OMEGA} + SCALE -1.0 + &END + &LDA_C_PMGB06 + SCALE -1.0 + _OMEGA ${MY_OMEGA} + &END + &LDA_C_PW + &END &END XC_FUNCTIONAL &HF FRACTION 1.0000000 diff --git a/tests/QS/regtest-rs-dhft/H2O-srLDAlrMP2.inp b/tests/QS/regtest-rs-dhft/H2O-rsLDAlrMP2.inp similarity index 88% rename from tests/QS/regtest-rs-dhft/H2O-srLDAlrMP2.inp rename to tests/QS/regtest-rs-dhft/H2O-rsLDAlrMP2.inp index 16c1e1c9b3..280ed19691 100644 --- a/tests/QS/regtest-rs-dhft/H2O-srLDAlrMP2.inp +++ b/tests/QS/regtest-rs-dhft/H2O-rsLDAlrMP2.inp @@ -32,10 +32,19 @@ ! ADDED_MOS 15000 15000 &END SCF &XC - &XC_FUNCTIONAL NONE - &SRLDA - OMEGA ${MY_OMEGA} - &END SRLDA + &XC_FUNCTIONAL + &LDA_X + &END LDA_X + &LDA_X_ERF + _OMEGA ${MY_OMEGA} + SCALE -1.0 + &END + &LDA_C_PMGB06 + SCALE -1.0 + _OMEGA ${MY_OMEGA} + &END + &LDA_C_PW + &END &END XC_FUNCTIONAL &HF FRACTION 1.0000000 diff --git a/tests/QS/regtest-rs-dhft/TEST_FILES b/tests/QS/regtest-rs-dhft/TEST_FILES index b6f0f4f7c3..9fb6661be4 100644 --- a/tests/QS/regtest-rs-dhft/TEST_FILES +++ b/tests/QS/regtest-rs-dhft/TEST_FILES @@ -1,3 +1,3 @@ -CH3-rsLDAlrMP2.inp 11 1e-8 -6.235769997143416 -H2O-srLDAlrMP2.inp 11 1e-8 -14.953293495365880 +CH3-rsLDAlrMP2.inp 11 1e-8 -7.724767496356234 +H2O-rsLDAlrMP2.inp 11 1e-8 -16.899255585036251 #EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index bb4a2f2915..98cf3a9f49 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -6,7 +6,7 @@ QS/regtest-grid QS/regtest-corr_dipm QS/regtest-admm-gapw libint -QS/regtest-rs-dhft libint +QS/regtest-rs-dhft libint libxc QS/regtest-sos-mp2-lr libint QS/regtest-rpa-lr libint QS/regtest-mp2-lr libint @@ -77,7 +77,7 @@ QS/regtest-hfx-block libint QS/regtest-ls-rtp QMMM/SE/regtest-force-mixing QS/regtest-xc -QS/regtest-hfx-wfn-fitting libint +QS/regtest-hfx-wfn-fitting libint libxc QS/regtest-tddfpt QS/regtest-tddfpt-stda libint QS/regtest-libxc libxc libint