From bd8365efd3eddc22f68c81605b30c6d4d18f45c0 Mon Sep 17 00:00:00 2001 From: Anna Hehn Date: Mon, 26 Jun 2023 11:41:26 +0200 Subject: [PATCH] GAPW triplet excitation energies --- src/qs_vxc_atom.F | 5 +- src/xc/xc.F | 10 ++- src/xc/xc_atom.F | 11 ++- .../QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp | 67 ++++++++++++++++ tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp | 76 +++++++++++++++++++ tests/QS/regtest-tddfpt/TEST_FILES | 2 + 6 files changed, 162 insertions(+), 9 deletions(-) create mode 100644 tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp create mode 100644 tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp diff --git a/src/qs_vxc_atom.F b/src/qs_vxc_atom.F index 42a8081f64..e390551408 100644 --- a/src/qs_vxc_atom.F +++ b/src/qs_vxc_atom.F @@ -546,6 +546,7 @@ CONTAINS CALL section_vals_val_get(input, "DFT%TDDFPT%RES_ETYPE", & i_val=res_etype) END IF + xc_fun_section => section_vals_get_subs_vals(xc_section, & "XC_FUNCTIONAL") IF (lsd) THEN @@ -710,11 +711,11 @@ CONTAINS CALL xc_2nd_deriv_of_r(xc_section=xc_section, & rho_set=rho_set_h, rho1_set=rho1_set_h, & deriv_set=deriv_set, & - w=weight, vxc=vxc_h, vxg=vxg_h) + w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=do_triplet) CALL xc_2nd_deriv_of_r(xc_section=xc_section, & rho_set=rho_set_s, rho1_set=rho1_set_s, & deriv_set=deriv_set, & - w=weight, vxc=vxc_s, vxg=vxg_s) + w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=do_triplet) CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss) IF (gradient_functional) THEN diff --git a/src/xc/xc.F b/src/xc/xc.F index 7f9ca0e2bb..842e3d8ddf 100644 --- a/src/xc/xc.F +++ b/src/xc/xc.F @@ -2623,7 +2623,7 @@ CONTAINS INTEGER, DIMENSION(2, 3) :: bo LOGICAL :: gradient_f, lsd, my_compute_virial, & my_gapw, tau_f, laplace_f, rho_f - REAL(KIND=dp) :: fac, gradient_cut, tmp + REAL(KIND=dp) :: fac, gradient_cut, tmp, factor2 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dr1dr, dra1dra, drb1drb REAL(kind=dp), DIMENSION(:, :, :), POINTER :: deriv_data, e_drhoa, e_drhob, & e_drho, norm_drho, norm_drhoa, & @@ -2661,7 +2661,9 @@ CONTAINS nspins = SIZE(v_xc) lsd = ASSOCIATED(rho_set%rhoa) fac = 0.0_dp + factor2 = 1.0_dp IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac + IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac bo = rho_set%local_bounds @@ -2938,12 +2940,12 @@ CONTAINS DO idir = 1, 3 !$OMP PARALLEL DO PRIVATE(ia,ir) DEFAULT(NONE) & -!$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir) COLLAPSE(2) +!$OMP SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) COLLAPSE(2) DO ia = bo(1, 1), bo(2, 1) DO ir = bo(1, 2), bo(2, 2) vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%cr3d(ia, ir, 1) IF (ASSOCIATED(e_drho)) THEN - vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1) + vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1) END IF END DO END DO @@ -3401,7 +3403,7 @@ CONTAINS SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac) REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & INTENT(OUT) :: dr1dr - TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b + TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN) :: drhoa, drhob, drho1a, drho1b REAL(KIND=dp), INTENT(IN) :: fac CHARACTER(len=*), PARAMETER :: routineN = 'prepare_dr1dr_ab' diff --git a/src/xc/xc_atom.F b/src/xc/xc_atom.F index c3fbef2dc8..a575f4d391 100644 --- a/src/xc/xc_atom.F +++ b/src/xc/xc_atom.F @@ -337,9 +337,10 @@ CONTAINS !> \param w ... !> \param vxc ... !> \param vxg ... +!> \param do_triplet ... ! ************************************************************************************************** SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, & - deriv_set, w, vxc, vxg) + deriv_set, w, vxc, vxg, do_triplet) ! As input of this routine one gets rho and drho on a one dimensional grid. ! The grid is the angular grid corresponding to a given point ir on the radial grid. @@ -354,12 +355,13 @@ CONTAINS REAL(dp), DIMENSION(:, :), POINTER :: w REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg + LOGICAL, INTENT(IN), OPTIONAL :: do_triplet CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r' INTEGER :: handle, ispin, nspins LOGICAL :: lsd - REAL(dp) :: drho_cutoff + REAL(dp) :: drho_cutoff, my_fac_triplet TYPE(cp_sll_xc_deriv_type), POINTER :: pos TYPE(pw_pool_type), POINTER :: pw_pool TYPE(pw_type), DIMENSION(:), POINTER :: vxc_pw, vxc_tau_pw @@ -373,6 +375,9 @@ CONTAINS IF (ASSOCIATED(rho_set%rhoa)) THEN lsd = .TRUE. END IF + my_fac_triplet = 1.0_dp + IF ((PRESENT(do_triplet)) .AND. (do_triplet)) my_fac_triplet = -1.0_dp + CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff) xc_fun_section => section_vals_get_subs_vals(xc_section, & "XC_FUNCTIONAL") @@ -402,7 +407,7 @@ CONTAINS NULLIFY (vxc_tau_pw) CALL xc_calc_2nd_deriv_analytical(vxc_pw, vxc_tau_pw, deriv_set, rho_set, rho1_set, pw_pool, & - xc_section, gapw=.TRUE., vxg=vxg) + xc_section, gapw=.TRUE., vxg=vxg, tddfpt_fac=my_fac_triplet) DEALLOCATE (vxc_pw) diff --git a/tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp b/tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp new file mode 100644 index 0000000000..6444292ce3 --- /dev/null +++ b/tests/QS/regtest-tddfpt/H2O_GAPW_1_triplet.inp @@ -0,0 +1,67 @@ +&GLOBAL + PROJECT H2O_GAPW + RUN_TYPE ENERGY + PRINT_LEVEL LOW +&END GLOBAL +&FORCE_EVAL + METHOD Quickstep + &PROPERTIES + &TDDFPT + NSTATES 3 + MAX_ITER 10 + MAX_KV 10 + CONVERGENCE 1.0e-5 + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + RKS_TRIPLETS T + &END TDDFPT + &END PROPERTIES + + &DFT + BASIS_SET_FILE_NAME EMSL_BASIS_SETS + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 200 + &END MGRID + &QS + METHOD GAPW + &END QS + &SCF + MAX_SCF 40 + SCF_GUESS ATOMIC + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC NONE + &END CELL + &COORD + O 0.000000 0.000000 -0.065587 H2O + H 0.000000 -0.757136 0.520545 H2O + H 0.000000 0.757136 0.520545 H2O + &END COORD + &TOPOLOGY + &CENTER_COORDINATES + &END + &END + &KIND H + BASIS_SET ORB 6-311Gxx + POTENTIAL ALL + &END KIND + &KIND O + BASIS_SET ORB 6-311Gxx + POTENTIAL ALL + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp b/tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp new file mode 100644 index 0000000000..2494b6a811 --- /dev/null +++ b/tests/QS/regtest-tddfpt/Ne_GAPW_triplet.inp @@ -0,0 +1,76 @@ +&GLOBAL + PROJECT Ne + PRINT_LEVEL LOW + RUN_TYPE ENERGY +&END GLOBAL +&FORCE_EVAL + &PROPERTIES + &TDDFPT + &DIPOLE_MOMENTS + DIPOLE_FORM LENGTH + &END DIPOLE_MOMENTS + NSTATES 1 + MAX_ITER 10 + MAX_KV 60 + CONVERGENCE 1.0e-5 + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + 2ND_DERIV_ANALYTICAL T + &END XC + RKS_TRIPLETS T + &END TDDFPT + &END PROPERTIES + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME EMSL_BASIS_SETS + BASIS_SET_FILE_NAME basis + POTENTIAL_FILE_NAME POTENTIAL + &QS + METHOD GAPW + EPS_DEFAULT 1.0E-17 + EPS_PGF_ORB 1.0E-20 + &END QS + &MGRID + CUTOFF 200 + REL_CUTOFF 80 + &END MGRID + &SCF + SCF_GUESS RESTART + &OT + PRECONDITIONER FULL_SINGLE_INVERSE + MINIMIZER DIIS + STEPSIZE 0.1 + &END + &OUTER_SCF + MAX_SCF 20 + EPS_SCF 1.0E-7 + &END + MAX_SCF 20 + EPS_SCF 1.0E-7 + &END SCF + &POISSON + PERIODIC NONE + PSOLVER WAVELET + &END + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + 2ND_DERIV_ANALYTICAL T + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + PERIODIC NONE + &END CELL + &COORD + Ne 2.5 2.5 2.5 + &END COORD + &KIND Ne + BASIS_SET Ahlrichs-def2-SVP + POTENTIAL ALL + &END KIND + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/QS/regtest-tddfpt/TEST_FILES b/tests/QS/regtest-tddfpt/TEST_FILES index 614903555f..24889f7a4f 100644 --- a/tests/QS/regtest-tddfpt/TEST_FILES +++ b/tests/QS/regtest-tddfpt/TEST_FILES @@ -24,4 +24,6 @@ H2O_GAPW_4.inp 37 4.0E-06 H2O_GAPW_XC_1.inp 37 4.0E-06 0.619451E+00 H2O_GAPW_XC_2.inp 37 4.0E-06 0.812642E+00 H2O_GAPW_XC_3.inp 37 4.0E-06 0.836577E+00 +H2O_GAPW_1_triplet.inp 37 4.0E-06 0.505991E+00 +Ne_GAPW_triplet.inp 37 4.0E-06 0.168075E+01 #EOF