From 2c8265a529330abef610916e7989fde72eb0bdf5 Mon Sep 17 00:00:00 2001 From: Ziwei Chai Date: Wed, 1 Jul 2026 21:18:53 +0800 Subject: [PATCH] Add Andreussi's Solvent Aware SCCS method (#5495) --- src/common/bibliography.F | 16 +- src/cp_control_types.F | 10 +- src/cp_control_utils.F | 59 +++- src/input_constants.F | 3 +- src/input_cp2k_dft.F | 115 ++++++- src/pw/pw_poisson_methods.F | 52 ++- src/qs_sccs.F | 399 ++++++++++++++++++++---- tests/QS/regtest-sasccs/H2_sasccs.inp | 74 +++++ tests/QS/regtest-sasccs/TEST_FILES.toml | 2 + tests/TEST_DIRS | 1 + 10 files changed, 651 insertions(+), 80 deletions(-) create mode 100644 tests/QS/regtest-sasccs/H2_sasccs.inp create mode 100644 tests/QS/regtest-sasccs/TEST_FILES.toml diff --git a/src/common/bibliography.F b/src/common/bibliography.F index c6807b418a..e36fb273ff 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -96,7 +96,8 @@ MODULE bibliography Hanasaki2025, Tan2025, Broyden1965, Johnson1988, Kerker1981, & FuHo1983, MethfesselPaxton1989, Marzari1999, dosSantos2023, Mermin1965, & KuhneHeskeProdan2020, Schreder2021, Schreder2024_1, Schreder2024_2, & - Shiga2022, Lindh1995, Chai2024a, Rullan2026, Sundararaman2017 + Shiga2022, Lindh1995, Chai2024a, Rullan2026, Sundararaman2017, Andreussi2019, & + Chai2025a CONTAINS @@ -2197,6 +2198,19 @@ CONTAINS source="J. Chem. Phys.", volume="146", pages="114104", & year=2017, doi="10.1063/1.4978411") + CALL add_reference(key=Andreussi2019, & + authors=s2a("O. Andreussi", "N. G. Hormann", "F. Nattino", "G. Fisicaro", "S. Goedecker", "N. Marzari"), & + title="Solvent-aware interfaces in continuum solvation", & + source="J. Chem. Theory Comput.", volume="15", pages="1996-2009", & + year=2019, doi="10.1021/acs.jctc.8b01174") + + CALL add_reference(key=Chai2025a, & + authors=s2a("Z. Chai", "S. Luber"), & + title="Functional analytic derivation and CP2K implementation of the SCCS model based "// & + "on the solvent-aware interface", & + source="Comput. Phys. Commun.", volume="311", pages="109563", & + year=2025, doi="10.1016/j.cpc.2025.109563") + END SUBROUTINE add_all_references END MODULE bibliography diff --git a/src/cp_control_types.F b/src/cp_control_types.F index d57974958f..47d3cb3bf7 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -554,7 +554,8 @@ MODULE cp_control_types ! \brief Control parameters for the SCCS models ! ************************************************************************************************** TYPE sccs_control_type - LOGICAL :: sccs_activated = .FALSE. + LOGICAL :: sccs_activated = .FALSE., & + saa_activated = .FALSE. INTEGER :: derivative_method = 0, & max_iter = 0, & method_id = 0 @@ -569,7 +570,12 @@ MODULE cp_control_types mixing = 0.0_dp, & rho_zero = 0.0_dp, & rho_max = 0.0_dp, & - rho_min = 0.0_dp + rho_min = 0.0_dp, & + f0 = 0.0_dp, & + delta_eta = 0.0_dp, & + alpha_zeta = 0.0_dp, & + delta_zeta = 0.0_dp, & + R_solv = 0.0_dp END TYPE sccs_control_type ! ************************************************************************************************** diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index c90741ffab..3a2478e19e 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -10,11 +10,11 @@ ! ************************************************************************************************** MODULE cp_control_utils USE bibliography, ONLY: & - Andreussi2012, Dewar1977, Dewar1985, Elstner1998, Fattebert2002, Grimme2017, Hu2007, & - Katbashev2025, Krack2000, Lippert1997, Lippert1999, Porezag1995, Pracht2019, Repasky2002, & - Rocha2006, Schenter2008, Seifert1996, Souza2002, Stengel2009, Stewart1989, Stewart2007, & - Thiel1992, Umari2002, VanVoorhis2015, VandeVondele2005a, VandeVondele2005b, Yin2017, & - Zhechkov2005, cite_reference + Andreussi2012, Andreussi2019, Chai2025a, Dewar1977, Dewar1985, Elstner1998, Fattebert2002, & + Grimme2017, Hu2007, Katbashev2025, Krack2000, Lippert1997, Lippert1999, Porezag1995, & + Pracht2019, Repasky2002, Rocha2006, Schenter2008, Seifert1996, Souza2002, Stengel2009, & + Stewart1989, Stewart2007, Thiel1992, Umari2002, VanVoorhis2015, VandeVondele2005a, & + VandeVondele2005b, Yin2017, Zhechkov2005, cite_reference USE cell_types, ONLY: cell_transform_input_cartesian,& cell_type USE cp_control_types, ONLY: & @@ -58,8 +58,8 @@ MODULE cp_control_utils gaussian_env, gfn1xtb, gfn_tblite, kg_tnadd_embed, kg_tnadd_embed_ri, no_admm_type, & numerical, ramp_env, real_time_propagation, rtp_method_bse, sccs_andreussi, & sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, & - sccs_fattebert_gygi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, & - sic_mauri_us, sic_none, slater, tblite_cli_born_kernel_auto, & + sccs_fattebert_gygi, sccs_saa_andreussi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, & + sic_mauri_spz, sic_mauri_us, sic_none, slater, tblite_cli_born_kernel_auto, & tblite_cli_solution_state_gsolv, tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, & tblite_cli_solvation_gb, tblite_cli_solvation_gbe, tblite_cli_solvation_gbsa, & tblite_guess_ceh, tblite_mixer_memory_inherit, tblite_scc_mixer_auto, & @@ -742,6 +742,40 @@ CONTAINS CALL section_vals_val_get(tmp_section, "RHO_ZERO", & r_val=dft_control%sccs_control%rho_zero) CALL cite_reference(Fattebert2002) + CASE (sccs_saa_andreussi) + tmp_section => section_vals_get_subs_vals(sccs_section, "SAA_ANDREUSSI") + CALL section_vals_get(tmp_section, explicit=is_present) + IF (.NOT. is_present) THEN + CALL cp_abort(__LOCATION__, & + "SCCS method SAA_ANDREUSSI requires the "// & + "SAA_ANDREUSSI section.") + END IF + IF (.NOT. ALL(cell%perd == 1)) THEN + CALL cp_abort(__LOCATION__, & + "SCCS method SAA_ANDREUSSI is only implemented for "// & + "3D periodic calculations.") + END IF + CALL section_vals_val_get(tmp_section, "RHO_MAX", & + r_val=dft_control%sccs_control%rho_max) + CALL section_vals_val_get(tmp_section, "RHO_MIN", & + r_val=dft_control%sccs_control%rho_min) + IF (dft_control%sccs_control%rho_max < dft_control%sccs_control%rho_min) THEN + CALL cp_abort(__LOCATION__, & + "The SCCS parameter RHO_MAX is smaller than RHO_MIN. "// & + "Please, check your input!") + END IF + CALL section_vals_val_get(tmp_section, "F0", & + r_val=dft_control%sccs_control%f0) + CALL section_vals_val_get(tmp_section, "DELTA_ETA", & + r_val=dft_control%sccs_control%delta_eta) + CALL section_vals_val_get(tmp_section, "DELTA_ZETA", & + r_val=dft_control%sccs_control%delta_zeta) + CALL section_vals_val_get(tmp_section, "R_SOLV", & + r_val=dft_control%sccs_control%r_solv) + CALL section_vals_val_get(tmp_section, "ALPHA_ZETA", & + r_val=dft_control%sccs_control%alpha_zeta) + CALL cite_reference(Andreussi2019) + CALL cite_reference(Chai2025a) CASE DEFAULT CPABORT("Invalid SCCS model specified. Please, check your input!") END SELECT @@ -2455,6 +2489,17 @@ CONTAINS "SCCS| Dielectric function proposed by Fattebert and Gygi", & "SCCS| beta", dft_control%sccs_control%beta, & "SCCS| rho_zero", dft_control%sccs_control%rho_zero + CASE (sccs_saa_andreussi) + WRITE (UNIT=output_unit, FMT="(T2,A,/,A,/,(T2,A,T61,ES20.6))") & + "SCCS| Dielectric function of the solvent aware algorithm", & + "SCCS| proposed by Andreussi et al.", & + "SCCS| rho_max", dft_control%sccs_control%rho_max, & + "SCCS| rho_min", dft_control%sccs_control%rho_min, & + "SCCS| f0", dft_control%sccs_control%f0, & + "SCCS| delta_eta", dft_control%sccs_control%delta_eta, & + "SCCS| alpha_zeta", dft_control%sccs_control%alpha_zeta, & + "SCCS| delta_zeta", dft_control%sccs_control%delta_zeta, & + "SCCS| R_solv", dft_control%sccs_control%R_solv CASE DEFAULT CPABORT("Invalid SCCS model specified. Please, check your input!") END SELECT diff --git a/src/input_constants.F b/src/input_constants.F index 26ba2b29f6..fdceaa510e 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1219,7 +1219,8 @@ MODULE input_constants ! Self-consistent continuum solvation (SCCS) models INTEGER, PARAMETER, PUBLIC :: sccs_andreussi = 0, & - sccs_fattebert_gygi = 1 + sccs_fattebert_gygi = 1, & + sccs_saa_andreussi = 2 ! Derivative method used by the SCCS INTEGER, PARAMETER, PUBLIC :: sccs_derivative_fft = 0, & diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 7b2cbd4984..082cc27f4b 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -15,9 +15,9 @@ MODULE input_cp2k_dft USE basis_set_types, ONLY: basis_sort_default,& basis_sort_zet USE bibliography, ONLY: & - Andermatt2016, Andreussi2012, Avezac2005, Bengtsson1999, Blochl1995, Brelaz1979, & - Chai2024a, Fattebert2002, Guidon2010, Iannuzzi2006, Kunert2003, Marek2025, Merlot2014, & - Perdew1981, VandeVondele2005b, Yin2017 + Andermatt2016, Andreussi2012, Andreussi2019, Avezac2005, Bengtsson1999, Blochl1995, & + Brelaz1979, Chai2024a, Chai2025a, Fattebert2002, Guidon2010, Iannuzzi2006, Kunert2003, & + Marek2025, Merlot2014, Perdew1981, VandeVondele2005b, Yin2017 USE cp_output_handling, ONLY: add_last_numeric,& cp_print_key_section_create,& high_print_level,& @@ -46,10 +46,10 @@ MODULE input_cp2k_dft rel_trans_full, rel_trans_molecule, rel_zora, rel_zora_full, rel_zora_mp, & rtp_bse_ham_g0w0, rtp_bse_ham_ks, rtp_method_bse, rtp_method_tddft, sccs_andreussi, & sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, & - sccs_fattebert_gygi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, & - sic_mauri_us, sic_none, slater, use_mom_ref_coac, use_mom_ref_com, use_mom_ref_user, & - use_mom_ref_zero, use_restart_wfn, use_rt_restart, use_scf_wfn, weight_type_mass, & - weight_type_unit + sccs_fattebert_gygi, sccs_saa_andreussi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, & + sic_mauri_spz, sic_mauri_us, sic_none, slater, use_mom_ref_coac, use_mom_ref_com, & + use_mom_ref_user, use_mom_ref_zero, use_restart_wfn, use_rt_restart, use_scf_wfn, & + weight_type_mass, weight_type_unit USE input_cp2k_almo, ONLY: create_almo_scf_section USE input_cp2k_as, ONLY: create_active_space_section USE input_cp2k_ec, ONLY: create_ec_section @@ -2192,7 +2192,7 @@ CONTAINS description="Define the parameters for self-consistent continuum solvation (SCCS) model", & citations=[Fattebert2002, Andreussi2012, Yin2017], & n_keywords=8, & - n_subsections=2, & + n_subsections=3, & repeats=.FALSE.) NULLIFY (keyword) @@ -2330,10 +2330,11 @@ CONTAINS description="Method used for the smoothing of the dielectric function", & usage="METHOD Fattebert-Gygi", & default_i_val=sccs_andreussi, & - enum_c_vals=s2a("ANDREUSSI", "FATTEBERT-GYGI"), & - enum_i_vals=[sccs_andreussi, sccs_fattebert_gygi], & + enum_c_vals=s2a("ANDREUSSI", "FATTEBERT-GYGI", "SAA_ANDREUSSI"), & + enum_i_vals=[sccs_andreussi, sccs_fattebert_gygi, sccs_saa_andreussi], & enum_desc=s2a("Smoothing function proposed by Andreussi et al.", & - "Smoothing function proposed by Fattebert and Gygi")) + "Smoothing function proposed by Fattebert and Gygi", & + "Smoothing function of the solvent aware algorithm")) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -2421,6 +2422,98 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL section_create(subsection, __LOCATION__, & + name="SAA_ANDREUSSI", & + description="Define the parameters of the dielectric smoothing function of "// & + "the solvent aware algorithm proposed by Andreussi et al.", & + citations=[Andreussi2019, Chai2025a], & + n_keywords=7, & + n_subsections=0, & + repeats=.FALSE.) + + CALL keyword_create(keyword, __LOCATION__, & + name="RHO_MAX", & + description="Maximum density value used for the smoothing of the dielectric function", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=0.0035_dp, & + usage="RHO_MAX 0.01") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="RHO_MIN", & + description="Minimum density value used for the smoothing of the dielectric function", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=0.0001_dp, & + usage="RHO_MIN 0.0003") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="F0", & + description="The threshold of the filled fraction that controls whether a point in space "// & + "needs to be removed from the continuum", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=0.65_dp, & + usage="F0 0.65") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="DELTA_ETA", & + description="Controls the softness of function t(r)", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=0.02_dp, & + usage="DELTA_ETA 0.02") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="DELTA_ZETA", & + description="Controls the softness of function u(r)", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + unit_str='bohr', & + default_r_val=0.5_dp, & + usage="DELTA_ZETA 0.5") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="ALPHA_ZETA", & + description="The scaling factor for R_SOLV", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + default_r_val=2.0_dp, & + usage="ALPHA_ZETA 2.0") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, & + name="R_SOLV", & + description="The solvent radius", & + repeats=.FALSE., & + n_var=1, & + type_of_var=real_t, & + unit_str='bohr', & + default_r_val=2.6_dp, & + usage="R_SOLV 2.6") + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_sccs_section ! ************************************************************************************************** diff --git a/src/pw/pw_poisson_methods.F b/src/pw/pw_poisson_methods.F index 68a9f2ece7..02f41935f0 100644 --- a/src/pw/pw_poisson_methods.F +++ b/src/pw/pw_poisson_methods.F @@ -65,7 +65,7 @@ MODULE pw_poisson_methods CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_poisson_methods' PUBLIC :: pw_poisson_rebuild, & - pw_poisson_solve, pw_poisson_set + pw_poisson_solve, pw_poisson_set, pw_func_u_convolution INTEGER, PARAMETER :: use_rs_grid = 0, & use_gs_grid = 1 @@ -1350,4 +1350,54 @@ CONTAINS END SUBROUTINE pw_poisson_set +! ************************************************************************************************** +!> \brief computes the convolution f * u using reciprocal-space multiplication +!> \par +!> Given two real-space functions f and u, this routine evaluates their +!> convolution by transforming both to G space, multiplying their Fourier +!> components, and transforming the product back to real space. +!> \param poisson_env poisson environment +!> \param func input function f in real space +!> \param convolution convolution result (f * u) in real space +!> \param u input function/kernel u in real space +!> \author Ziwei Chai +!> \date 12.03.2026 +! ************************************************************************************************** + SUBROUTINE pw_func_u_convolution(poisson_env, func, convolution, u) + + TYPE(pw_poisson_type), POINTER :: poisson_env + TYPE(pw_r3d_rs_type), INTENT(IN) :: func, u + TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: convolution + CHARACTER(len=*), PARAMETER :: routineN = 'pw_func_u_convolution' + INTEGER :: handle + TYPE(pw_grid_type), POINTER :: pw_grid + TYPE(pw_pool_type), POINTER :: pw_pool + TYPE(pw_c1d_gs_type) :: funcg, ug + + CALL timeset(routineN, handle) + + CPASSERT(ASSOCIATED(poisson_env)) + CALL pw_poisson_rebuild(poisson_env, func) + + ! point pw + pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool + pw_grid => pw_pool%pw_grid + ! func and u in G space + CALL pw_pool%create_pw(funcg) + CALL pw_pool%create_pw(ug) + + CALL pw_transfer(func, funcg) + CALL pw_transfer(u, ug) + + funcg%array(:) = funcg%array(:)*ug%array(:) + + CALL pw_transfer(funcg, convolution) + + CALL pw_pool%give_back_pw(funcg) + CALL pw_pool%give_back_pw(ug) + + CALL timestop(handle) + + END SUBROUTINE pw_func_u_convolution + END MODULE pw_poisson_methods diff --git a/src/qs_sccs.F b/src/qs_sccs.F index 535b461e42..639207d7b8 100644 --- a/src/qs_sccs.F +++ b/src/qs_sccs.F @@ -47,7 +47,8 @@ MODULE qs_sccs sccs_derivative_cd5,& sccs_derivative_cd7,& sccs_derivative_fft,& - sccs_fattebert_gygi + sccs_fattebert_gygi,& + sccs_saa_andreussi USE input_section_types, ONLY: section_get_ivals,& section_get_lval,& section_vals_get_subs_vals,& @@ -57,6 +58,7 @@ MODULE qs_sccs dp,& int_8 USE mathconstants, ONLY: fourpi,& + pi,& twopi USE message_passing, ONLY: mp_para_env_type USE particle_list_types, ONLY: particle_list_type @@ -67,9 +69,11 @@ MODULE qs_sccs pw_derive,& pw_integral_ab,& pw_integrate_function,& + pw_scale,& pw_transfer,& pw_zero - USE pw_poisson_methods, ONLY: pw_poisson_solve + USE pw_poisson_methods, ONLY: pw_func_u_convolution,& + pw_poisson_solve USE pw_poisson_types, ONLY: pw_poisson_analytic,& pw_poisson_mt,& pw_poisson_type @@ -138,10 +142,10 @@ CONTAINS INTEGER, DIMENSION(3) :: lb, ub LOGICAL :: append_cube, calculate_stress_tensor, & do_kpoints, mpi_io, should_output - REAL(KIND=dp) :: cavity_surface, cavity_volume, cell_volume, dphi2, dvol, e_tot, & - epsilon_solvent, f, polarisation_charge, rho_delta, rho_delta_avg, rho_delta_max, & - rho_iter_new, tot_rho_elec, tot_rho_solute - REAL(KIND=dp), DIMENSION(3) :: abc + REAL(KIND=dp) :: alpha_zeta, cavity_surface, cavity_volume, cell_volume, delta_eta, & + delta_zeta, dphi2, dvol, e_tot, epsilon_solvent, f, f0, polarisation_charge, r, R_solv, & + rho_delta, rho_delta_avg, rho_delta_max, rho_iter_new, tot_rho_elec, tot_rho_solute + REAL(KIND=dp), DIMENSION(3) :: abc, lxyz, uxyz TYPE(cp_logger_type), POINTER :: logger TYPE(cp_subsys_type), POINTER :: cp_subsys TYPE(dft_control_type), POINTER :: dft_control @@ -151,7 +155,9 @@ CONTAINS TYPE(pw_poisson_type), POINTER :: poisson_env TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools TYPE(pw_pool_type), POINTER :: auxbas_pw_pool - TYPE(pw_r3d_rs_type) :: deps_elec, eps_elec + TYPE(pw_r3d_rs_type) :: d_s_rhoel, deps_elec, dtf, eps_elec, ff, & + p_e_interface, s, t, tem_convolution, & + tem_func, u TYPE(pw_r3d_rs_type), DIMENSION(3) :: dln_eps_elec, dphi_tot TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_pw_r TYPE(pw_r3d_rs_type), POINTER :: rho_pw_r_sccs @@ -181,8 +187,8 @@ CONTAINS ! Load data from Quickstep environment CALL get_qs_env(qs_env=qs_env, & cp_subsys=cp_subsys, & - dft_control=dft_control, & do_kpoints=do_kpoints, & + dft_control=dft_control, & energy=energy, & input=input, & para_env=para_env, & @@ -290,6 +296,20 @@ CONTAINS ! Calculate the dielectric (smoothed) function of rho_elec in r-space CALL auxbas_pw_pool%create_pw(eps_elec) CALL auxbas_pw_pool%create_pw(deps_elec) + IF (sccs_control%method_id == sccs_saa_andreussi) THEN + CALL auxbas_pw_pool%create_pw(s) + CALL auxbas_pw_pool%create_pw(d_s_rhoel) + CALL auxbas_pw_pool%create_pw(u) + CALL auxbas_pw_pool%create_pw(ff) + CALL auxbas_pw_pool%create_pw(t) + CALL auxbas_pw_pool%create_pw(dtf) + f0 = sccs_control%f0 + delta_eta = sccs_control%delta_eta + alpha_zeta = sccs_control%alpha_zeta + R_solv = sccs_control%R_solv + delta_zeta = sccs_control%delta_zeta + END IF + ! Relative permittivity or dielectric constant of the solvent (medium) epsilon_solvent = sccs_control%epsilon_solvent SELECT CASE (sccs_control%method_id) @@ -299,37 +319,42 @@ CONTAINS CASE (sccs_fattebert_gygi) CALL fattebert_gygi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%beta, & sccs_control%rho_zero) + CASE (sccs_saa_andreussi) + CALL sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, sccs_control%rho_max, & + sccs_control%rho_min, s, d_s_rhoel) CASE DEFAULT CPABORT("Invalid method specified for SCCS model") END SELECT - ! Optional output of the dielectric function in cube file format - filename = "DIELECTRIC_FUNCTION" - cube_path = TRIM(print_path)//"%"//TRIM(filename) - IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), & - cp_p_file)) THEN - append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND") - my_pos_cube = "REWIND" - IF (append_cube) my_pos_cube = "APPEND" - mpi_io = .TRUE. - cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), & - extension=".cube", middle_name=TRIM(filename), & - file_position=my_pos_cube, log_filename=.FALSE., & - mpi_io=mpi_io, fout=mpi_filename) - IF (output_unit > 0) THEN - IF (.NOT. mpi_io) THEN - INQUIRE (UNIT=cube_unit, NAME=filename) - ELSE - filename = mpi_filename + IF (sccs_control%method_id /= sccs_saa_andreussi) THEN + ! Optional output of the dielectric function in cube file format + filename = "DIELECTRIC_FUNCTION" + cube_path = TRIM(print_path)//"%"//TRIM(filename) + IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), & + cp_p_file)) THEN + append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND") + my_pos_cube = "REWIND" + IF (append_cube) my_pos_cube = "APPEND" + mpi_io = .TRUE. + cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., & + mpi_io=mpi_io, fout=mpi_filename) + IF (output_unit > 0) THEN + IF (.NOT. mpi_io) THEN + INQUIRE (UNIT=cube_unit, NAME=filename) + ELSE + filename = mpi_filename + END IF + WRITE (UNIT=output_unit, FMT="(T3,A)") & + "SCCS| The dielectric function is written in cube file format to the file:", & + "SCCS| "//TRIM(filename) END IF - WRITE (UNIT=output_unit, FMT="(T3,A)") & - "SCCS| The dielectric function is written in cube file format to the file:", & - "SCCS| "//TRIM(filename) + CALL cp_pw_to_cube(eps_elec, cube_unit, TRIM(filename), particles=particles, & + stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), & + mpi_io=mpi_io) + CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) END IF - CALL cp_pw_to_cube(eps_elec, cube_unit, TRIM(filename), particles=particles, & - stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), & - mpi_io=mpi_io) - CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) END IF ! Calculate the (quantum) volume and surface of the solute cavity @@ -425,6 +450,7 @@ CONTAINS CALL surface_fattebert_gygi(rho_elec, norm_drho_elec, theta, epsilon_solvent, & sccs_control%beta, sccs_control%rho_zero, & sccs_control%delta_rho) + CASE (sccs_saa_andreussi) CASE DEFAULT CPABORT("Invalid method specified for SCCS model") END SELECT @@ -466,35 +492,136 @@ CONTAINS ! Reassign work storage to rho_tot_zero, because rho_elec is no longer needed CALL auxbas_pw_pool%create_pw(rho_tot_zero) - ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space - ! eps_elec <- ln(eps_elec) + IF (sccs_control%method_id /= sccs_saa_andreussi) THEN + + ! Build the initial (rho_iter = 0) total charge density (solute plus polarisation) in r-space + ! eps_elec <- ln(eps_elec) !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,k) & !$OMP SHARED(eps_elec,lb,message,output_unit,para_env,ub) & !$OMP SHARED(rho_solute,rho_tot_zero) - DO k = lb(3), ub(3) - DO j = lb(2), ub(2) - DO i = lb(1), ub(1) - IF (eps_elec%array(i, j, k) < 1.0_dp) THEN - WRITE (UNIT=message, FMT="(A,ES12.3,A,3(I0,A))") & - "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), & - " encountered at grid point (", i, ",", j, ",", k, ")" - CPABORT(message) - END IF - rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k) - eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k)) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + IF (eps_elec%array(i, j, k) < 1.0_dp) THEN + WRITE (UNIT=message, FMT="(A,ES12.3,A,3(I0,A))") & + "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), & + " encountered at grid point (", i, ",", j, ",", k, ")" + CPABORT(message) + END IF + rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k) + eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k)) + END DO END DO END DO - END DO !$OMP END PARALLEL DO + ELSE + + lxyz(:) = REAL(u%pw_grid%bounds(1, :), dp) + uxyz(:) = REAL(u%pw_grid%bounds(2, :), dp) +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k,r) & +!$OMP SHARED(lb,message,output_unit,para_env,ub,dtf,uxyz,lxyz) & +!$OMP SHARED(u, alpha_zeta, R_solv, delta_zeta, eps_elec) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + IF (eps_elec%array(i, j, k) < 1.0_dp) THEN + WRITE (UNIT=message, FMT="(A,ES12.3,A,3(I0,A))") & + "SCCS| Invalid dielectric function value ", eps_elec%array(i, j, k), & + " encountered at grid point (", i, ",", j, ",", k, ")" + CPABORT(message) + END IF + IF (i > -1) r = ((REAL(uxyz(1), dp) - REAL(i, dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp + IF (i < 0) r = ((REAL(i, dp) - REAL(lxyz(1), dp) + 0.5_dp)*u%pw_grid%dr(1))**2.0_dp + IF (j > -1) r = ((REAL(uxyz(2), dp) - REAL(j, dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r + IF (j < 0) r = ((REAL(j, dp) - REAL(lxyz(2), dp) + 0.5_dp)*u%pw_grid%dr(2))**2.0_dp + r + IF (k > -1) r = ((REAL(uxyz(3), dp) - REAL(k, dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r + IF (k < 0) r = ((REAL(k, dp) - REAL(lxyz(3), dp) + 0.5_dp)*u%pw_grid%dr(3))**2.0_dp + r + r = SQRT(r) + u%array(i, j, k) = (1.0_dp/2.0_dp)*erfc((r - alpha_zeta*R_solv)/delta_zeta) + dtf%array(i, j, k) = 1.0_dp + END DO + END DO + END DO +!$OMP END PARALLEL DO + + CALL pw_func_u_convolution(poisson_env=poisson_env, func=dtf, convolution=ff, u=u) + CALL pw_scale(u, 1.0_dp/ff%array(0, 0, 0)) + CALL pw_func_u_convolution(poisson_env=poisson_env, func=s, convolution=ff, u=u) + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k) & +!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) & +!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + t%array(i, j, k) = 0.5_dp*(1.0_dp + erf((ff%array(i, j, k) - f0)/delta_eta)) + dtf%array(i, j, k) = (1.0_dp/SQRT(pi))*EXP(-(ff%array(i, j, k) - f0)**2.0/(delta_eta**2.0)) + eps_elec%array(i, j, k) = EXP(LOG(epsilon_solvent)*(1 - (s%array(i, j, k) & + + (1 - s%array(i, j, k))*t%array(i, j, k)))) + rho_tot_zero%array(i, j, k) = rho_solute%array(i, j, k)/eps_elec%array(i, j, k) +! eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k)) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + ! Optional output of the dielectric function in cube file format (solvent aware algorithm case) + filename = "DIELECTRIC_FUNCTION" + cube_path = TRIM(print_path)//"%"//TRIM(filename) + IF (BTEST(cp_print_key_should_output(logger%iter_info, input, TRIM(cube_path)), & + cp_p_file)) THEN + append_cube = section_get_lval(input, TRIM(cube_path)//"%APPEND") + my_pos_cube = "REWIND" + IF (append_cube) my_pos_cube = "APPEND" + mpi_io = .TRUE. + cube_unit = cp_print_key_unit_nr(logger, input, TRIM(cube_path), & + extension=".cube", middle_name=TRIM(filename), & + file_position=my_pos_cube, log_filename=.FALSE., & + mpi_io=mpi_io, fout=mpi_filename) + IF (output_unit > 0) THEN + IF (.NOT. mpi_io) THEN + INQUIRE (UNIT=cube_unit, NAME=filename) + ELSE + filename = mpi_filename + END IF + WRITE (UNIT=output_unit, FMT="(T3,A)") & + "SCCS| The dielectric function is written in cube file format to the file:", & + "SCCS| "//TRIM(filename) + END IF + CALL cp_pw_to_cube(eps_elec, cube_unit, TRIM(filename), particles=particles, & + stride=section_get_ivals(input, TRIM(cube_path)//"%STRIDE"), & + mpi_io=mpi_io) + CALL cp_print_key_finished_output(cube_unit, logger, input, TRIM(cube_path), mpi_io=mpi_io) + END IF + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k) & +!$OMP SHARED(lb,para_env,ub, delta_eta, f0, epsilon_solvent) & +!$OMP SHARED(ff, t, dtf, eps_elec, rho_tot_zero, rho_solute, s) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + eps_elec%array(i, j, k) = LOG(eps_elec%array(i, j, k)) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + END IF + ! Build the derivative of LOG(eps_elec) DO i = 1, 3 CALL auxbas_pw_pool%create_pw(dln_eps_elec(i)) CALL pw_zero(dln_eps_elec(i)) END DO CALL derive(eps_elec, dln_eps_elec, sccs_control%derivative_method, pw_env, input) - CALL auxbas_pw_pool%give_back_pw(eps_elec) + IF (sccs_control%method_id /= sccs_saa_andreussi) THEN + CALL auxbas_pw_pool%give_back_pw(eps_elec) + END IF ! Print header for the SCCS cycle IF (should_output .AND. (output_unit > 0)) THEN @@ -534,6 +661,11 @@ CONTAINS CALL auxbas_pw_pool%create_pw(phi_tot) CALL pw_zero(phi_tot) + IF (sccs_control%method_id == sccs_saa_andreussi) THEN + CALL auxbas_pw_pool%create_pw(p_e_interface) + CALL pw_zero(p_e_interface) + END IF + ! Main SCCS iteration loop iter = 0 @@ -620,6 +752,46 @@ CONTAINS END DO iter_loop + IF (sccs_control%method_id == sccs_saa_andreussi) THEN + CALL auxbas_pw_pool%create_pw(tem_convolution) + CALL auxbas_pw_pool%create_pw(tem_func) + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k,dphi2) & +!$OMP SHARED(lb,ub, tem_func, s, dtf) & +!$OMP SHARED(dphi_tot, p_e_interface, eps_elec, epsilon_solvent) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + p_e_interface%array(i, j, k) = 1.0_dp/(4.0_dp*twopi)*(EXP(eps_elec%array(i, j, k))*LOG(epsilon_solvent)) + dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + & + dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + & + dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k) + p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*dphi2 + tem_func%array(i, j, k) = (1.0_dp - s%array(i, j, k))*p_e_interface%array(i, j, k)*dtf%array(i, j, k) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + CALL pw_func_u_convolution(poisson_env=poisson_env, func=tem_func, convolution=tem_convolution, u=u) + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k) & +!$OMP SHARED(lb,ub) & +!$OMP SHARED(p_e_interface, t, tem_convolution, phi_tot, d_s_rhoel) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + p_e_interface%array(i, j, k) = p_e_interface%array(i, j, k)*(1.0_dp - t%array(i, j, k)) & + + tem_convolution%array(i, j, k) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + END IF + ! Release work storage which is no longer needed CALL auxbas_pw_pool%give_back_pw(rho_tot_zero) DO i = 1, 3 @@ -756,6 +928,16 @@ CONTAINS CALL auxbas_pw_pool%give_back_pw(rho_solute) CALL auxbas_pw_pool%give_back_pw(phi_tot) CALL auxbas_pw_pool%give_back_pw(rho_tot) + IF (sccs_control%method_id == sccs_saa_andreussi) THEN + CALL auxbas_pw_pool%give_back_pw(s) + CALL auxbas_pw_pool%give_back_pw(u) + CALL auxbas_pw_pool%give_back_pw(ff) + CALL auxbas_pw_pool%give_back_pw(t) + CALL auxbas_pw_pool%give_back_pw(dtf) + CALL auxbas_pw_pool%give_back_pw(tem_convolution) + CALL auxbas_pw_pool%give_back_pw(tem_func) + CALL auxbas_pw_pool%give_back_pw(eps_elec) + END IF END BLOCK ! Calculate additional solvation terms @@ -777,23 +959,45 @@ CONTAINS END IF ! Calculate SCCS contribution to the Kohn-Sham potential - f = -0.5_dp*dvol/fourpi + IF (sccs_control%method_id /= sccs_saa_andreussi) THEN + + f = -0.5_dp*dvol/fourpi !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(dphi2,i,j,k) & !$OMP SHARED(f,deps_elec,dphi_tot) & !$OMP SHARED(lb,ub,v_sccs) - DO k = lb(3), ub(3) - DO j = lb(2), ub(2) - DO i = lb(1), ub(1) - dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + & - dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + & - dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k) - v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2 + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + dphi2 = dphi_tot(1)%array(i, j, k)*dphi_tot(1)%array(i, j, k) + & + dphi_tot(2)%array(i, j, k)*dphi_tot(2)%array(i, j, k) + & + dphi_tot(3)%array(i, j, k)*dphi_tot(3)%array(i, j, k) + v_sccs%array(i, j, k) = v_sccs%array(i, j, k) + f*deps_elec%array(i, j, k)*dphi2 + END DO END DO END DO - END DO !$OMP END PARALLEL DO + ELSE + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(i,j,k) & +!$OMP SHARED(lb,ub, dvol) & +!$OMP SHARED(p_E_interface, v_sccs, d_s_rhoel) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + v_sccs%array(i, j, k) = dvol*d_s_rhoel%array(i, j, k)*p_E_interface%array(i, j, k) + END DO + END DO + END DO +!$OMP END PARALLEL DO + + CALL auxbas_pw_pool%give_back_pw(d_s_rhoel) + CALL auxbas_pw_pool%give_back_pw(p_E_interface) + + END IF + CALL auxbas_pw_pool%give_back_pw(deps_elec) DO i = 1, 3 CALL auxbas_pw_pool%give_back_pw(dphi_tot(i)) @@ -946,6 +1150,87 @@ CONTAINS END SUBROUTINE fattebert_gygi +! ************************************************************************************************** +!> \brief Calculate the smoothed dielectric function of solvent aware algorithm +!> \param rho_elec ... +!> \param eps_elec ... +!> \param deps_elec ... +!> \param epsilon_solvent ... +!> \param rho_max ... +!> \param rho_min ... +!> \param s ... +!> \param d_s_rhoel ... +!> \author Ziwei Chai +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE sa_andreussi(rho_elec, eps_elec, deps_elec, epsilon_solvent, rho_max, & + rho_min, s, d_s_rhoel) + + TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_elec, eps_elec, deps_elec + REAL(KIND=dp), INTENT(IN) :: epsilon_solvent, rho_max, rho_min + TYPE(pw_r3d_rs_type), INTENT(IN) :: s, d_s_rhoel + + CHARACTER(LEN=*), PARAMETER :: routineN = 'sa_andreussi' + REAL(KIND=dp), PARAMETER :: rhotol = 1.0E-12_dp + + INTEGER :: handle, i, j, k + INTEGER, DIMENSION(3) :: lb, ub + REAL(KIND=dp) :: diff, dq, dt, f, ln_rho_max, ln_rho_min, & + q, rho, t, x, y + + CALL timeset(routineN, handle) + + f = LOG(epsilon_solvent)/twopi + diff = rho_max - rho_min + IF (diff < SQRT(rhotol)) CPABORT("SCCS: Difference between rho(min) and rho(max) is too small") + IF (rho_min >= rhotol) THEN + ln_rho_max = LOG(rho_max) + ln_rho_min = LOG(rho_min) + q = twopi/(ln_rho_max - ln_rho_min) + dq = -f*q + END IF + + lb(1:3) = rho_elec%pw_grid%bounds_local(1, 1:3) + ub(1:3) = rho_elec%pw_grid%bounds_local(2, 1:3) + + ! Calculate the dielectric function and its derivative +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP PRIVATE(dt,i,j,k,rho,t,x,y) & +!$OMP SHARED(deps_elec,dq,eps_elec,epsilon_solvent,f,lb,ub) & +!$OMP SHARED(ln_rho_max,rho_elec,q,rho_max,rho_min, s, d_s_rhoel) + DO k = lb(3), ub(3) + DO j = lb(2), ub(2) + DO i = lb(1), ub(1) + rho = rho_elec%array(i, j, k) + IF (rho < rho_min) THEN + eps_elec%array(i, j, k) = epsilon_solvent + s%array(i, j, k) = 0.0_dp + deps_elec%array(i, j, k) = 0.0_dp + d_s_rhoel%array(i, j, k) = 0.0_dp + ELSE IF (rho <= rho_max) THEN + x = LOG(rho) + y = q*(ln_rho_max - x) + t = f*(y - SIN(y)) + s%array(i, j, k) = 1 - t/(f*twopi) + eps_elec%array(i, j, k) = EXP(t) + dt = dq*(1.0_dp - COS(y)) + deps_elec%array(i, j, k) = eps_elec%array(i, j, k)*dt/rho + d_s_rhoel%array(i, j, k) = deps_elec%array(i, j, k)/(-eps_elec%array(i, j, k)*LOG(epsilon_solvent)) + ELSE + eps_elec%array(i, j, k) = 1.0_dp + deps_elec%array(i, j, k) = 0.0_dp + s%array(i, j, k) = 1.0_dp + d_s_rhoel%array(i, j, k) = 0.0_dp + END IF + END DO + END DO + END DO +!$OMP END PARALLEL DO + + CALL timestop(handle) + + END SUBROUTINE sa_andreussi + ! ************************************************************************************************** !> \brief Build the numerical derivative of a function on realspace grid !> \param f ... diff --git a/tests/QS/regtest-sasccs/H2_sasccs.inp b/tests/QS/regtest-sasccs/H2_sasccs.inp new file mode 100644 index 0000000000..a3929ca41f --- /dev/null +++ b/tests/QS/regtest-sasccs/H2_sasccs.inp @@ -0,0 +1,74 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT H2_sasccs + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_SET + LSD + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 100 + REL_CUTOFF 40 + &END MGRID + &PLANAR_COUNTER_CHARGE ON + DIST_EDGE 1.0 + GAU_C 0.2 + PARALLEL_PLANE XZ + &END PLANAR_COUNTER_CHARGE + &QS + EPS_DEFAULT 1.0E-10 + &END QS + &SCCS ON + ALPHA [N*m^-1] 0.0 + BETA [kbar] 0.0 + DELTA_RHO 2.0E-5 + DERIVATIVE_METHOD FFT + DIELECTRIC_CONSTANT 78.3 + EPS_SCCS 1.0E-7 + EPS_SCF 0.3 + GAMMA [mN/m] 0.0 + MAX_ITER 100 + METHOD ANDREUSSI + MIXING 0.4 + &SAA_ANDREUSSI + ALPHA_ZETA 2.0 + DELTA_ETA 0.02 + DELTA_ZETA [Bohr] 0.5 + F0 0.65 + RHO_MAX 0.1 + RHO_MIN 0.01 + R_SOLV [Bohr] 2.6 + &END SAA_ANDREUSSI + &END SCCS + &SCF + EPS_SCF 1.0E-6 + MAX_SCF 50 + SCF_GUESS ATOMIC + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_SINGLE_INVERSE + &END OT + &END SCF + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 8.0 9.0 10.0 + &END CELL + &COORD + H 3.60 4.50 5.00 + H 4.40 4.50 5.00 + &END COORD + &KIND H + BASIS_SET DZV-GTH-PADE + POTENTIAL GTH-PADE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-sasccs/TEST_FILES.toml b/tests/QS/regtest-sasccs/TEST_FILES.toml new file mode 100644 index 0000000000..2d489338d7 --- /dev/null +++ b/tests/QS/regtest-sasccs/TEST_FILES.toml @@ -0,0 +1,2 @@ +# runs are executed in the same order as in this file +"H2_sasccs.inp" = [{matcher="E_total", tol=1e-7, ref=-1.114289814572276}] diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 464a34eb85..1a0e56c8f3 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -13,6 +13,7 @@ QS/regtest-dft-d3bj-lib-ref libint s_dftd3 QS/regtest-pcc QS/regtest-new-pulay-mixing QS/regtest-gcscf +QS/regtest-sasccs QS/regtest-kubo-transport QS/regtest-gw-realspace libint !ifx QS/regtest-moments-kpoints