comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add Andreussi's Solvent Aware SCCS method (#5495)

Browse source
This commit is contained in:
Ziwei Chai 2026-07-01 21:18:53 +08:00 • committed by GitHub
parent 9fc3e74bfe
commit 2c8265a529
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
10 changed files with 651 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

16
src/common/bibliography.F
View file

@ -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

10
src/cp_control_types.F
View file

@ -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
! **************************************************************************************************

59
src/cp_control_utils.F
View file

@ -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

3
src/input_constants.F
View file

@ -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, &

115
src/input_cp2k_dft.F
View file

@ -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
! **************************************************************************************************

52
src/pw/pw_poisson_methods.F
View file

@ -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

399
src/qs_sccs.F
View file

@ -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 ...

74
tests/QS/regtest-sasccs/H2_sasccs.inp Normal file
View file

@ -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

2
tests/QS/regtest-sasccs/TEST_FILES.toml Normal file
View file

@ -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}]

1
tests/TEST_DIRS
View file

@ -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