diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f7db346ef2..4583fdfe19 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -820,6 +820,7 @@ list( task_list_types.F tblite_interface.F tblite_ks_matrix.F + tblite_scc_mixer.F tblite_types.F tip_scan_methods.F tip_scan_types.F diff --git a/src/common/bibliography.F b/src/common/bibliography.F index 949cc94460..b8e4f5811c 100644 --- a/src/common/bibliography.F +++ b/src/common/bibliography.F @@ -77,7 +77,7 @@ MODULE bibliography Stoychev2016, Futera2017, Bailey2006, Papior2017, Lehtola2018, & Brieuc2016, Barca2018, Scheiber2018, Huang2011, Heaton_Burgess2007, & Schuett2018, Holmberg2018, Togo2018, Staub2019, Grimme2013, Grimme2016, & - Grimme2017, Asgeirsson2017, Bannwarth2019, Kondov2007, Clabaut2020, & + Grimme2017, Asgeirsson2017, Bannwarth2019, Katbashev2025, Kondov2007, Clabaut2020, & Clabaut2021, Ren2011, Ren2013, Cohen2000, Rogers2002, Filippetti2000, & Limpanuparb2011, Martin2003, Yin2017, Goerigk2017, & Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, Wilhelm2021, Lass2018, & @@ -785,6 +785,14 @@ CONTAINS source="J. Chem. Theory Comput.", volume="15", pages="1652-1671", & year=2019, doi="10.1021/acs.jctc.8b01176") + CALL add_reference(key=Katbashev2025, & + authors=s2a("A. Katbashev", "M. Stahn", "T. Rose", "V. Alizadeh", & + "M. Friede", "C. Plett", "P. Steinbach", "S. Ehlert"), & + title="Overview on Building Blocks and Applications of Efficient and Robust "// & + "Extended Tight Binding", & + source="J. Phys. Chem. A", volume="129", pages="2667-2682", & + year=2025, doi="10.1021/acs.jpca.4c08263") + CALL add_reference(key=Branduardi2007, & authors=s2a("D. Branduardi", "F. L. Gervasio", "M. Parrinello"), & title="From A to B in free energy space", & diff --git a/src/cp_control_types.F b/src/cp_control_types.F index ae2cee758c..26553d13fa 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -13,11 +13,13 @@ MODULE cp_control_types USE cp_fm_types, ONLY: cp_fm_release,& cp_fm_type USE eeq_input, ONLY: eeq_solver_type - USE input_constants, ONLY: do_full_density,& - rtp_bse_ham_G0W0,& - rtp_method_tddft,& - tblite_scc_mixer_auto,& - use_mom_ref_coac + USE input_constants, ONLY: & + do_full_density, rtp_bse_ham_G0W0, rtp_method_tddft, tblite_cli_born_kernel_auto, & + tblite_cli_solution_state_gsolv, tblite_cli_solvation_none, tblite_guess_ceh, & + tblite_guess_sad, tblite_mixer_damping_default, tblite_mixer_iterations_default, & + tblite_mixer_max_weight_default, tblite_mixer_min_weight_default, & + tblite_mixer_omega0_default, tblite_mixer_weight_factor_default, tblite_scc_mixer_auto, & + tblite_solver_gvd, use_mom_ref_coac USE kinds, ONLY: default_path_length,& default_string_length,& dp @@ -152,10 +154,17 @@ MODULE cp_control_types LOGICAL :: dispersion = .FALSE. INTEGER :: tblite_scc_mixer = tblite_scc_mixer_auto INTEGER :: dispersion_type = 0 + INTEGER :: tblite_mixer_iterations = tblite_mixer_iterations_default + INTEGER :: tblite_mixer_memory = tblite_mixer_iterations_default + INTEGER :: tblite_mixer_solver = tblite_solver_gvd LOGICAL :: dftb3_diagonal = .FALSE. LOGICAL :: hb_sr_damp = .FALSE. REAL(KIND=dp) :: hb_sr_para = 0.0_dp - REAL(KIND=dp) :: tblite_mixer_damping = 0.4_dp + REAL(KIND=dp) :: tblite_mixer_damping = tblite_mixer_damping_default + REAL(KIND=dp) :: tblite_mixer_max_weight = tblite_mixer_max_weight_default + REAL(KIND=dp) :: tblite_mixer_min_weight = tblite_mixer_min_weight_default + REAL(KIND=dp) :: tblite_mixer_omega0 = tblite_mixer_omega0_default + REAL(KIND=dp) :: tblite_mixer_weight_factor = tblite_mixer_weight_factor_default REAL(KIND=dp) :: eps_disp = 0.0_dp REAL(KIND=dp) :: epscn = 0.0_dp REAL(KIND=dp) :: exp_pre = 0.0_dp @@ -175,6 +184,42 @@ MODULE cp_control_types ! ************************************************************************************************** ! \brief Control parameters for xTB calculations ! ************************************************************************************************** + TYPE xtb_reference_cli_guess_type + LOGICAL :: enabled = .FALSE. + LOGICAL :: efield_active = .FALSE. + LOGICAL :: grad = .FALSE. + INTEGER :: method = tblite_guess_ceh + INTEGER :: solver = tblite_solver_gvd + REAL(KIND=dp), DIMENSION(3) :: efield = 0.0_dp + REAL(KIND=dp) :: electronic_temperature_guess = 0.0_dp + CHARACTER(LEN=default_path_length) :: input_file = "" + CHARACTER(LEN=default_path_length) :: json_file = "" + CHARACTER(LEN=default_string_length) :: input_format = "gen" + END TYPE xtb_reference_cli_guess_type + + TYPE xtb_reference_cli_param_type + LOGICAL :: enabled = .FALSE. + LOGICAL :: method_explicit = .FALSE. + INTEGER :: method = 0 + CHARACTER(LEN=default_path_length) :: input_file = "" + CHARACTER(LEN=default_path_length) :: output_file = "" + END TYPE xtb_reference_cli_param_type + + TYPE xtb_reference_cli_fit_type + LOGICAL :: enabled = .FALSE. + LOGICAL :: dry_run = .FALSE. + CHARACTER(LEN=default_path_length) :: copy_file = "" + CHARACTER(LEN=default_path_length) :: input_file = "" + CHARACTER(LEN=default_path_length) :: param_file = "" + END TYPE xtb_reference_cli_fit_type + + TYPE xtb_reference_cli_tagdiff_type + LOGICAL :: enabled = .FALSE. + LOGICAL :: fit = .FALSE. + CHARACTER(LEN=default_path_length) :: actual_file = "" + CHARACTER(LEN=default_path_length) :: reference_file = "" + END TYPE xtb_reference_cli_tagdiff_type + TYPE xtb_reference_cli_type LOGICAL :: enabled = .FALSE. LOGICAL :: check_energy = .TRUE. @@ -182,12 +227,29 @@ MODULE cp_control_types LOGICAL :: check_virial = .FALSE. LOGICAL :: keep_files = .FALSE. LOGICAL :: stop_on_error = .FALSE. - INTEGER :: iterations = 250 - REAL(KIND=dp) :: accuracy = 1.0E-8_dp + LOGICAL :: efield_active = .FALSE. + LOGICAL :: solvation_active = .FALSE. + INTEGER :: guess = tblite_guess_sad + INTEGER :: solvation_born_kernel = tblite_cli_born_kernel_auto + INTEGER :: solvation_model = tblite_cli_solvation_none + INTEGER :: solvation_state = tblite_cli_solution_state_gsolv + REAL(KIND=dp), DIMENSION(3) :: efield = 0.0_dp + REAL(KIND=dp) :: electronic_temperature_guess = 0.0_dp REAL(KIND=dp) :: error_limit = 1.0E-8_dp CHARACTER(LEN=default_path_length) :: program_name = "tblite" + CHARACTER(LEN=default_path_length) :: grad_file = "" + CHARACTER(LEN=default_path_length) :: json_file = "" + CHARACTER(LEN=default_path_length) :: post_processing = "" + CHARACTER(LEN=default_path_length) :: post_processing_output_file = "" + CHARACTER(LEN=default_path_length) :: restart_file = "" + CHARACTER(LEN=default_path_length) :: solvation_solvent = "" CHARACTER(LEN=default_path_length) :: work_directory = "." CHARACTER(LEN=default_path_length) :: prefix = "tblite-reference" + CHARACTER(LEN=default_string_length) :: input_format = "gen" + TYPE(xtb_reference_cli_guess_type) :: guess_cli + TYPE(xtb_reference_cli_param_type) :: param_cli + TYPE(xtb_reference_cli_fit_type) :: fit_cli + TYPE(xtb_reference_cli_tagdiff_type) :: tagdiff_cli END TYPE xtb_reference_cli_type TYPE xtb_control_type @@ -201,8 +263,16 @@ MODULE cp_control_types INTEGER :: h_sto_ng = 0 INTEGER :: tblite_method = 0 INTEGER :: tblite_scc_mixer = tblite_scc_mixer_auto + INTEGER :: tblite_mixer_iterations = tblite_mixer_iterations_default + INTEGER :: tblite_mixer_memory = tblite_mixer_iterations_default + INTEGER :: tblite_mixer_solver = tblite_solver_gvd REAL(KIND=dp) :: tblite_accuracy = 1.0_dp - REAL(KIND=dp) :: tblite_mixer_damping = 0.4_dp + REAL(KIND=dp) :: tblite_mixer_damping = tblite_mixer_damping_default + REAL(KIND=dp) :: tblite_mixer_max_weight = tblite_mixer_max_weight_default + REAL(KIND=dp) :: tblite_mixer_min_weight = tblite_mixer_min_weight_default + REAL(KIND=dp) :: tblite_mixer_omega0 = tblite_mixer_omega0_default + REAL(KIND=dp) :: tblite_mixer_weight_factor = tblite_mixer_weight_factor_default + CHARACTER(LEN=default_path_length) :: tblite_param_file = "" ! INTEGER :: vdw_type = -1 CHARACTER(LEN=default_path_length) :: parameter_file_path = "" diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 22db8af4ca..8afee457cd 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, Asgeirsson2017, Bannwarth2019, Caldeweyher2017, Caldeweyher2020, Dewar1977, & - Dewar1985, Elstner1998, Fattebert2002, Grimme2017, Hu2007, Krack2000, Lippert1997, & - Lippert1999, Porezag1995, Pracht2019, Repasky2002, Rocha2006, Schenter2008, Seifert1996, & - Souza2002, Stengel2009, Stewart1989, Stewart2007, Thiel1992, Umari2002, VanVoorhis2015, & - VandeVondele2005a, VandeVondele2005b, Yin2017, Zhechkov2005, cite_reference + 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 USE cp_control_types, ONLY: & admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, & dft_control_create, dft_control_type, efield_type, expot_control_create, & @@ -57,9 +57,13 @@ MODULE cp_control_utils 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_scc_mixer_auto, tblite_scc_mixer_cp2k, & - tblite_scc_mixer_none, tblite_scc_mixer_tblite, tddfpt_dipole_length, tddfpt_kernel_stda, & - use_mom_ref_user, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none + 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, & + tblite_scc_mixer_cp2k, tblite_scc_mixer_none, tblite_scc_mixer_tblite, tblite_solver_gvd, & + tblite_solver_gvr, tddfpt_dipole_length, tddfpt_kernel_stda, use_mom_ref_user, & + xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none USE input_cp2k_check, ONLY: xc_functionals_expand USE input_cp2k_dft, ONLY: create_dft_section USE input_enumeration_types, ONLY: enum_i2c,& @@ -877,16 +881,17 @@ CONTAINS REAL(dp) :: tmp, tmpsqrt, value REAL(dp), POINTER :: scal(:) TYPE(section_vals_type), POINTER :: cdft_control_section, ddapc_restraint_section, & - dftb_parameter, dftb_section, eeq_section, genpot_section, lri_optbas_section, & - mull_section, nonbonded_section, s2_restraint_section, se_section, xtb_parameter, & - xtb_section, xtb_tblite, xtb_tblite_ref_cli + dftb_parameter, dftb_section, dftb_tblite_mixer, eeq_section, genpot_section, & + lri_optbas_section, mull_section, nonbonded_section, s2_restraint_section, se_section, & + xtb_parameter, xtb_section, xtb_tblite, xtb_tblite_mixer, xtb_tblite_ref_cli CALL timeset(routineN, handle) was_present = .FALSE. NULLIFY (mull_section, ddapc_restraint_section, s2_restraint_section, & se_section, dftb_section, xtb_section, dftb_parameter, xtb_parameter, lri_optbas_section, & - cdft_control_section, genpot_section, eeq_section, xtb_tblite_ref_cli) + cdft_control_section, genpot_section, eeq_section, dftb_tblite_mixer, & + xtb_tblite_mixer, xtb_tblite_ref_cli) mull_section => section_vals_get_subs_vals(qs_section, "MULLIKEN_RESTRAINT") ddapc_restraint_section => section_vals_get_subs_vals(qs_section, "DDAPC_RESTRAINT") @@ -895,12 +900,14 @@ CONTAINS dftb_section => section_vals_get_subs_vals(qs_section, "DFTB") xtb_section => section_vals_get_subs_vals(qs_section, "xTB") dftb_parameter => section_vals_get_subs_vals(dftb_section, "PARAMETER") + dftb_tblite_mixer => section_vals_get_subs_vals(dftb_section, "TBLITE_MIXER") xtb_parameter => section_vals_get_subs_vals(xtb_section, "PARAMETER") eeq_section => section_vals_get_subs_vals(xtb_section, "EEQ") lri_optbas_section => section_vals_get_subs_vals(qs_section, "OPTIMIZE_LRI_BASIS") cdft_control_section => section_vals_get_subs_vals(qs_section, "CDFT") nonbonded_section => section_vals_get_subs_vals(xtb_section, "NONBONDED") genpot_section => section_vals_get_subs_vals(nonbonded_section, "GENPOT") + xtb_tblite_mixer => section_vals_get_subs_vals(xtb_section, "TBLITE_MIXER") xtb_tblite => section_vals_get_subs_vals(xtb_section, "TBLITE") xtb_tblite_ref_cli => section_vals_get_subs_vals(xtb_tblite, "REFERENCE_CLI") @@ -1296,8 +1303,18 @@ CONTAINS l_val=qs_control%dftb_control%hb_sr_damp) CALL section_vals_val_get(dftb_section, "SCC_MIXER", & i_val=qs_control%dftb_control%tblite_scc_mixer) - CALL section_vals_val_get(dftb_section, "TBLITE_MIXER_DAMPING", & - r_val=qs_control%dftb_control%tblite_mixer_damping) + CALL read_tblite_mixer_section(dftb_tblite_mixer, & + qs_control%dftb_control%tblite_mixer_iterations, & + qs_control%dftb_control%tblite_mixer_memory, & + qs_control%dftb_control%tblite_mixer_solver, & + qs_control%dftb_control%tblite_mixer_damping, & + qs_control%dftb_control%tblite_mixer_omega0, & + qs_control%dftb_control%tblite_mixer_min_weight, & + qs_control%dftb_control%tblite_mixer_max_weight, & + qs_control%dftb_control%tblite_mixer_weight_factor, & + "DFTB/TBLITE_MIXER") + IF (qs_control%dftb_control%tblite_mixer_damping <= 0.0_dp) & + CPABORT("DFTB/TBLITE_MIXER/DAMPING must be positive") CALL section_vals_val_get(dftb_section, "EPS_DISP", & r_val=qs_control%dftb_control%eps_disp) CALL section_vals_val_get(dftb_section, "DO_EWALD", explicit=explicit) @@ -1366,8 +1383,18 @@ CONTAINS END IF CALL section_vals_val_get(xtb_section, "SCC_MIXER", & i_val=qs_control%xtb_control%tblite_scc_mixer) - CALL section_vals_val_get(xtb_section, "TBLITE_MIXER_DAMPING", & - r_val=qs_control%xtb_control%tblite_mixer_damping) + CALL read_tblite_mixer_section(xtb_tblite_mixer, & + qs_control%xtb_control%tblite_mixer_iterations, & + qs_control%xtb_control%tblite_mixer_memory, & + qs_control%xtb_control%tblite_mixer_solver, & + qs_control%xtb_control%tblite_mixer_damping, & + qs_control%xtb_control%tblite_mixer_omega0, & + qs_control%xtb_control%tblite_mixer_min_weight, & + qs_control%xtb_control%tblite_mixer_max_weight, & + qs_control%xtb_control%tblite_mixer_weight_factor, & + "XTB/TBLITE_MIXER") + IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) & + CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive") CALL section_vals_val_get(xtb_section, "DO_EWALD", explicit=explicit) IF (explicit) THEN CALL section_vals_val_get(xtb_section, "DO_EWALD", & @@ -1683,14 +1710,14 @@ CONTAINS IF (qs_control%xtb_control%do_tblite) THEN CALL section_vals_val_get(xtb_tblite, "METHOD", & i_val=qs_control%xtb_control%tblite_method) + CALL section_vals_val_get(xtb_tblite, "PARAM", & + c_val=qs_control%xtb_control%tblite_param_file) CALL section_vals_val_get(xtb_tblite, "ACCURACY", & r_val=qs_control%xtb_control%tblite_accuracy) IF (qs_control%xtb_control%tblite_accuracy <= 0.0_dp) & CPABORT("XTB/TBLITE/ACCURACY must be positive") - CALL section_vals_val_get(xtb_tblite, "SCC_MIXER", & - i_val=qs_control%xtb_control%tblite_scc_mixer) - CALL section_vals_val_get(xtb_tblite, "TBLITE_MIXER_DAMPING", & - r_val=qs_control%xtb_control%tblite_mixer_damping) + IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) & + CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive") CALL section_vals_val_get(xtb_tblite, "REFERENCE_CLI", l_val=tblite_reference_cli) CALL section_vals_get(xtb_tblite_ref_cli, explicit=tblite_reference_cli_section) IF (tblite_reference_cli .AND. (.NOT. tblite_reference_cli_section)) & @@ -1699,11 +1726,7 @@ CONTAINS CALL read_xtb_reference_cli_section(xtb_tblite_ref_cli, qs_control%xtb_control%reference_cli) qs_control%xtb_control%reference_cli%enabled = .TRUE. END IF - CALL cite_reference(Caldeweyher2017) - CALL cite_reference(Caldeweyher2020) - CALL cite_reference(Asgeirsson2017) - CALL cite_reference(Grimme2017) - CALL cite_reference(Bannwarth2019) + CALL cite_reference(Katbashev2025) ! tblite handles periodic long-range terms internally from the CP2K cell periodicity. ! Keep xtb_control%do_ewald as read above from XTB/DO_EWALD or SUBSYS/CELL/PERIODIC, ! matching the DFTB and CP2K-internal xTB setup. @@ -1712,6 +1735,62 @@ CONTAINS CALL timestop(handle) END SUBROUTINE read_qs_section +! ************************************************************************************************** +!> \brief Read a TBLITE_MIXER section. +!> \param mixer_section input section +!> \param iterations tblite SCC iteration limit +!> \param memory Broyden history length +!> \param solver native tblite electronic solver id +!> \param damping mixer damping parameter +!> \param omega0 Broyden regularization weight +!> \param min_weight minimum dynamic Broyden weight +!> \param max_weight maximum dynamic Broyden weight +!> \param weight_factor residual-to-weight scaling factor +!> \param section_name diagnostic section name +! ************************************************************************************************** + SUBROUTINE read_tblite_mixer_section(mixer_section, iterations, memory, solver, damping, omega0, min_weight, & + max_weight, weight_factor, section_name) + TYPE(section_vals_type), POINTER :: mixer_section + INTEGER, INTENT(INOUT) :: iterations, memory, solver + REAL(KIND=dp), INTENT(INOUT) :: damping, omega0, min_weight, max_weight, & + weight_factor + CHARACTER(LEN=*), INTENT(IN) :: section_name + + LOGICAL :: explicit, memory_explicit + + CALL section_vals_get(mixer_section, explicit=explicit) + IF (.NOT. explicit) RETURN + + CALL section_vals_val_get(mixer_section, "ITERATIONS", i_val=iterations) + CALL section_vals_val_get(mixer_section, "MEMORY", explicit=memory_explicit) + IF (memory_explicit) CALL section_vals_val_get(mixer_section, "MEMORY", i_val=memory) + IF (.NOT. memory_explicit .OR. memory == tblite_mixer_memory_inherit) THEN + memory = iterations + END IF + CALL section_vals_val_get(mixer_section, "SOLVER", i_val=solver) + CALL section_vals_val_get(mixer_section, "DAMPING", r_val=damping) + CALL section_vals_val_get(mixer_section, "OMEGA0", r_val=omega0) + CALL section_vals_val_get(mixer_section, "MIN_WEIGHT", r_val=min_weight) + CALL section_vals_val_get(mixer_section, "MAX_WEIGHT", r_val=max_weight) + CALL section_vals_val_get(mixer_section, "WEIGHT_FACTOR", r_val=weight_factor) + + IF (iterations < 1) CPABORT(TRIM(section_name)//"/ITERATIONS must be positive") + IF (memory < 1) CPABORT(TRIM(section_name)//"/MEMORY must be positive") + SELECT CASE (solver) + CASE (tblite_solver_gvd, tblite_solver_gvr) + CASE DEFAULT + CPABORT(TRIM(section_name)//"/SOLVER must be GVD or GVR") + END SELECT + IF (damping <= 0.0_dp) CPABORT(TRIM(section_name)//"/DAMPING must be positive") + IF (omega0 <= 0.0_dp) CPABORT(TRIM(section_name)//"/OMEGA0 must be positive") + IF (min_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MIN_WEIGHT must be positive") + IF (max_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MAX_WEIGHT must be positive") + IF (max_weight < min_weight) & + CPABORT(TRIM(section_name)//"/MAX_WEIGHT must not be smaller than MIN_WEIGHT") + IF (weight_factor <= 0.0_dp) CPABORT(TRIM(section_name)//"/WEIGHT_FACTOR must be positive") + + END SUBROUTINE read_tblite_mixer_section + ! ************************************************************************************************** !> \brief Read native tblite CLI reference options. !> \param ref_cli_section input section @@ -1721,18 +1800,113 @@ CONTAINS TYPE(section_vals_type), POINTER :: ref_cli_section TYPE(xtb_reference_cli_type), INTENT(INOUT) :: ref_cli + REAL(KIND=dp), DIMENSION(:), POINTER :: efield + TYPE(section_vals_type), POINTER :: fit_section, guess_section, & + param_section, solvation_section, & + tagdiff_section + CALL section_vals_val_get(ref_cli_section, "_SECTION_PARAMETERS_", l_val=ref_cli%enabled) CALL section_vals_val_get(ref_cli_section, "PROGRAM_NAME", c_val=ref_cli%program_name) + CALL section_vals_val_get(ref_cli_section, "GUESS", i_val=ref_cli%guess) CALL section_vals_val_get(ref_cli_section, "WORK_DIRECTORY", c_val=ref_cli%work_directory) CALL section_vals_val_get(ref_cli_section, "PREFIX", c_val=ref_cli%prefix) + CALL section_vals_val_get(ref_cli_section, "INPUT_FORMAT", c_val=ref_cli%input_format) + CALL section_vals_val_get(ref_cli_section, "RESTART", c_val=ref_cli%restart_file) + CALL section_vals_val_get(ref_cli_section, "GRAD", c_val=ref_cli%grad_file) + CALL section_vals_val_get(ref_cli_section, "JSON", c_val=ref_cli%json_file) + CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING", c_val=ref_cli%post_processing) + CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING_OUTPUT", & + c_val=ref_cli%post_processing_output_file) + CALL section_vals_val_get(ref_cli_section, "EFIELD", explicit=ref_cli%efield_active) + IF (ref_cli%efield_active) THEN + NULLIFY (efield) + CALL section_vals_val_get(ref_cli_section, "EFIELD", r_vals=efield) + ref_cli%efield = efield(1:3) + END IF + solvation_section => section_vals_get_subs_vals(ref_cli_section, "IMPLICIT_SOLVATION") + CALL section_vals_get(solvation_section, explicit=ref_cli%solvation_active) + IF (ref_cli%solvation_active) THEN + CALL section_vals_val_get(solvation_section, "MODEL", i_val=ref_cli%solvation_model) + CALL section_vals_val_get(solvation_section, "SOLVENT", c_val=ref_cli%solvation_solvent) + CALL section_vals_val_get(solvation_section, "BORN_KERNEL", i_val=ref_cli%solvation_born_kernel) + CALL section_vals_val_get(solvation_section, "SOLUTION_STATE", i_val=ref_cli%solvation_state) + IF (LEN_TRIM(ref_cli%solvation_solvent) == 0) & + CPABORT("REFERENCE_CLI implicit solvation needs SOLVENT") + IF (ref_cli%solvation_model == tblite_cli_solvation_cpcm .AND. & + ref_cli%solvation_born_kernel /= tblite_cli_born_kernel_auto) & + CPABORT("BORN_KERNEL is invalid with MODEL CPCM") + IF (ref_cli%solvation_state /= tblite_cli_solution_state_gsolv) THEN + SELECT CASE (ref_cli%solvation_model) + CASE (tblite_cli_solvation_alpb, tblite_cli_solvation_gbsa) + ! Native tblite supports solution-state shifts for parametrized named-solvent ALPB/GBSA. + CASE (tblite_cli_solvation_gbe, tblite_cli_solvation_gb, tblite_cli_solvation_cpcm) + CPABORT("SOLUTION_STATE is valid only for ALPB/GBSA") + END SELECT + END IF + END IF + CALL section_vals_val_get(ref_cli_section, "ELECTRONIC_TEMPERATURE_GUESS", & + r_val=ref_cli%electronic_temperature_guess) + IF (ref_cli%electronic_temperature_guess < 0.0_dp) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative") + IF (ref_cli%electronic_temperature_guess > 0.0_dp .AND. ref_cli%guess /= tblite_guess_ceh) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS requires GUESS CEH") + guess_section => section_vals_get_subs_vals(ref_cli_section, "GUESS_CLI") + CALL section_vals_get(guess_section, explicit=ref_cli%guess_cli%enabled) + IF (ref_cli%guess_cli%enabled) THEN + CALL section_vals_val_get(guess_section, "METHOD", i_val=ref_cli%guess_cli%method) + CALL section_vals_val_get(guess_section, "ELECTRONIC_TEMPERATURE_GUESS", & + r_val=ref_cli%guess_cli%electronic_temperature_guess) + IF (ref_cli%guess_cli%electronic_temperature_guess < 0.0_dp) & + CPABORT("REFERENCE_CLI/GUESS_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative") + CALL section_vals_val_get(guess_section, "SOLVER", i_val=ref_cli%guess_cli%solver) + CALL section_vals_val_get(guess_section, "EFIELD", explicit=ref_cli%guess_cli%efield_active) + IF (ref_cli%guess_cli%efield_active) THEN + NULLIFY (efield) + CALL section_vals_val_get(guess_section, "EFIELD", r_vals=efield) + ref_cli%guess_cli%efield = efield(1:3) + END IF + CALL section_vals_val_get(guess_section, "GRAD", l_val=ref_cli%guess_cli%grad) + CALL section_vals_val_get(guess_section, "JSON", c_val=ref_cli%guess_cli%json_file) + CALL section_vals_val_get(guess_section, "INPUT_FORMAT", c_val=ref_cli%guess_cli%input_format) + CALL section_vals_val_get(guess_section, "INPUT_FILE", c_val=ref_cli%guess_cli%input_file) + END IF + param_section => section_vals_get_subs_vals(ref_cli_section, "PARAM_CLI") + CALL section_vals_get(param_section, explicit=ref_cli%param_cli%enabled) + IF (ref_cli%param_cli%enabled) THEN + CALL section_vals_val_get(param_section, "METHOD", explicit=ref_cli%param_cli%method_explicit, & + i_val=ref_cli%param_cli%method) + CALL section_vals_val_get(param_section, "OUTPUT", c_val=ref_cli%param_cli%output_file) + CALL section_vals_val_get(param_section, "INPUT_FILE", c_val=ref_cli%param_cli%input_file) + END IF + fit_section => section_vals_get_subs_vals(ref_cli_section, "FIT_CLI") + CALL section_vals_get(fit_section, explicit=ref_cli%fit_cli%enabled) + IF (ref_cli%fit_cli%enabled) THEN + CALL section_vals_val_get(fit_section, "PARAM_FILE", c_val=ref_cli%fit_cli%param_file) + CALL section_vals_val_get(fit_section, "INPUT_FILE", c_val=ref_cli%fit_cli%input_file) + CALL section_vals_val_get(fit_section, "DRY_RUN", l_val=ref_cli%fit_cli%dry_run) + CALL section_vals_val_get(fit_section, "COPY", c_val=ref_cli%fit_cli%copy_file) + IF (LEN_TRIM(ref_cli%fit_cli%param_file) == 0) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs PARAM_FILE") + IF (LEN_TRIM(ref_cli%fit_cli%input_file) == 0) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs INPUT_FILE") + END IF + tagdiff_section => section_vals_get_subs_vals(ref_cli_section, "TAGDIFF_CLI") + CALL section_vals_get(tagdiff_section, explicit=ref_cli%tagdiff_cli%enabled) + IF (ref_cli%tagdiff_cli%enabled) THEN + CALL section_vals_val_get(tagdiff_section, "ACTUAL", c_val=ref_cli%tagdiff_cli%actual_file) + CALL section_vals_val_get(tagdiff_section, "REFERENCE", c_val=ref_cli%tagdiff_cli%reference_file) + CALL section_vals_val_get(tagdiff_section, "FIT", l_val=ref_cli%tagdiff_cli%fit) + IF (LEN_TRIM(ref_cli%tagdiff_cli%actual_file) == 0) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs ACTUAL") + IF (LEN_TRIM(ref_cli%tagdiff_cli%reference_file) == 0) & + CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs REFERENCE") + END IF CALL section_vals_val_get(ref_cli_section, "KEEP_FILES", l_val=ref_cli%keep_files) CALL section_vals_val_get(ref_cli_section, "ERROR_LIMIT", r_val=ref_cli%error_limit) CALL section_vals_val_get(ref_cli_section, "STOP_ON_ERROR", l_val=ref_cli%stop_on_error) CALL section_vals_val_get(ref_cli_section, "CHECK_ENERGY", l_val=ref_cli%check_energy) CALL section_vals_val_get(ref_cli_section, "CHECK_FORCES", l_val=ref_cli%check_forces) CALL section_vals_val_get(ref_cli_section, "CHECK_VIRIAL", l_val=ref_cli%check_virial) - CALL section_vals_val_get(ref_cli_section, "ACCURACY", r_val=ref_cli%accuracy) - CALL section_vals_val_get(ref_cli_section, "ITERATIONS", i_val=ref_cli%iterations) END SUBROUTINE read_xtb_reference_cli_section @@ -2362,7 +2536,7 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'write_xtb_control' - CHARACTER(LEN=16) :: scc_mixer_name + CHARACTER(LEN=16) :: scc_mixer_name, solver_name INTEGER :: handle, output_unit TYPE(cp_logger_type), POINTER :: logger @@ -2432,12 +2606,38 @@ CONTAINS CASE DEFAULT CPABORT("Unknown tblite SCC mixer") END SELECT + SELECT CASE (xtb_control%tblite_mixer_solver) + CASE (tblite_solver_gvd) + solver_name = "GVD" + CASE (tblite_solver_gvr) + solver_name = "GVR" + CASE DEFAULT + CPABORT("Unknown tblite SCC mixer solver") + END SELECT WRITE (UNIT=output_unit, FMT="(T2,A,T72,A)") & "xTB| SCC mixer:", TRIM(scc_mixer_name) WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") & "xTB| tblite accuracy:", xtb_control%tblite_accuracy + IF (LEN_TRIM(xtb_control%tblite_param_file) > 0) THEN + WRITE (UNIT=output_unit, FMT="(T2,A,T33,A)") & + "xTB| tblite parameter file:", TRIM(xtb_control%tblite_param_file) + END IF WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & "xTB| tblite SCC mixer damping:", xtb_control%tblite_mixer_damping + WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & + "xTB| tblite SCC mixer iterations:", xtb_control%tblite_mixer_iterations + WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") & + "xTB| tblite SCC mixer memory:", xtb_control%tblite_mixer_memory + WRITE (UNIT=output_unit, FMT="(T2,A,T72,A)") & + "xTB| tblite SCC mixer solver:", TRIM(solver_name) + WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") & + "xTB| tblite SCC mixer omega0:", xtb_control%tblite_mixer_omega0 + WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") & + "xTB| tblite SCC mixer min weight:", xtb_control%tblite_mixer_min_weight + WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") & + "xTB| tblite SCC mixer max weight:", xtb_control%tblite_mixer_max_weight + WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") & + "xTB| tblite SCC mixer weight factor:", xtb_control%tblite_mixer_weight_factor WRITE (UNIT=output_unit, FMT="(/)") END IF diff --git a/src/input_constants.F b/src/input_constants.F index 3a683f126a..5eee5476e0 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1361,6 +1361,40 @@ MODULE input_constants tblite_scc_mixer_cp2k = 2, & tblite_scc_mixer_none = 3 + ! tblite SCC mixer defaults matching the tblite library + INTEGER, PARAMETER, PUBLIC :: tblite_mixer_iterations_default = 250, & + tblite_mixer_memory_inherit = 0 + REAL(KIND=dp), PARAMETER, PUBLIC :: tblite_mixer_damping_default = 0.4_dp, & + tblite_mixer_max_weight_default = 100000.0_dp, & + tblite_mixer_min_weight_default = 1.0_dp, & + tblite_mixer_omega0_default = 0.01_dp, & + tblite_mixer_weight_factor_default = 0.01_dp + + ! tblite reference CLI initial guess selection + INTEGER, PARAMETER, PUBLIC :: tblite_guess_sad = 0, & + tblite_guess_eeq = 1, & + tblite_guess_ceh = 2 + + ! tblite native electronic solver selection + INTEGER, PARAMETER, PUBLIC :: tblite_solver_gvd = 0, & + tblite_solver_gvr = 1 + + ! tblite reference CLI implicit solvation selection + INTEGER, PARAMETER, PUBLIC :: tblite_cli_solvation_none = 0, & + tblite_cli_solvation_alpb = 1, & + tblite_cli_solvation_gbsa = 2, & + tblite_cli_solvation_gbe = 3, & + tblite_cli_solvation_gb = 4, & + tblite_cli_solvation_cpcm = 5 + + INTEGER, PARAMETER, PUBLIC :: tblite_cli_born_kernel_auto = 0, & + tblite_cli_born_kernel_p16 = 1, & + tblite_cli_born_kernel_still = 2 + + INTEGER, PARAMETER, PUBLIC :: tblite_cli_solution_state_gsolv = 0, & + tblite_cli_solution_state_bar1mol = 1, & + tblite_cli_solution_state_reference = 2 + ! SMEAGOL interface INTEGER, PARAMETER, PUBLIC :: smeagol_bulklead_leftright = 1, & smeagol_bulklead_left = 2, & diff --git a/src/input_cp2k_tb.F b/src/input_cp2k_tb.F index 54ba18b40f..7c62a26197 100644 --- a/src/input_cp2k_tb.F +++ b/src/input_cp2k_tb.F @@ -12,14 +12,27 @@ !> \author fawzi ! ************************************************************************************************** MODULE input_cp2k_tb - USE bibliography, ONLY: & - Asgeirsson2017, Bannwarth2019, Caldeweyher2017, Caldeweyher2020, Elstner1998, Grimme2017, & - Hu2007, Porezag1995, Seifert1996, Zhechkov2005 + USE bibliography, ONLY: Elstner1998,& + Grimme2017,& + Hu2007,& + Katbashev2025,& + Porezag1995,& + Seifert1996,& + Zhechkov2005 + USE cp_units, ONLY: cp_unit_to_cp2k USE eeq_input, ONLY: create_eeq_control_section USE input_constants, ONLY: & dispersion_d2, dispersion_d3, dispersion_d3bj, dispersion_uff, gfn1xtb, gfn2xtb, & - gfn_tblite, ipea1xtb, slater, tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, & - tblite_scc_mixer_none, tblite_scc_mixer_tblite + gfn_tblite, ipea1xtb, slater, tblite_cli_born_kernel_auto, tblite_cli_born_kernel_p16, & + tblite_cli_born_kernel_still, tblite_cli_solution_state_bar1mol, & + tblite_cli_solution_state_gsolv, tblite_cli_solution_state_reference, & + tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, tblite_cli_solvation_gb, & + tblite_cli_solvation_gbe, tblite_cli_solvation_gbsa, tblite_guess_ceh, tblite_guess_eeq, & + tblite_guess_sad, tblite_mixer_damping_default, tblite_mixer_iterations_default, & + tblite_mixer_max_weight_default, tblite_mixer_memory_inherit, & + tblite_mixer_min_weight_default, tblite_mixer_omega0_default, & + tblite_mixer_weight_factor_default, tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, & + tblite_scc_mixer_none, tblite_scc_mixer_tblite, tblite_solver_gvd, tblite_solver_gvr USE input_cp2k_mm, ONLY: create_GENPOT_section USE input_keyword_types, ONLY: keyword_create,& keyword_release,& @@ -29,7 +42,8 @@ MODULE input_cp2k_tb section_create,& section_release,& section_type - USE input_val_types, ONLY: char_t + USE input_val_types, ONLY: char_t,& + real_t USE kinds, ONLY: dp USE string_utilities, ONLY: s2a #include "./base/base_uses.f90" @@ -64,6 +78,11 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + NULLIFY (subsection) + CALL create_tblite_mixer_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="self_consistent", & description="Use self-consistent method", & @@ -150,6 +169,10 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_tblite_mixer_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + CALL create_xtb_tblite_section(subsection) CALL section_add_subsection(section, subsection) CALL section_release(subsection) @@ -169,8 +192,9 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL add_tb_scc_mixer_keywords(section, "CP2K-internal xTB", & - "AUTO keeps the CP2K charge mixer for CP2K-internal xTB.") + CALL add_tb_scc_mixer_keywords(section, "xTB", & + "AUTO keeps the CP2K charge mixer for CP2K-internal xTB and uses "// & + "tblite's native SCC mixer for XTB/GFN_TYPE TBLITE.") CALL keyword_create(keyword, __LOCATION__, name="DO_EWALD", & description="Deprecated debug override for the xTB Ewald treatment. "// & @@ -579,50 +603,85 @@ CONTAINS description="Section used to specify options for an xTB computation using tblite. "// & "This section is required when XTB/GFN_TYPE is set to TBLITE. "// & "Periodicity follows SUBSYS/CELL/PERIODIC; the deprecated XTB/DO_EWALD keyword "// & - "is read consistently with CP2K-internal xTB for diagnostics/debugging.", & - n_keywords=1, n_subsections=1, repeats=.FALSE., citations=[Caldeweyher2017, Caldeweyher2020, & - Asgeirsson2017, Grimme2017, Bannwarth2019]) + "is read consistently with CP2K-internal xTB for diagnostics/debugging. "// & + "Initial densities and restarts are controlled by the regular CP2K SCF settings, "// & + "for example DFT/SCF/SCF_GUESS, not by a separate tblite guess in this section. "// & + "Spin-polarized CP2K/tblite calculations are selected with the regular DFT/UKS keyword "// & + "or its LSD/SPIN_POLARIZED aliases; DFT/MULTIPLICITY sets the number of unpaired "// & + "electrons. UKS uses tblite's spin-polarized Hamiltonian; AUTO selects the native tblite "// & + "SCC mixer for XTB/GFN_TYPE TBLITE, while SCC_MIXER CP2K mixes both total-charge and "// & + "magnetization SCC variables through CP2K's TB charge-mixing infrastructure. Native "// & + "tblite SCC convergence follows ACCURACY; SCC_MIXER CP2K follows the regular CP2K "// & + "SCF convergence settings such as DFT/SCF/EPS_SCF.", & + n_keywords=1, n_subsections=1, repeats=.FALSE., citations=[Katbashev2025]) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & - description="Marks the tblite options section as present. The tblite backend is selected "// & - "with XTB/GFN_TYPE TBLITE.", & + description="Marks the tblite options section as present. The actual switch to the tblite "// & + "backend is XTB/GFN_TYPE TBLITE; this section then collects the numerical model, SCC "// & + "mixing, and optional native-reference settings used by that backend.", & lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="METHOD", & - description="Selection of the method used in tblite.", & + description="Selects the xTB Hamiltonian model used by tblite. Without PARAM this chooses "// & + "one of tblite's built-in parametrizations, including the basis, reference occupations, "// & + "repulsion, dispersion, electrostatic, and coordination-number dependent terms. When PARAM "// & + "is given, METHOD remains the declared model family that the custom parametrization is "// & + "expected to implement and that CP2K uses for method-dependent setup and diagnostics.", & usage="METHOD (GFN1|GFN2|IPEA1)", & enum_c_vals=s2a("GFN1", "GFN2", "IPEA1"), & enum_i_vals=[gfn1xtb, gfn2xtb, ipea1xtb], & enum_desc=s2a("Uses the GFN1-XTB method by Grimme.", & "Uses the GFN2-XTB method by Grimme.", & - "Uses the IEPEA1 method by Grimme."), & + "Uses the IPEA1 method by Grimme."), & default_i_val=gfn2xtb) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL add_tb_scc_mixer_keywords(section, "CP2K/tblite xTB", & - "AUTO uses tblite's internal SCC mixer.") + CALL keyword_create(keyword, __LOCATION__, name="PARAM", & + description="Custom tblite parametrization file. The file is read by tblite and used to "// & + "construct the xTB calculator instead of the built-in parameter table for METHOD. It may "// & + "therefore change model constants such as Hamiltonian, repulsion, dispersion, electrostatic, "// & + "and coordination-number response parameters. It does not change the CP2K charge, spin, "// & + "cell, k-point, or SCF setup. The parametrization must be compatible with the selected "// & + "METHOD; CP2K does not reinterpret or validate the physical content of the file.", & + usage="PARAM filename", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="ACCURACY", & - description="Numerical accuracy factor for the CP2K/tblite backend, matching native "// & - "tblite's --acc option. The default matches native tblite. Smaller values tighten, "// & - "larger values loosen, the SCC convergence target and integral cutoff.", & + description="Numerical accuracy factor for the CP2K/tblite backend. This scales tblite's "// & + "internal accuracy-dependent thresholds. With the native tblite SCC mixer, this is the "// & + "user-visible convergence control and corresponds to native tblite's --acc setting. With "// & + "SCC_MIXER CP2K, SCF convergence instead follows DFT/SCF/EPS_SCF; ACCURACY still controls "// & + "the tblite numerical setup, such as the distance cutoff used to build the xTB basis and "// & + "neighbor data. "// & + "The default matches native tblite. Smaller values make the calculation tighter and usually "// & + "more expensive; larger values loosen the numerical setup and can affect energies, forces, "// & + "and stress. Native REFERENCE_CLI diagnostics inherit this value so the production and "// & + "external reference calculations use the same accuracy setting.", & usage="ACCURACY 1.0", default_r_val=1.0_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="REFERENCE_CLI", & - description="Enable the native tblite CLI reference check. If true, the "// & - "XTB/TBLITE/REFERENCE_CLI section must be present. The reference check is also "// & - "enabled whenever that section is present. This diagnostic path is skipped for "// & - "CP2K multi-k-point calculations, because the native tblite CLI does not reproduce "// & - "CP2K KPOINTS sampling.", & + description="Enable an optional diagnostic comparison against an external native tblite "// & + "command line calculation. CP2K writes the current geometry, charge, spin, method, optional "// & + "PARAM file, electronic temperature, and requested reference settings to temporary files, "// & + "runs the external program, and compares the resulting energy, gradient, and optionally "// & + "virial with the CP2K/tblite values. This check is meant for debugging backend consistency, "// & + "not for production sampling. The native tblite CLI exposes only an electronic temperature "// & + "for the occupation broadening; if CP2K uses another smearing method, the reference run "// & + "prints a warning and uses Fermi-Dirac smearing with the CP2K electronic temperature. "// & + "If true, the XTB/TBLITE/REFERENCE_CLI section must be present. "// & + "The check is also enabled whenever that section is present. It is skipped for CP2K "// & + "multi-k-point calculations, because the native tblite CLI does not reproduce CP2K KPOINTS "// & + "sampling.", & usage="REFERENCE_CLI", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -634,6 +693,121 @@ CONTAINS END SUBROUTINE create_xtb_tblite_section +! ************************************************************************************************** +!> \brief Creates the common &TBLITE_MIXER section for TB SCC-variable mixing. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_tblite_mixer_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="TBLITE_MIXER", & + description="Controls tblite's modified Broyden SCC-variable mixer when SCC_MIXER is set "// & + "to TBLITE, either explicitly or through AUTO where AUTO selects tblite. The mixer "// & + "variables are the self-consistent "// & + "xTB/DFTB populations and, for methods that require them, atomic multipoles. ITERATIONS "// & + "limits the number of tblite SCC-mixer updates, while MEMORY controls how many previous "// & + "Broyden update vectors are retained. These are independent quantities when MEMORY is set "// & + "explicitly; otherwise MEMORY follows ITERATIONS, matching the tblite library default. "// & + "DAMPING controls the size of the retained population/multipole update. OMEGA0, "// & + "MIN_WEIGHT, MAX_WEIGHT, and WEIGHT_FACTOR expose tblite's modified-Broyden weighting "// & + "constants. ITERATIONS is the tblite-side SCC-mixer update limit and defaults to the "// & + "native tblite value of 250. When SCC_MIXER is explicitly set to TBLITE, CP2K also "// & + "sets the regular DFT/SCF/MAX_SCF limit to ITERATIONS so that the tblite SCC cycle is "// & + "not truncated by CP2K's default SCF limit. With SCC_MIXER AUTO or CP2K, DFT/SCF/MAX_SCF "// & + "keeps its normal CP2K meaning. This section configures the tblite SCC mixer, not the "// & + "regular CP2K charge mixer. With SCC_MIXER CP2K, the CP2K-side SCC-variable mixer is "// & + "controlled by DFT/SCF/MIXING. Native tblite "// & + "REFERENCE_CLI diagnostics inherit this section. ITERATIONS and SOLVER are representable "// & + "by tblite run directly; CP2K reports inherited MEMORY, DAMPING, or Broyden-weight settings "// & + "that native tblite CLI cannot reproduce.", & + n_keywords=1, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="ITERATIONS", & + description="Maximum number of tblite SCC-mixer updates. This is the tblite-side SCC "// & + "iteration limit, matching native tblite's --iterations setting for reference CLI "// & + "diagnostics. For explicit SCC_MIXER TBLITE runs, CP2K synchronizes DFT/SCF/MAX_SCF to "// & + "this value. For SCC_MIXER AUTO or CP2K, DFT/SCF/MAX_SCF remains the normal CP2K "// & + "SCF-loop iteration limit. ITERATIONS is independent of MEMORY, which controls the "// & + "Broyden history length.", & + usage="ITERATIONS 250", default_i_val=tblite_mixer_iterations_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MEMORY", & + description="Number of previous modified-Broyden update vectors kept by the tblite SCC "// & + "mixer. This is a history length, not an iteration limit. Larger values retain more "// & + "history and use more memory; smaller values reduce the history and can be more robust for "// & + "difficult SCC cycles. A value of zero means library default: inherit ITERATIONS, as native "// & + "tblite does when it constructs the Broyden mixer from calc%max_iter.", & + usage="MEMORY 250", default_i_val=tblite_mixer_memory_inherit) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="SOLVER", & + description="Electronic eigensolver selected for native tblite SCC calculations where "// & + "tblite owns the electronic solve, and forwarded to native REFERENCE_CLI diagnostics as "// & + "tblite run --solver. GVD is tblite's default LAPACK divide-and-conquer generalized "// & + "symmetric eigensolver. GVR selects LAPACK's relatively robust generalized symmetric "// & + "eigensolver. Production CP2K/tblite calculations still use CP2K's electronic "// & + "diagonalization machinery.", & + usage="SOLVER (GVD|GVR)", & + enum_c_vals=s2a("GVD", "GVR"), & + enum_i_vals=[tblite_solver_gvd, tblite_solver_gvr], & + enum_desc=s2a("Use LAPACK's divide-and-conquer generalized symmetric eigensolver.", & + "Use LAPACK's relatively robust generalized symmetric eigensolver."), & + default_i_val=tblite_solver_gvd) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DAMPING", & + description="Damping parameter for the tblite modified Broyden SCC mixer. It controls the "// & + "size of the population/multipole update retained by the native tblite mixer between SCC "// & + "iterations. Smaller values damp oscillatory SCC cycles more strongly but may slow "// & + "convergence; larger values can accelerate well-behaved cases but may destabilize difficult "// & + "ones.", & + usage="DAMPING 0.4", default_r_val=tblite_mixer_damping_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="OMEGA0", & + description="Regularization weight added to the diagonal of tblite's modified-Broyden "// & + "linear system. The default matches the hard-coded value in native tblite and damps nearly "// & + "singular history updates without changing the first simple-damping step.", & + usage="OMEGA0 0.01", default_r_val=tblite_mixer_omega0_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MIN_WEIGHT", & + description="Lower bound for the dynamic Broyden history weight derived from the current "// & + "SCC residual. The default matches native tblite and prevents tiny weights from removing "// & + "history vectors from the update.", & + usage="MIN_WEIGHT 1.0", default_r_val=tblite_mixer_min_weight_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="MAX_WEIGHT", & + description="Upper bound for the dynamic Broyden history weight derived from the current "// & + "SCC residual. The default matches native tblite and caps the influence of very small "// & + "residual norms.", & + usage="MAX_WEIGHT 100000.0", default_r_val=tblite_mixer_max_weight_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="WEIGHT_FACTOR", & + description="Residual-to-weight scaling factor used by tblite's modified Broyden mixer. "// & + "For a residual norm above WEIGHT_FACTOR/MAX_WEIGHT, the dynamic weight is approximately "// & + "WEIGHT_FACTOR divided by that residual norm; otherwise MAX_WEIGHT is used before applying "// & + "MIN_WEIGHT. The default matches native tblite.", & + usage="WEIGHT_FACTOR 0.01", default_r_val=tblite_mixer_weight_factor_default) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_tblite_mixer_section + ! ************************************************************************************************** !> \brief Adds common SCC-mixer keywords for TB methods. !> \param section section receiving the keywords @@ -649,8 +823,21 @@ CONTAINS NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="SCC_MIXER", & description="Selects the SCC variable mixer for "//TARGET//" calculations. "// & - auto_desc//" TBLITE forces the tblite modified Broyden SCC mixer, CP2K forces CP2K "// & - "charge mixing, and NONE updates SCC variables without mixing.", & + auto_desc//" The SCC variables are the model populations and, where required by the "// & + "selected method, atomic multipoles that feed the self-consistent xTB potential. TBLITE "// & + "uses tblite's modified Broyden mixer for those variables inside the CP2K SCF cycle. CP2K "// & + "maps the variables onto CP2K's charge-mixing infrastructure. When TBLITE is active, the "// & + "native SCC convergence follows XTB/TBLITE/ACCURACY and the tblite-side SCC update limit "// & + "follows TBLITE_MIXER/ITERATIONS. If SCC_MIXER is explicitly set to TBLITE, CP2K also "// & + "sets DFT/SCF/MAX_SCF to TBLITE_MIXER/ITERATIONS. With SCC_MIXER AUTO or CP2K, MAX_SCF "// & + "keeps its normal CP2K meaning. When CP2K is active, convergence follows the normal CP2K SCF "// & + "controls, in particular DFT/SCF/EPS_SCF and DFT/SCF/MIXING. For GFN2/tblite, the CP2K "// & + "mixer acts on the full SCC-variable set: shell charges plus atomic dipole and quadrupole "// & + "variables. For CP2K-internal GFN1 and DFTB/DFTB3, the SCC variables are only charge or "// & + "shell-charge like. The "// & + "initial SCC variables are seeded from CP2K's current density, so DFT/SCF/SCF_GUESS and "// & + "DFT/WFN_RESTART_FILE_NAME define the starting point. NONE performs direct updates without "// & + "mixing and is intended for debugging, since it can be unstable.", & usage="SCC_MIXER (AUTO|TBLITE|CP2K|NONE)", & enum_c_vals=s2a("AUTO", "TBLITE", "CP2K", "NONE"), & enum_i_vals=[tblite_scc_mixer_auto, tblite_scc_mixer_tblite, & @@ -663,14 +850,6 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="TBLITE_MIXER_DAMPING", & - description="Damping parameter for the tblite modified Broyden SCC mixer. "// & - "Only active when the effective SCC mixer is TBLITE, either explicitly or via AUTO. "// & - "The iteration limit is taken from DFT/SCF/MAX_SCF.", & - usage="TBLITE_MIXER_DAMPING 0.4", default_r_val=0.4_dp) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - END SUBROUTINE add_tb_scc_mixer_keywords ! ************************************************************************************************** @@ -681,88 +860,603 @@ CONTAINS TYPE(section_type), POINTER :: section TYPE(keyword_type), POINTER :: keyword + TYPE(section_type), POINTER :: subsection CPASSERT(.NOT. ASSOCIATED(section)) CALL section_create(section, __LOCATION__, name="REFERENCE_CLI", & description="Run the native tblite command line interface as a reference "// & "for CP2K/tblite energies, nuclear gradients, and optionally virials. "// & - "The section is a lone section: if it is present, the reference check is enabled.", & + "The section is a lone section: if it is present, the reference check is enabled. By default "// & + "the native CLI calculation inherits the CP2K/tblite method, PARAM file, ACCURACY, charge, "// & + "spin state, spin-polarized Hamiltonian setting, electronic temperature, "// & + "TBLITE_MIXER/ITERATIONS, and TBLITE_MIXER/SOLVER from the surrounding input. CP2K forwards "// & + "DFT/UKS, including its LSD and SPIN_POLARIZED aliases, as native tblite "// & + "--spin-polarized; DFT/MULTIPLICITY is forwarded separately as --spin. "// & + "REFERENCE_CLI/EFIELD can add a homogeneous electric field to the "// & + "external native tblite command only. REFERENCE_CLI/IMPLICIT_SOLVATION can add native "// & + "tblite implicit-solvent settings to the external reference command only. "// & + "Native tblite run does not expose independent options for mixer MEMORY, DAMPING, OMEGA0, "// & + "MIN_WEIGHT, MAX_WEIGHT, or WEIGHT_FACTOR; CP2K warns if these inherited production settings "// & + "cannot be reproduced and aborts when STOP_ON_ERROR is enabled. "// & + "Native CLI verbosity follows GLOBAL/PRINT_LEVEL: SILENT is forwarded as --silent, and "// & + "HIGH or DEBUG are forwarded as --verbose. "// & + "It remains a diagnostic one-shot calculation from temporary input files and is independent "// & + "of CP2K's SCF history, diagonalization, and k-point machinery. If RESTART "// & + "is set, the named native tblite restart file is passed to the external command. Non-Fermi "// & + "CP2K smearing methods cannot be represented by native tblite and are replaced by a "// & + "Fermi-Dirac electronic temperature in the reference command with an explicit warning.", & n_keywords=1, n_subsections=0, repeats=.FALSE.) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & - description="Activates the native tblite CLI reference check when the section is present.", & + description="Activates the native tblite CLI reference check when the section is present. "// & + "The check compares CP2K/tblite against a freshly executed external tblite calculation.", & lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="PROGRAM_NAME", & - description="Executable name or path for the native tblite command line program.", & + description="Executable name or path for the native tblite command line program used for "// & + "the reference calculation. This allows testing the CP2K/tblite implementation against a "// & + "specific tblite build without changing the CP2K executable.", & usage="PROGRAM_NAME tblite", n_var=1, type_of_var=char_t, default_c_val="tblite") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="GUESS", & + description="Initial population guess used only for the native tblite reference command. "// & + "SAD uses the superposition of atomic densities, matching the default of tblite run. EEQ "// & + "uses electronegativity equilibration charges, while CEH uses the charge-equilibration "// & + "Hamiltonian guess. This keyword does not affect the production CP2K/tblite calculation; "// & + "there, initial densities and restarts remain controlled by CP2K's SCF machinery.", & + usage="GUESS (SAD|EEQ|CEH)", & + enum_c_vals=s2a("SAD", "EEQ", "CEH"), & + enum_i_vals=[tblite_guess_sad, tblite_guess_eeq, tblite_guess_ceh], & + enum_desc=s2a("Use the superposition of atomic densities guess.", & + "Use the electronegativity equilibration guess.", & + "Use the charge-equilibration Hamiltonian guess."), & + default_i_val=tblite_guess_sad) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="WORK_DIRECTORY", & - description="Directory used for temporary native tblite CLI input and output files.", & + description="Directory used for the temporary native tblite CLI input and output files. "// & + "The CP2K/tblite result is not taken from this directory; the files are only used to run and "// & + "inspect the independent reference calculation.", & usage="WORK_DIRECTORY .", n_var=1, type_of_var=char_t, default_c_val=".") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="PREFIX", & - description="Prefix for temporary native tblite CLI input and output files.", & + description="Prefix for the temporary native tblite CLI input and output files. Changing "// & + "the prefix is useful when several reference checks are run in the same directory or when "// & + "KEEP_FILES is used to inspect the generated geometry, gradient, JSON, or log files.", & usage="PREFIX tblite-reference", n_var=1, type_of_var=char_t, & default_c_val="tblite-reference") CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="INPUT_FORMAT", & + description="Input format hint passed to native tblite run as --input. The default is GEN, "// & + "matching the temporary DFTB+ GEN geometry written by CP2K for the reference calculation. "// & + "Only change this together with a compatible externally inspected or modified input file "// & + "workflow; the normal CP2K-generated reference geometry is GEN.", & + usage="INPUT_FORMAT gen", n_var=1, type_of_var=char_t, default_c_val="gen") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="RESTART", & + description="Native tblite restart file passed to the external reference command. When this "// & + "keyword is absent, CP2K passes --no-restart so each diagnostic reference is a clean one-shot "// & + "run from the generated geometry and selected CLI options. When a filename is given, CP2K "// & + "passes --restart with that file. This affects only the external native tblite run; it does "// & + "not read or write CP2K wavefunction or density restart data and is therefore independent of "// & + "DFT/SCF/SCF_GUESS RESTART.", & + usage="RESTART filename", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="POST_PROCESSING", & + description="Additional native tblite post-SCF analysis passed as tblite run "// & + "--post-processing. The value can be a built-in post-processing module name or a TOML file. "// & + "Native tblite already enables Mayer-Wiberg bond-order analysis internally by default. "// & + "The additional built-in module names currently exposed by tblite include molmom for "// & + "molecular moments and xtbml for xTB-ML descriptors. If the value names a TOML file, tblite "// & + "reads its [post-processing] configuration and can select the post-processing modules more "// & + "finely. This keyword affects only the external REFERENCE_CLI diagnostic run; CP2K does not "// & + "parse the generated post-processing data for the energy/force/virial comparison.", & + usage="POST_PROCESSING molmom", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="POST_PROCESSING_OUTPUT", & + description="Filename passed to native tblite as --post-processing-output. Native tblite "// & + "documents this as the NPZ output file for post-processing data, with tblite-data.npz as its "// & + "CLI default. In the locally checked tblite 0.5.0 source this option appears in the help text "// & + "but is not wired in the command-line parser, so CP2K emits it only when explicitly requested "// & + "and prints a warning before running the reference command. Use this keyword only with tblite "// & + "CLI builds that actually accept the option. Unless KEEP_FILES is enabled, CP2K removes this "// & + "explicit output file after the reference check.", & + usage="POST_PROCESSING_OUTPUT tblite-data.npz", n_var=1, type_of_var=char_t, & + default_c_val="", lone_keyword_c_val="tblite-data.npz") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="GRAD", & + description="Filename passed to native tblite as --grad for the tagged ASCII gradient output. "// & + "The native file contains the total energy, energy contributions, analytical molecular "// & + "gradient, and virial. CP2K reads this file for the REFERENCE_CLI energy, force, and virial "// & + "checks. If the keyword is omitted, CP2K writes a temporary PREFIX.grad file in "// & + "WORK_DIRECTORY. If the keyword is present without a filename, CP2K uses tblite's CLI default "// & + "name tblite.txt. Unless KEEP_FILES is enabled, CP2K removes the selected file after the "// & + "reference check.", & + usage="GRAD tblite.txt", n_var=1, type_of_var=char_t, default_c_val="", & + lone_keyword_c_val="tblite.txt") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="JSON", & + description="Filename passed to native tblite as --json for JSON result output. Without "// & + "--grad the native JSON usually contains the energy and energy contributions; because CP2K's "// & + "REFERENCE_CLI path also requests --grad, native tblite can additionally include derivative "// & + "data such as gradient and virial when supported by the CLI version. CP2K currently reads the "// & + "tagged GRAD file for comparisons and keeps the JSON file as an inspectable native artifact. "// & + "If the keyword is omitted, CP2K writes a temporary PREFIX.json file in WORK_DIRECTORY. If the "// & + "keyword is present without a filename, CP2K uses tblite's CLI default name tblite.json. "// & + "Unless KEEP_FILES is enabled, CP2K removes the selected file after the reference check.", & + usage="JSON tblite.json", n_var=1, type_of_var=char_t, default_c_val="", & + lone_keyword_c_val="tblite.json") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ELECTRONIC_TEMPERATURE_GUESS", & + description="Electronic temperature used only for the native tblite reference command's "// & + "initial CEH population guess. CP2K forwards positive values as tblite run "// & + "--etemp-guess. The value does not affect CP2K/tblite production calculations and is "// & + "independent of DFT/SCF/SMEAR/ELECTRONIC_TEMPERATURE, which controls the final electronic "// & + "occupation temperature passed as --etemp. A value of zero leaves the native tblite default "// & + "for the guess temperature untouched. Positive values are valid only together with GUESS CEH; "// & + "CP2K aborts for SAD or EEQ because native tblite uses this option for the CEH guess.", & + usage="ELECTRONIC_TEMPERATURE_GUESS [K] 4000", & + default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="K"), unit_str="K") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EFIELD", & + description="Homogeneous electric field passed only to the native tblite reference command "// & + "as tblite run --efield. The three components are given in V/Angstrom, matching the native "// & + "tblite CLI convention. This keyword does not affect the production CP2K/tblite "// & + "calculation because that path currently does not construct tblite's external electric "// & + "field container. If enabled together with energy, force, or virial checks, the comparison "// & + "therefore intentionally diagnoses the native field response rather than matching a CP2K "// & + "production field calculation.", & + usage="EFIELD 0.0 0.0 0.0", repeats=.FALSE., n_var=3, type_of_var=real_t, & + default_r_vals=[0.0_dp, 0.0_dp, 0.0_dp]) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="KEEP_FILES", & - description="Keep the temporary native tblite CLI input and output files.", & + description="Keep the temporary native tblite CLI input and output files after the "// & + "comparison. This is mainly a debugging aid for inspecting the generated GEN geometry, the "// & + "native tblite gradient/virial file, the JSON output, and the CLI log.", & usage="KEEP_FILES T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="ERROR_LIMIT", & - description="Tolerance for reporting reference CLI deviations.", & + description="Tolerance for reporting reference CLI deviations. The value is applied to the "// & + "absolute CP2K-minus-native differences printed by the enabled energy, force, and virial "// & + "checks. Exceeding this limit produces a diagnostic warning, or an abort when STOP_ON_ERROR "// & + "is enabled.", & usage="ERROR_LIMIT 1.0E-8", default_r_val=1.0E-8_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="STOP_ON_ERROR", & - description="Abort if a checked reference CLI deviation exceeds ERROR_LIMIT.", & + description="Abort the CP2K run if any enabled reference CLI comparison exceeds "// & + "ERROR_LIMIT. When disabled, CP2K prints the deviation but continues, which is useful for "// & + "exploratory diagnostics and regression-test probes.", & usage="STOP_ON_ERROR T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="CHECK_ENERGY", & - description="Compare CP2K/tblite energy against the native tblite CLI energy.", & + description="Compare the CP2K/tblite total xTB energy with the energy obtained from the "// & + "native tblite reference run. This checks the Hamiltonian setup, SCC solution, static "// & + "repulsion/dispersion terms, charge and spin state, electronic temperature, and optional "// & + "PARAM file consistency for the one-shot reference geometry.", & usage="CHECK_ENERGY T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="CHECK_FORCES", & - description="Compare CP2K/tblite nuclear gradients against native tblite CLI gradients.", & + description="Compare CP2K/tblite nuclear gradients with the gradients written by the "// & + "native tblite reference run. This is a sensitive check of derivative terms, including "// & + "repulsion, dispersion, electrostatic response, and coordination-number response. The "// & + "comparison uses the same geometry and model setup as the energy reference.", & usage="CHECK_FORCES T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="CHECK_VIRIAL", & - description="Compare CP2K/tblite virial against the native tblite CLI virial.", & + description="Compare CP2K/tblite virial with the virial written by the native tblite "// & + "reference run. This is a diagnostic for cell derivatives and stress-related terms. It is "// & + "disabled by default because virial conventions, finite periodicity, and the absence of "// & + "native CLI k-point sampling make this comparison more restrictive than the energy and "// & + "force checks.", & usage="CHECK_VIRIAL T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="ACCURACY", & - description="Value passed to native tblite as --acc.", & - usage="ACCURACY 1.0E-8", default_r_val=1.0E-8_dp) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) + NULLIFY (subsection) + CALL create_xtb_reference_cli_guess_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) - CALL keyword_create(keyword, __LOCATION__, name="ITERATIONS", & - description="Value passed to native tblite as --iterations.", & - usage="ITERATIONS 250", default_i_val=250) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) + NULLIFY (subsection) + CALL create_xtb_reference_cli_param_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + NULLIFY (subsection) + CALL create_xtb_reference_cli_fit_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + NULLIFY (subsection) + CALL create_xtb_reference_cli_tagdiff_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + NULLIFY (subsection) + CALL create_xtb_reference_cli_solvation_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) END SUBROUTINE create_xtb_reference_cli_section +! ************************************************************************************************** +!> \brief Creates the &REFERENCE_CLI%GUESS_CLI section for native tblite guess. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_xtb_reference_cli_guess_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="GUESS_CLI", & + description="Run the native tblite guess subcommand as an auxiliary REFERENCE_CLI "// & + "diagnostic. The command executes tblite guess on the generated GEN geometry by default and "// & + "prints the initializing charges and dipole moments. CP2K forwards the surrounding charge "// & + "and multiplicity as native --charge and --spin. GLOBAL/PRINT_LEVEL controls native "// & + "--silent or --verbose exactly as for the main REFERENCE_CLI run command. This section is "// & + "diagnostic only; CP2K does not use the generated guess data to seed the production "// & + "calculation.", & + n_keywords=1, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Activates the native tblite guess auxiliary command when the section is "// & + "present.", lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="METHOD", & + description="Initial-population method passed to native tblite guess as --method. SAD uses "// & + "superposed atomic densities, EEQ uses electronegativity equilibration, and CEH uses the "// & + "charge-extended Huckel guess. The default matches native tblite guess.", & + usage="METHOD (SAD|EEQ|CEH)", & + enum_c_vals=s2a("SAD", "EEQ", "CEH"), & + enum_i_vals=[tblite_guess_sad, tblite_guess_eeq, tblite_guess_ceh], & + enum_desc=s2a("Use the superposition of atomic densities guess.", & + "Use the electronegativity equilibration guess.", & + "Use the charge-extended Huckel guess."), & + default_i_val=tblite_guess_ceh) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ELECTRONIC_TEMPERATURE_GUESS", & + description="Electronic temperature for native tblite guess --etemp-guess. This option is "// & + "used by the CEH guess in native tblite. A value of zero leaves tblite's own default "// & + "untouched.", usage="ELECTRONIC_TEMPERATURE_GUESS [K] 4000", & + default_r_val=cp_unit_to_cp2k(value=0.0_dp, unit_str="K"), unit_str="K") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="SOLVER", & + description="Electronic solver passed to native tblite guess as --solver. GVD is the "// & + "default generalized eigenvalue decomposition solver; GVR selects the alternative solver "// & + "exposed by tblite.", usage="SOLVER (GVD|GVR)", & + enum_c_vals=s2a("GVD", "GVR"), enum_i_vals=[tblite_solver_gvd, tblite_solver_gvr], & + enum_desc=s2a("Use the GVD electronic solver.", "Use the GVR electronic solver."), & + default_i_val=tblite_solver_gvd) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EFIELD", & + description="Homogeneous electric field passed to native tblite guess as --efield. The "// & + "three components are given in V/Angstrom, matching native tblite.", & + usage="EFIELD 0.0 0.0 0.0", repeats=.FALSE., n_var=3, type_of_var=real_t, & + default_r_vals=[0.0_dp, 0.0_dp, 0.0_dp]) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="GRAD", & + description="Request native tblite guess --grad, evaluating the analytic gradient of the "// & + "guessed charges.", usage="GRAD T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="JSON", & + description="Filename passed to native tblite guess as --json. If the keyword is present "// & + "without a filename, CP2K uses tblite's CLI default name tblite.json. If omitted, CP2K does "// & + "not request JSON output for the auxiliary guess command.", & + usage="JSON tblite.json", n_var=1, type_of_var=char_t, default_c_val="", & + lone_keyword_c_val="tblite.json") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="INPUT_FORMAT", & + description="Input format hint passed to native tblite guess as --input. The default is GEN, "// & + "matching the temporary DFTB+ GEN geometry written by CP2K for the enclosing reference "// & + "calculation.", & + usage="INPUT_FORMAT gen", n_var=1, type_of_var=char_t, default_c_val="gen") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="INPUT_FILE", & + description="Positional input structure for native tblite guess. The file format must match "// & + "GUESS_CLI/INPUT_FORMAT. If omitted, CP2K uses the temporary GEN geometry written for the "// & + "enclosing REFERENCE_CLI calculation.", & + usage="INPUT_FILE structure.gen", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_reference_cli_guess_section + +! ************************************************************************************************** +!> \brief Creates the &REFERENCE_CLI%PARAM_CLI section for native tblite param. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_xtb_reference_cli_param_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="PARAM_CLI", & + description="Run the native tblite param subcommand as an auxiliary REFERENCE_CLI "// & + "diagnostic. The command inspects, imports, exports, or converts tblite parameter data. "// & + "If no INPUT_FILE is supplied, CP2K selects the surrounding XTB/TBLITE/METHOD as native "// & + "--method. If METHOD is specified here, it overrides only this auxiliary command.", & + n_keywords=1, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Activates the native tblite param auxiliary command when the section is "// & + "present.", lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="METHOD", & + description="Base parametrization passed to native tblite param as --method. If omitted, "// & + "CP2K uses the surrounding XTB/TBLITE/METHOD.", & + usage="METHOD (GFN1|GFN2|IPEA1)", & + enum_c_vals=s2a("GFN1", "GFN2", "IPEA1"), & + enum_i_vals=[gfn1xtb, gfn2xtb, ipea1xtb], & + enum_desc=s2a("Use the GFN1-xTB base parametrization.", & + "Use the GFN2-xTB base parametrization.", & + "Use the IPEA1 base parametrization."), & + default_i_val=gfn2xtb) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="OUTPUT", & + description="Output parameter file passed to native tblite param as --output. If omitted, "// & + "native tblite writes to its default stream or file for the selected operation.", & + usage="OUTPUT parameters.toml", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="INPUT_FILE", & + description="Optional positional input parameter file for native tblite param. If omitted, "// & + "the command operates on the selected built-in METHOD.", & + usage="INPUT_FILE parameters.toml", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_reference_cli_param_section + +! ************************************************************************************************** +!> \brief Creates the &REFERENCE_CLI%FIT_CLI section for native tblite fit. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_xtb_reference_cli_fit_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="FIT_CLI", & + description="Run the native tblite fit subcommand as an auxiliary REFERENCE_CLI diagnostic. "// & + "The command optimizes or dry-runs tblite parameters from a starting PARAM_FILE and a fit "// & + "settings INPUT_FILE. Because native fitting can be expensive and writes its own outputs, "// & + "the section is inactive unless explicitly present.", & + n_keywords=1, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Activates the native tblite fit auxiliary command when the section is present.", & + lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="PARAM_FILE", & + description="Starting parameter file passed as the first positional argument to native "// & + "tblite fit.", usage="PARAM_FILE start.toml", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="INPUT_FILE", & + description="Fit-settings input file passed as the second positional argument to native "// & + "tblite fit.", usage="INPUT_FILE fit.toml", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DRY_RUN", & + description="Pass native tblite fit --dry-run. This expands and checks the fit setup "// & + "without starting the optimization.", usage="DRY_RUN T", default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="COPY", & + description="Filename passed to native tblite fit as --copy. Native tblite writes the full "// & + "representation of the fit input to this file, with defaults filled in and masks expanded.", & + usage="COPY expanded-fit.toml", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_reference_cli_fit_section + +! ************************************************************************************************** +!> \brief Creates the &REFERENCE_CLI%TAGDIFF_CLI section for native tblite tagdiff. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_xtb_reference_cli_tagdiff_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="TAGDIFF_CLI", & + description="Run the native tblite tagdiff auxiliary command. It compares an ACTUAL tagged "// & + "data file against a REFERENCE tagged data file and reports differences using the keys from "// & + "the reference file.", n_keywords=1, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Activates the native tblite tagdiff auxiliary command when the section is "// & + "present.", lone_keyword_l_val=.TRUE., default_l_val=.FALSE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ACTUAL", & + description="Actual tagged data file passed as the first positional argument to native "// & + "tblite tagdiff.", usage="ACTUAL actual.tag", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="REFERENCE", & + description="Reference tagged data file passed as the second positional argument to native "// & + "tblite tagdiff. Native tblite uses only keys present in this reference file.", & + usage="REFERENCE reference.tag", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="FIT", & + description="Pass native tblite tagdiff --fit, producing output suitable for tblite fit.", & + usage="FIT T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_reference_cli_tagdiff_section + +! ************************************************************************************************** +!> \brief Creates the &REFERENCE_CLI%IMPLICIT_SOLVATION section for native tblite. +!> \param section section to create +! ************************************************************************************************** + SUBROUTINE create_xtb_reference_cli_solvation_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="IMPLICIT_SOLVATION", & + description="Native tblite implicit-solvation settings for the external REFERENCE_CLI "// & + "diagnostic command. This section does not affect the production CP2K/tblite calculation. "// & + "The section is active when present; no additional enable keyword is needed. "// & + "CP2K aborts for fully three-dimensional periodic cells and prints a warning for one- or "// & + "two-dimensional periodic cells because native tblite's implicit-solvent models are mainly "// & + "molecular continuum-solvation diagnostics. "// & + "MODEL selects the native tblite run solvation switch: ALPB maps to --alpb, GBSA to --gbsa, "// & + "GBE to --gbe, GB to --gb, and CPCM to --cpcm. SOLVENT is the argument passed to that switch. "// & + "For ALPB and GBSA it must be a named solvent with tblite parameters; for GBE, GB, and CPCM "// & + "it can be either a named solvent or a dielectric constant written as a token such as 78.36. "// & + "BORN_KERNEL maps to --born-kernel and applies only to ALPB, GBSA, GBE, and GB; AUTO leaves "// & + "the native tblite model-dependent default untouched. SOLUTION_STATE maps to --solv-state "// & + "and is only supported for the parametrized named solvent ALPB/GBSA models; GSOLV leaves "// & + "the native tblite default untouched.", & + n_keywords=4, repeats=.FALSE.) + + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="MODEL", & + description="Implicit-solvation model passed to native tblite run. ALPB selects the "// & + "parametrized analytical linearized Poisson-Boltzmann model with CDS and solution-state "// & + "shift support. GBSA selects the parametrized generalized Born surface-area model with CDS "// & + "and solution-state shift support. GBE selects finite-epsilon generalized Born/ALPB "// & + "electrostatics without the named-solvent CDS and shift terms. GB selects generalized Born "// & + "electrostatics without the named-solvent CDS and shift terms. CPCM selects tblite's "// & + "conductor-like polarizable continuum model.", & + usage="MODEL (ALPB|GBSA|GBE|GB|CPCM)", & + enum_c_vals=s2a("ALPB", "GBSA", "GBE", "GB", "CPCM"), & + enum_i_vals=[tblite_cli_solvation_alpb, tblite_cli_solvation_gbsa, & + tblite_cli_solvation_gbe, tblite_cli_solvation_gb, & + tblite_cli_solvation_cpcm], & + enum_desc=s2a("Use native tblite --alpb.", & + "Use native tblite --gbsa.", & + "Use native tblite --gbe.", & + "Use native tblite --gb.", & + "Use native tblite --cpcm."), & + default_i_val=tblite_cli_solvation_alpb) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="SOLVENT", & + description="Solvent argument passed to the selected native tblite solvation switch. "// & + "For MODEL ALPB and MODEL GBSA this must be a named solvent supported by tblite's "// & + "parametrized solvation data, for example water, acetone, methanol, thf, or toluene. "// & + "For MODEL GBE, MODEL GB, and MODEL CPCM this can also be a dielectric constant, for "// & + "example 78.36. Native tblite normalizes names by ignoring case, spaces, commas, and "// & + "hyphens before lookup.", & + usage="SOLVENT water", n_var=1, type_of_var=char_t, default_c_val="") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="BORN_KERNEL", & + description="Born interaction kernel passed to native tblite as --born-kernel. AUTO leaves "// & + "the native default untouched: tblite uses p16 for ALPB/GBE and still for GB/GBSA. P16 "// & + "selects the Lange p16 kernel, while STILL selects the classical Still kernel. This keyword "// & + "is not valid with MODEL CPCM because CPCM does not use the Born-kernel machinery.", & + usage="BORN_KERNEL (AUTO|P16|STILL)", & + enum_c_vals=s2a("AUTO", "P16", "STILL"), & + enum_i_vals=[tblite_cli_born_kernel_auto, tblite_cli_born_kernel_p16, & + tblite_cli_born_kernel_still], & + enum_desc=s2a("Use tblite's model-dependent default kernel.", & + "Use the Lange p16 Born kernel.", & + "Use the classical Still Born kernel."), & + default_i_val=tblite_cli_born_kernel_auto) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="SOLUTION_STATE", & + description="Solution-state convention passed to native tblite as --solv-state. GSOLV is "// & + "the default solvation free-energy convention. BAR1MOL applies the 1 bar ideal-gas to "// & + "1 mol/L liquid-solution shift. REFERENCE requests tblite's reference-state convention. "// & + "Native tblite supports non-default solution-state shifts only for named parametrized "// & + "ALPB/GBSA solvents, not for GBE, GB, or CPCM.", & + usage="SOLUTION_STATE (GSOLV|BAR1MOL|REFERENCE)", & + enum_c_vals=s2a("GSOLV", "BAR1MOL", "REFERENCE"), & + enum_i_vals=[tblite_cli_solution_state_gsolv, tblite_cli_solution_state_bar1mol, & + tblite_cli_solution_state_reference], & + enum_desc=s2a("Use tblite's default solvation free-energy convention.", & + "Apply the 1 bar ideal-gas to 1 mol/L solution shift.", & + "Use tblite's reference-state convention."), & + default_i_val=tblite_cli_solution_state_gsolv) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_reference_cli_solvation_section + END MODULE input_cp2k_tb diff --git a/src/qs_charge_mixing.F b/src/qs_charge_mixing.F index a898e19d70..a8ec820d55 100644 --- a/src/qs_charge_mixing.F +++ b/src/qs_charge_mixing.F @@ -10,9 +10,15 @@ MODULE qs_charge_mixing #if defined(__TBLITE) USE mctc_env, ONLY: error_type - USE tblite_scf, ONLY: new_mixer + USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer #endif - USE input_constants, ONLY: tblite_scc_mixer_auto, & + USE input_constants, ONLY: tblite_mixer_damping_default, & + tblite_mixer_iterations_default, & + tblite_mixer_max_weight_default, & + tblite_mixer_min_weight_default, & + tblite_mixer_omega0_default, & + tblite_mixer_weight_factor_default, & + tblite_scc_mixer_auto, & tblite_scc_mixer_cp2k, & tblite_scc_mixer_none, & tblite_scc_mixer_tblite @@ -46,28 +52,41 @@ CONTAINS !> \param para_env ... !> \param iter_count ... !> \param scc_mixer ... +!> \param tblite_mixer_iterations ... !> \param tblite_mixer_damping ... !> \param tblite_mixer_memory ... +!> \param tblite_mixer_omega0 ... +!> \param tblite_mixer_min_weight ... +!> \param tblite_mixer_max_weight ... +!> \param tblite_mixer_weight_factor ... !> \par History !> \author JGH ! ************************************************************************************************** SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, & - scc_mixer, tblite_mixer_damping, tblite_mixer_memory) + scc_mixer, tblite_mixer_iterations, tblite_mixer_damping, & + tblite_mixer_memory, tblite_mixer_omega0, tblite_mixer_min_weight, & + tblite_mixer_max_weight, tblite_mixer_weight_factor) INTEGER, INTENT(IN) :: mixing_method TYPE(mixing_storage_type), POINTER :: mixing_store REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges TYPE(mp_para_env_type), POINTER :: para_env INTEGER, INTENT(IN) :: iter_count - INTEGER, INTENT(IN), OPTIONAL :: scc_mixer + INTEGER, INTENT(IN), OPTIONAL :: scc_mixer, tblite_mixer_iterations REAL(KIND=dp), INTENT(IN), OPTIONAL :: tblite_mixer_damping INTEGER, INTENT(IN), OPTIONAL :: tblite_mixer_memory + REAL(KIND=dp), INTENT(IN), OPTIONAL :: tblite_mixer_omega0, & + tblite_mixer_min_weight, & + tblite_mixer_max_weight, & + tblite_mixer_weight_factor CHARACTER(len=*), PARAMETER :: routineN = 'charge_mixing' INTEGER :: effective_scc_mixer, handle, ia, ii, & - imin, inow, mixer_memory, nbuffer, ns, & - nvec - REAL(dp) :: alpha, mixer_damping + imin, inow, mixer_iterations, & + mixer_memory, nbuffer, ns, nvec + REAL(dp) :: alpha, mixer_damping, mixer_max_weight, & + mixer_min_weight, mixer_omega0, & + mixer_weight_factor CALL timeset(routineN, handle) @@ -81,12 +100,33 @@ CONTAINS CASE (tblite_scc_mixer_tblite) CPASSERT(ASSOCIATED(mixing_store)) #if defined(__TBLITE) - mixer_damping = 0.4_dp + mixer_damping = tblite_mixer_damping_default IF (PRESENT(tblite_mixer_damping)) mixer_damping = tblite_mixer_damping + IF (mixer_damping <= 0.0_dp) CPABORT("tblite SCC mixer DAMPING must be positive") + mixer_omega0 = tblite_mixer_omega0_default + IF (PRESENT(tblite_mixer_omega0)) mixer_omega0 = tblite_mixer_omega0 + IF (mixer_omega0 <= 0.0_dp) CPABORT("tblite SCC mixer OMEGA0 must be positive") + mixer_min_weight = tblite_mixer_min_weight_default + IF (PRESENT(tblite_mixer_min_weight)) mixer_min_weight = tblite_mixer_min_weight + IF (mixer_min_weight <= 0.0_dp) CPABORT("tblite SCC mixer MIN_WEIGHT must be positive") + mixer_max_weight = tblite_mixer_max_weight_default + IF (PRESENT(tblite_mixer_max_weight)) mixer_max_weight = tblite_mixer_max_weight + IF (mixer_max_weight <= 0.0_dp) CPABORT("tblite SCC mixer MAX_WEIGHT must be positive") + IF (mixer_max_weight < mixer_min_weight) & + CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT") + mixer_weight_factor = tblite_mixer_weight_factor_default + IF (PRESENT(tblite_mixer_weight_factor)) mixer_weight_factor = tblite_mixer_weight_factor + IF (mixer_weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive") + mixer_iterations = tblite_mixer_iterations_default + IF (PRESENT(tblite_mixer_iterations)) mixer_iterations = tblite_mixer_iterations + IF (mixer_iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive") + IF (iter_count > mixer_iterations) CPABORT("tblite SCC mixer exceeded ITERATIONS") mixer_memory = MAX(1, mixing_store%nbuffer) - IF (PRESENT(tblite_mixer_memory)) mixer_memory = MAX(1, tblite_mixer_memory) + IF (PRESENT(tblite_mixer_memory)) mixer_memory = tblite_mixer_memory + IF (mixer_memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive") CALL tblite_charge_mixing(mixing_store, charges, para_env, iter_count, & - mixer_damping, mixer_memory) + mixer_damping, mixer_memory, mixer_omega0, mixer_min_weight, & + mixer_max_weight, mixer_weight_factor) CALL timestop(handle) RETURN #else @@ -158,7 +198,7 @@ CONTAINS END SUBROUTINE charge_mixing ! ************************************************************************************************** -!> \brief Return the tblite SCC-mixer residual on the CP2K EPS_SCF scale. +!> \brief Return the CP2K-side tblite SCC-mixer residual on the CP2K EPS_SCF scale. !> \param mixing_store ... !> \param eps_scf ... !> \return ... @@ -188,13 +228,19 @@ CONTAINS !> \param iter_count ... !> \param damping ... !> \param memory ... +!> \param omega0 ... +!> \param min_weight ... +!> \param max_weight ... +!> \param weight_factor ... ! ************************************************************************************************** - SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory) + SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory, omega0, & + min_weight, max_weight, weight_factor) TYPE(mixing_storage_type), POINTER :: mixing_store REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges TYPE(mp_para_env_type), POINTER :: para_env INTEGER, INTENT(IN) :: iter_count, memory - REAL(KIND=dp), INTENT(IN) :: damping + REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, & + weight_factor #if defined(__TBLITE) TYPE(error_type), ALLOCATABLE :: error @@ -219,7 +265,8 @@ CONTAINS IF (reset_mixer) THEN IF (ALLOCATED(mixing_store%tb_scc_mixer)) DEALLOCATE (mixing_store%tb_scc_mixer) IF (on_source) THEN - CALL new_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping) + CALL new_cp2k_tblite_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping, omega0, & + min_weight, max_weight, weight_factor) CALL mixing_store%tb_scc_mixer%set(qvec) END IF mixing_store%tb_scc_mixer_natom = natom @@ -252,6 +299,10 @@ CONTAINS MARK_USED(iter_count) MARK_USED(damping) MARK_USED(memory) + MARK_USED(omega0) + MARK_USED(min_weight) + MARK_USED(max_weight) + MARK_USED(weight_factor) CPABORT("SCC_MIXER TBLITE requires CP2K to be built with tblite") #endif @@ -304,8 +355,7 @@ CONTAINS TYPE(mp_para_env_type), POINTER :: para_env INTEGER :: i, ia, ii, imin, j, nbuffer, nv - REAL(KIND=dp) :: alpha, maxw, minw, omega0, rskip, wdf, & - wfac + REAL(KIND=dp) :: alpha, broy_w0, rskip, wdf, wprod REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now @@ -338,20 +388,8 @@ CONTAINS END DO END IF - omega0 = 0.01_dp - minw = 1.0_dp - maxw = 100000.0_dp - wfac = 0.01_dp - - mixing_store%wbroy(nv) = SUM(dq_now(:, 1:ns)**2) - CALL para_env%sum(mixing_store%wbroy(nv)) - mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv)) - IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN - mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv) - ELSE - mixing_store%wbroy(nv) = maxw - END IF - IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw + broy_w0 = mixing_store%broy_w0 + mixing_store%wbroy(nv) = 1.0_dp ! dfbroy mixing_store%dfbroy(:, :, nv) = 0.0_dp @@ -363,25 +401,25 @@ CONTAINS ! abroy matrix DO i = 1, nv - wfac = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv)) - CALL para_env%sum(wfac) - mixing_store%abroy(i, nv) = wfac - mixing_store%abroy(nv, i) = wfac + wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv)) + CALL para_env%sum(wprod) + mixing_store%abroy(i, nv) = wprod + mixing_store%abroy(nv, i) = wprod END DO ! broyden matrices ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv)) DO i = 1, nv - wfac = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns)) - CALL para_env%sum(wfac) - cvec(i) = mixing_store%wbroy(i)*wfac + wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns)) + CALL para_env%sum(wprod) + cvec(i) = mixing_store%wbroy(i)*wprod END DO DO i = 1, nv DO j = 1, nv beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i) END DO - beta(i, i) = beta(i, i) + omega0*omega0 + beta(i, i) = beta(i, i) + broy_w0 END DO rskip = 1.e-12_dp diff --git a/src/qs_dftb_matrices.F b/src/qs_dftb_matrices.F index 678920bde7..e0b09c9297 100644 --- a/src/qs_dftb_matrices.F +++ b/src/qs_dftb_matrices.F @@ -585,8 +585,17 @@ CONTAINS CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, & mix_charge, para_env, scf_env%iter_count, & scc_mixer=dft_control%qs_control%dftb_control%tblite_scc_mixer, & + tblite_mixer_iterations= & + dft_control%qs_control%dftb_control%tblite_mixer_iterations, & tblite_mixer_damping=dft_control%qs_control%dftb_control%tblite_mixer_damping, & - tblite_mixer_memory=qs_env%scf_control%max_scf) + tblite_mixer_memory=dft_control%qs_control%dftb_control%tblite_mixer_memory, & + tblite_mixer_omega0=dft_control%qs_control%dftb_control%tblite_mixer_omega0, & + tblite_mixer_min_weight= & + dft_control%qs_control%dftb_control%tblite_mixer_min_weight, & + tblite_mixer_max_weight= & + dft_control%qs_control%dftb_control%tblite_mixer_max_weight, & + tblite_mixer_weight_factor= & + dft_control%qs_control%dftb_control%tblite_mixer_weight_factor) mcharge(:) = mix_charge(:, 1) DEALLOCATE (mix_charge) END IF diff --git a/src/qs_environment.F b/src/qs_environment.F index 8395345828..6fbd644067 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -96,9 +96,9 @@ MODULE qs_environment do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rigpw, & do_method_rm1, do_method_xtb, do_qmmm_gauss, do_qmmm_swave, driver_run, ehrenfest, & general_roks, geo_opt_run, hden_atomic, kg_tnadd_embed_ri, linear_response_run, & - mol_dyn_run, rel_none, rel_trans_atom, smear_fermi_dirac, vdw_pairpot_dftd2, & - vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, wfi_gext_proj_nr, & - wfi_gext_proj_qtr_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, & + mol_dyn_run, rel_none, rel_trans_atom, smear_fermi_dirac, tblite_scc_mixer_tblite, & + vdw_pairpot_dftd2, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, vdw_pairpot_dftd4, & + wfi_gext_proj_nr, wfi_gext_proj_qtr_nr, wfi_linear_ps_method_nr, wfi_linear_wf_method_nr, & wfi_use_prev_wf_method_nr, xc_vdw_fun_none, xc_vdw_fun_nonloc, xc_vdw_fun_pairpot, & xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none USE input_section_types, ONLY: section_get_ivals,& @@ -1225,7 +1225,7 @@ CONTAINS CALL tb_init_geometry(qs_env, qs_env%tb_tblite) ! select tblite method CALL tb_set_calculator(qs_env%tb_tblite, xtb_control%tblite_method, & - xtb_control%tblite_accuracy) + xtb_control%tblite_accuracy, xtb_control%tblite_param_file) !set up wave function CALL tb_init_wf(qs_env%tb_tblite, dft_control) !get basis set @@ -1630,6 +1630,14 @@ CONTAINS END IF CALL scf_c_create(scf_control) CALL scf_c_read_parameters(scf_control, dft_section) + SELECT CASE (dft_control%qs_control%method_id) + CASE (do_method_dftb) + IF (dft_control%qs_control%dftb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) & + scf_control%max_scf = dft_control%qs_control%dftb_control%tblite_mixer_iterations + CASE (do_method_xtb) + IF (dft_control%qs_control%xtb_control%tblite_scc_mixer == tblite_scc_mixer_tblite) & + scf_control%max_scf = dft_control%qs_control%xtb_control%tblite_mixer_iterations + END SELECT ! Allocate the data structure for Quickstep energies CALL allocate_qs_energy(energy) diff --git a/src/qs_mixing_utils.F b/src/qs_mixing_utils.F index b4545931aa..a4dadcd868 100644 --- a/src/qs_mixing_utils.F +++ b/src/qs_mixing_utils.F @@ -40,7 +40,7 @@ MODULE qs_mixing_utils CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils' INTEGER, PARAMETER, PRIVATE :: max_dftb_scc_vars = 1, & - max_xtb_scc_vars = 14 + max_xtb_scc_vars = 28 PUBLIC :: mixing_allocate, mixing_init, charge_mixing_init, & self_consistency_check diff --git a/src/tblite_interface.F b/src/tblite_interface.F index 0d39a9cb24..e4d072ea79 100644 --- a/src/tblite_interface.F +++ b/src/tblite_interface.F @@ -25,7 +25,7 @@ MODULE tblite_interface USE tblite_data_spin, ONLY: get_spin_constant USE tblite_integral_multipole, ONLY: multipole_cgto, multipole_grad_cgto, maxl, msao USE tblite_integral_type, ONLY: integral_type, new_integral - USE tblite_scf, ONLY: get_mixer_dimension, new_mixer + USE tblite_scf, ONLY: get_mixer_dimension USE tblite_scf_info, ONLY: scf_info, atom_resolved, shell_resolved, & orbital_resolved, not_used USE tblite_scf_potential, ONLY: potential_type, new_potential, add_pot_to_h1 @@ -59,11 +59,24 @@ MODULE tblite_interface USE cp_log_handling, ONLY: cp_get_default_logger, & cp_logger_type, cp_logger_get_default_io_unit USE cp_output_handling, ONLY: cp_print_key_should_output, & - cp_print_key_unit_nr, cp_print_key_finished_output + cp_print_key_unit_nr, cp_print_key_finished_output, & + debug_print_level, high_print_level, silent_print_level USE cp_units, ONLY: cp_unit_from_cp2k - USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb, & + USE input_constants, ONLY: gfn1xtb, gfn2xtb, ipea1xtb, smear_energy_window, & + smear_fermi_dirac, smear_gaussian, smear_list, smear_mp, smear_mv, & + tblite_cli_born_kernel_auto, tblite_cli_born_kernel_p16, & + tblite_cli_born_kernel_still, tblite_cli_solution_state_bar1mol, & + tblite_cli_solution_state_gsolv, tblite_cli_solution_state_reference, & + tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, & + tblite_cli_solvation_gb, tblite_cli_solvation_gbe, & + tblite_cli_solvation_gbsa, & + tblite_guess_ceh, tblite_guess_eeq, tblite_guess_sad, & + tblite_mixer_damping_default, & + tblite_mixer_max_weight_default, tblite_mixer_min_weight_default, & + tblite_mixer_omega0_default, tblite_mixer_weight_factor_default, & tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, & - tblite_scc_mixer_none, tblite_scc_mixer_tblite + tblite_scc_mixer_none, tblite_scc_mixer_tblite, & + tblite_solver_gvd, tblite_solver_gvr USE input_section_types, ONLY: section_vals_val_get USE kinds, ONLY: default_path_length, default_string_length, dp USE kpoint_types, ONLY: get_kpoint_info, kpoint_type @@ -90,6 +103,9 @@ MODULE tblite_interface USE scf_control_types, ONLY: scf_control_type USE input_section_types, ONLY: section_vals_get_subs_vals, section_vals_type USE string_utilities, ONLY: integer_to_string +#if defined(__TBLITE) + USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer +#endif USE tblite_types, ONLY: tblite_type, deallocate_tblite_type, allocate_tblite_type USE virial_types, ONLY: virial_type USE xtb_types, ONLY: get_xtb_atom_param, xtb_atom_type @@ -403,9 +419,10 @@ CONTAINS info = tb%calc%variable_info() IF (ALLOCATED(tb%mixer)) DEALLOCATE (tb%mixer) - CALL new_mixer(tb%mixer, tb%calc%max_iter, & - tb%wfn%nspin*get_mixer_dimension(tb%mol, tb%calc%bas, info), & - tb%calc%mixer_damping) + CALL new_cp2k_tblite_mixer(tb%mixer, tb%mixer_memory, & + tb%wfn%nspin*get_mixer_dimension(tb%mol, tb%calc%bas, info), & + tb%mixer_damping, tb%mixer_omega0, tb%mixer_min_weight, & + tb%mixer_max_weight, tb%mixer_weight_factor) #else MARK_USED(tb) @@ -417,27 +434,60 @@ CONTAINS ! ************************************************************************************************** !> \brief Configure tblite's internal SCC mixer from CP2K input. !> \param tb ... -!> \param max_scf ... +!> \param iterations ... +!> \param memory ... !> \param damping ... +!> \param omega0 ... +!> \param min_weight ... +!> \param max_weight ... +!> \param weight_factor ... +!> \param solver ... ! ************************************************************************************************** - SUBROUTINE tb_configure_mixer(tb, max_scf, damping) + SUBROUTINE tb_configure_mixer(tb, iterations, memory, damping, omega0, min_weight, max_weight, & + weight_factor, solver) TYPE(tblite_type), POINTER :: tb - INTEGER, INTENT(IN) :: max_scf - REAL(KIND=dp), INTENT(IN) :: damping + INTEGER, INTENT(IN) :: iterations, memory, solver + REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, & + weight_factor #if defined(__TBLITE) - IF (max_scf < 1) CPABORT("tblite SCC mixer MAX_SCF must be positive") + IF (iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive") + IF (memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive") IF (damping <= 0.0_dp) CPABORT("tblite SCC mixer damping must be positive") + IF (omega0 <= 0.0_dp) CPABORT("tblite SCC mixer OMEGA0 must be positive") + IF (min_weight <= 0.0_dp) CPABORT("tblite SCC mixer MIN_WEIGHT must be positive") + IF (max_weight <= 0.0_dp) CPABORT("tblite SCC mixer MAX_WEIGHT must be positive") + IF (max_weight < min_weight) & + CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT") + IF (weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive") + SELECT CASE (solver) + CASE (tblite_solver_gvd, tblite_solver_gvr) + CASE DEFAULT + CPABORT("Unknown tblite SCC mixer SOLVER") + END SELECT - tb%calc%max_iter = max_scf + tb%calc%max_iter = iterations + tb%mixer_memory = memory + tb%mixer_solver = solver + tb%mixer_damping = damping tb%calc%mixer_damping = damping + tb%mixer_omega0 = omega0 + tb%mixer_min_weight = min_weight + tb%mixer_max_weight = max_weight + tb%mixer_weight_factor = weight_factor #else MARK_USED(tb) - MARK_USED(max_scf) + MARK_USED(iterations) + MARK_USED(memory) MARK_USED(damping) + MARK_USED(omega0) + MARK_USED(min_weight) + MARK_USED(max_weight) + MARK_USED(weight_factor) + MARK_USED(solver) CPABORT("Built without TBLITE") #endif @@ -470,10 +520,10 @@ CONTAINS END FUNCTION tb_native_scc_mixer_active ! ************************************************************************************************** -!> \brief Return the native tblite SCC mixer residual for CP2K SCF convergence. +!> \brief Return the native tblite SCC mixer residual on the CP2K iter_delta scale. !> \param dft_control ... !> \param tb ... -!> \param eps_scf ... +!> \param eps_scf CP2K reporting scale for the native residual. !> \return ... ! ************************************************************************************************** FUNCTION tb_scf_mixer_error(dft_control, tb, eps_scf) RESULT(mixer_error) @@ -512,27 +562,38 @@ CONTAINS !> \param tb ... !> \param typ ... !> \param accuracy ... +!> \param param_file ... ! ************************************************************************************************** - SUBROUTINE tb_set_calculator(tb, typ, accuracy) + SUBROUTINE tb_set_calculator(tb, typ, accuracy, param_file) TYPE(tblite_type), POINTER :: tb INTEGER :: typ REAL(KIND=dp), INTENT(IN) :: accuracy + CHARACTER(LEN=*), INTENT(IN) :: param_file #if defined(__TBLITE) TYPE(error_type), ALLOCATABLE :: error - SELECT CASE (typ) - CASE default - CPABORT("Unknown xtb type") - CASE (gfn1xtb) - CALL new_gfn1_calculator(tb%calc, tb%mol, error) - CASE (gfn2xtb) - CALL new_gfn2_calculator(tb%calc, tb%mol, error) - CASE (ipea1xtb) - CALL new_ipea1_calculator(tb%calc, tb%mol, error) - END SELECT + IF (ALLOCATED(tb%param)) DEALLOCATE (tb%param) + IF (LEN_TRIM(param_file) > 0) THEN + ALLOCATE (tb%param) + CALL tb%param%load(TRIM(param_file), error) + IF (ALLOCATED(error)) CPABORT("Could not load tblite PARAM file: "//TRIM(param_file)) + CALL new_xtb_calculator(tb%calc, tb%mol, tb%param, error) + ELSE + SELECT CASE (typ) + CASE default + CPABORT("Unknown xtb type") + CASE (gfn1xtb) + CALL new_gfn1_calculator(tb%calc, tb%mol, error) + CASE (gfn2xtb) + CALL new_gfn2_calculator(tb%calc, tb%mol, error) + CASE (ipea1xtb) + CALL new_ipea1_calculator(tb%calc, tb%mol, error) + END SELECT + END IF + IF (ALLOCATED(error)) CPABORT("tblite calculator setup failed") tb%accuracy = accuracy @@ -540,6 +601,7 @@ CONTAINS MARK_USED(tb) MARK_USED(typ) MARK_USED(accuracy) + MARK_USED(param_file) CPABORT("Built without TBLITE") #endif @@ -1279,9 +1341,11 @@ CONTAINS LOGICAL :: discard_mixed_output, do_combined_mixing, do_dipole, do_quadrupole, & native_sign_mixing, phased_combined_mixing, skip_charge_mixing, & reuse_native_input, skip_scf_dispersion, skip_scf_dispersion_energy, & + seed_native_from_rho, & skip_scf_dispersion_gradient, skip_scf_dispersion_potential, & use_native_mixer, use_no_mixer - REAL(KIND=dp) :: mix_alpha, norm, moment_alpha, new_charge, new_moment, pao + REAL(KIND=dp) :: mix_alpha, native_seed_charge, norm, moment_alpha, & + new_charge, new_moment, pao #if defined(__TBLITE_DEBUG_DIAGNOSTICS) INTEGER :: debug_status CHARACTER(LEN=32) :: debug_value @@ -1358,10 +1422,20 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set, & atomic_kind_set=atomic_kind_set, matrix_s_kp=matrix_s, rho=rho, para_env=para_env, & scf_control=scf_control) + IF (use_native_mixer .AND. use_rho .AND. (.NOT. calculate_forces) .AND. & + scf_env%iter_count > dft_control%qs_control%xtb_control%tblite_mixer_iterations) THEN + CPABORT("tblite SCC mixer exceeded TBLITE_MIXER/ITERATIONS") + END IF IF (use_native_mixer .AND. use_rho .AND. (.NOT. calculate_forces) .AND. & scf_env%iter_count == 1) THEN - CALL tb_configure_mixer(tb, scf_control%max_scf, & - dft_control%qs_control%xtb_control%tblite_mixer_damping) + CALL tb_configure_mixer(tb, dft_control%qs_control%xtb_control%tblite_mixer_iterations, & + dft_control%qs_control%xtb_control%tblite_mixer_memory, & + dft_control%qs_control%xtb_control%tblite_mixer_damping, & + dft_control%qs_control%xtb_control%tblite_mixer_omega0, & + dft_control%qs_control%xtb_control%tblite_mixer_min_weight, & + dft_control%qs_control%xtb_control%tblite_mixer_max_weight, & + dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, & + dft_control%qs_control%xtb_control%tblite_mixer_solver) CALL tb_reset_mixer(tb) END IF nspin = dft_control%nspins @@ -1377,8 +1451,6 @@ CONTAINS ELSE matrix_p => scf_env%p_mix_new END IF - IF (nspin > 1 .AND. (.NOT. use_native_mixer)) & - CPABORT("tblite UKS/LSD currently requires SCC_MIXER AUTO or TBLITE") IF (nspin > 1 .AND. mp_test_mode /= 0) & CPABORT("tblite UKS/LSD debug multipole modes are not implemented") IF (nspin > 1 .AND. (.NOT. calculate_forces)) CALL tb_store_density_ref(tb, matrix_p) @@ -1518,6 +1590,9 @@ CONTAINS END DO END DO END DO + native_seed_charge = -SUM(ch_atom(:, 1)) + seed_native_from_rho = SUM(ABS(aocg_spin)) > 1.0E-10_dp .AND. & + ABS(native_seed_charge - REAL(dft_control%charge, dp)) < 1.0E-5_dp DEALLOCATE (aocg, aocg_spin) raw_state_status = 1 @@ -1653,10 +1728,39 @@ CONTAINS mp_test_mode == 232 .OR. mp_test_mode == 233) .AND. & use_rho .AND. (.NOT. calculate_forces)) THEN IF (.NOT. reuse_native_input) THEN - tb%wfn%qsh(:, :) = 0.0_dp - tb%wfn%qat(:, :) = 0.0_dp - IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp - IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp + IF (seed_native_from_rho) THEN + ! Seed the native tblite mixer from CP2K's current density so SCF_GUESS/RESTART + ! controls the initial SCC state. + DO iatom = 1, n_atom + ii = tb%calc%bas%ish_at(iatom) + DO ispin = 1, nspin + DO is = 1, tb%calc%bas%nsh_at(iatom) + tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin) + END DO + tb%wfn%qat(iatom, ispin) = & + SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin)) + END DO + END DO + IF (do_dipole) THEN + DO iatom = 1, n_atom + DO ispin = 1, nspin + tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin) + END DO + END DO + END IF + IF (do_quadrupole) THEN + DO iatom = 1, n_atom + DO ispin = 1, nspin + tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin) + END DO + END DO + END IF + ELSE + tb%wfn%qsh(:, :) = 0.0_dp + tb%wfn%qat(:, :) = 0.0_dp + IF (do_dipole) tb%wfn%dpat(:, :, :) = 0.0_dp + IF (do_quadrupole) tb%wfn%qpat(:, :, :) = 0.0_dp + END IF END IF IF (do_dipole) DEALLOCATE (ao_dip) IF (do_quadrupole) DEALLOCATE (ao_quad) @@ -1701,8 +1805,66 @@ CONTAINS END IF ELSE ! charge mixing + native_sign_mixing = .FALSE. IF (dft_control%qs_control%do_ls_scf) THEN - ! + IF (nspin > 1) CPABORT("CP2K/tblite UKS/LSD with LS_SCF is not implemented") + ELSEIF (nspin > 1) THEN + n_mix_cols = nspin*max_shell + IF (do_dipole) n_mix_cols = n_mix_cols + nspin*dip_n + IF (do_quadrupole) n_mix_cols = n_mix_cols + nspin*quad_n + ALLOCATE (mix_vars(n_atom, n_mix_cols)) + mix_vars = 0.0_dp + + mix_offset = 0 + DO ispin = 1, nspin + mix_vars(:, mix_offset + 1:mix_offset + max_shell) = & + -ch_shell_spin(:, 1:max_shell, ispin) + mix_offset = mix_offset + max_shell + END DO + IF (do_dipole) THEN + DO ispin = 1, nspin + mix_vars(:, mix_offset + 1:mix_offset + dip_n) = & + -ao_dip_spin(:, 1:dip_n, ispin) + mix_offset = mix_offset + dip_n + END DO + END IF + IF (do_quadrupole) THEN + DO ispin = 1, nspin + mix_vars(:, mix_offset + 1:mix_offset + quad_n) = & + -ao_quad_spin(:, 1:quad_n, ispin) + mix_offset = mix_offset + quad_n + END DO + END IF + + IF (.NOT. use_no_mixer) THEN + CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, & + mix_vars, para_env, scf_env%iter_count) + END IF + + mix_offset = 0 + DO ispin = 1, nspin + ch_shell_spin(:, 1:max_shell, ispin) = & + -mix_vars(:, mix_offset + 1:mix_offset + max_shell) + mix_offset = mix_offset + max_shell + END DO + ch_shell(:, 1:max_shell) = ch_shell_spin(:, 1:max_shell, 1) + IF (do_dipole) THEN + DO ispin = 1, nspin + ao_dip_spin(:, 1:dip_n, ispin) = & + -mix_vars(:, mix_offset + 1:mix_offset + dip_n) + mix_offset = mix_offset + dip_n + END DO + ao_dip(:, 1:dip_n) = ao_dip_spin(:, 1:dip_n, 1) + END IF + IF (do_quadrupole) THEN + DO ispin = 1, nspin + ao_quad_spin(:, 1:quad_n, ispin) = & + -mix_vars(:, mix_offset + 1:mix_offset + quad_n) + mix_offset = mix_offset + quad_n + END DO + ao_quad(:, 1:quad_n) = ao_quad_spin(:, 1:quad_n, 1) + END IF + DEALLOCATE (mix_vars) ELSE ! CP2K SCC mixing has to see all tblite SCC variables, not only shell charges. do_combined_mixing = do_dipole .OR. do_quadrupole .OR. & @@ -1819,55 +1981,84 @@ CONTAINS moment_alpha = scf_env%p_mix_alpha IF (mp_test_mode == 61 .OR. mp_test_mode == 91) moment_alpha = 0.05_dp IF (mp_test_mode == 62 .OR. mp_test_mode == 92) moment_alpha = 0.5_dp - DO iatom = 1, n_atom - ii = tb%calc%bas%ish_at(iatom) - DO is = 1, tb%calc%bas%nsh_at(iatom) - new_charge = -ch_shell(iatom, is) - IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & - scf_env%iter_count > 1) THEN - new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qsh(ii + is, 1) - END IF - tb%wfn%qsh(ii + is, 1) = new_charge + IF (nspin > 1) THEN + DO iatom = 1, n_atom + ii = tb%calc%bas%ish_at(iatom) + DO ispin = 1, nspin + DO is = 1, tb%calc%bas%nsh_at(iatom) + tb%wfn%qsh(ii + is, ispin) = -ch_shell_spin(iatom, is, ispin) + END DO + tb%wfn%qat(iatom, ispin) = & + SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), ispin)) + END DO END DO - IF (native_sign_mixing .OR. mp_test_mode == 122) THEN - tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1)) - ELSE - new_charge = -ch_atom(iatom, 1) - IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & - scf_env%iter_count > 1) THEN - new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qat(iatom, 1) - END IF - tb%wfn%qat(iatom, 1) = new_charge + IF (do_dipole) THEN + DO iatom = 1, n_atom + DO ispin = 1, nspin + tb%wfn%dpat(:, iatom, ispin) = -ao_dip_spin(iatom, :, ispin) + END DO + END DO + DEALLOCATE (ao_dip) END IF - END DO + IF (do_quadrupole) THEN + DO iatom = 1, n_atom + DO ispin = 1, nspin + tb%wfn%qpat(:, iatom, ispin) = -ao_quad_spin(iatom, :, ispin) + END DO + END DO + DEALLOCATE (ao_quad) + END IF + ELSE + DO iatom = 1, n_atom + ii = tb%calc%bas%ish_at(iatom) + DO is = 1, tb%calc%bas%nsh_at(iatom) + new_charge = -ch_shell(iatom, is) + IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & + scf_env%iter_count > 1) THEN + new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qsh(ii + is, 1) + END IF + tb%wfn%qsh(ii + is, 1) = new_charge + END DO + IF (native_sign_mixing .OR. mp_test_mode == 122) THEN + tb%wfn%qat(iatom, 1) = SUM(tb%wfn%qsh(ii + 1:ii + tb%calc%bas%nsh_at(iatom), 1)) + ELSE + new_charge = -ch_atom(iatom, 1) + IF ((mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & + scf_env%iter_count > 1) THEN + new_charge = moment_alpha*new_charge + (1.0_dp - moment_alpha)*tb%wfn%qat(iatom, 1) + END IF + tb%wfn%qat(iatom, 1) = new_charge + END IF + END DO - IF (do_dipole) THEN - DO iatom = 1, n_atom - DO im = 1, dip_n - new_moment = -ao_dip(iatom, im) - IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. & - mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & - scf_env%iter_count > 1) THEN - new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%dpat(im, iatom, 1) - END IF - tb%wfn%dpat(im, iatom, 1) = new_moment + IF (do_dipole) THEN + DO iatom = 1, n_atom + DO im = 1, dip_n + new_moment = -ao_dip(iatom, im) + IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. & + mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & + scf_env%iter_count > 1) THEN + new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%dpat(im, iatom, 1) + END IF + tb%wfn%dpat(im, iatom, 1) = new_moment + END DO END DO - END DO - DEALLOCATE (ao_dip) - END IF - IF (do_quadrupole) THEN - DO iatom = 1, n_atom - DO im = 1, quad_n - new_moment = -ao_quad(iatom, im) - IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. & - mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & - scf_env%iter_count > 1) THEN - new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%qpat(im, iatom, 1) - END IF - tb%wfn%qpat(im, iatom, 1) = new_moment + DEALLOCATE (ao_dip) + END IF + IF (do_quadrupole) THEN + DO iatom = 1, n_atom + DO im = 1, quad_n + new_moment = -ao_quad(iatom, im) + IF ((mp_test_mode == 60 .OR. mp_test_mode == 61 .OR. mp_test_mode == 62 .OR. & + mp_test_mode == 90 .OR. mp_test_mode == 91 .OR. mp_test_mode == 92) .AND. & + scf_env%iter_count > 1) THEN + new_moment = moment_alpha*new_moment + (1.0_dp - moment_alpha)*tb%wfn%qpat(im, iatom, 1) + END IF + tb%wfn%qpat(im, iatom, 1) = new_moment + END DO END DO - END DO - DEALLOCATE (ao_quad) + DEALLOCATE (ao_quad) + END IF END IF END IF IF (mp_test_mode == 2 .OR. mp_test_mode == 7 .OR. mp_test_mode == 8 .OR. & @@ -4212,22 +4403,29 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'tb_reference_cli_compare' - CHARACTER(LEN=32) :: acc_str, charge_str, iter_str, spin_str, & - spinpol_str - CHARACTER(LEN=8) :: method + CHARACTER(LEN=16) :: solvation_model_name + CHARACTER(LEN=32) :: acc_str, charge_str, efield_x_str, efield_y_str, efield_z_str, & + etemp_guess_val_str, etemp_str, iter_str, spin_str, spinpol_str + CHARACTER(LEN=4*default_path_length+16) :: efield_str, etemp_guess_str, param_str, & + post_processing_output_str, post_processing_str, restart_str, solvation_str, verbosity_str + CHARACTER(LEN=8) :: guess, method, solver CHARACTER(LEN=8*default_path_length) :: command CHARACTER(LEN=default_path_length) :: file_base, gen_file, grad_file, & - json_file, log_file + json_file, log_file, & + post_processing_output_file INTEGER :: cmdstat, exitstat, handle, iounit, & - natom, nkp, spin + n_periodic, natom, nkp, & + reference_iterations, spin + INTEGER, DIMENSION(3) :: periodic LOGICAL :: do_kpoints, have_energy, have_gradient, & have_virial, too_large, & unsupported_kpoints - REAL(KIND=dp) :: cli_energy, cp_energy, ediff, fmax, & - fsum, vmax, vsum + REAL(KIND=dp) :: cli_energy, cp_energy, ediff, etemp, & + etemp_guess, fmax, fsum, vmax, vsum REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cli_gradient, cli_virial, cp_gradient REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell TYPE(cp_logger_type), POINTER :: logger TYPE(dft_control_type), POINTER :: dft_control TYPE(kpoint_type), POINTER :: kpoints @@ -4235,17 +4433,18 @@ CONTAINS TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_energy_type), POINTER :: energy TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(scf_control_type), POINTER :: scf_control TYPE(virial_type), POINTER :: virial TYPE(xtb_reference_cli_type) :: ref CALL timeset(routineN, handle) - NULLIFY (atomic_kind_set, dft_control, energy, force, kpoints, logger, para_env, particle_set, & - virial, xkp) - CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, & + NULLIFY (atomic_kind_set, cell, dft_control, energy, force, kpoints, logger, para_env, particle_set, & + scf_control, virial, xkp) + CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, & dft_control=dft_control, energy=energy, force=force, & - para_env=para_env, particle_set=particle_set, virial=virial, & - do_kpoints=do_kpoints, kpoints=kpoints) + para_env=para_env, particle_set=particle_set, scf_control=scf_control, & + virial=virial, do_kpoints=do_kpoints, kpoints=kpoints) ref = dft_control%qs_control%xtb_control%reference_cli IF (.NOT. ref%enabled) THEN @@ -4259,6 +4458,31 @@ CONTAINS logger => cp_get_default_logger() iounit = cp_logger_get_default_io_unit(logger) + verbosity_str = "" + IF (logger%iter_info%print_level == silent_print_level) THEN + verbosity_str = " --silent" + ELSE IF (logger%iter_info%print_level == high_print_level .OR. & + logger%iter_info%print_level == debug_print_level) THEN + verbosity_str = " --verbose" + END IF + IF (ref%solvation_active) THEN + periodic = 0 + IF (ASSOCIATED(cell)) CALL get_cell(cell=cell, periodic=periodic) + n_periodic = COUNT(periodic == 1) + IF (n_periodic == 3) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "tblite reference CLI implicit solvation is not supported for PERIODIC XYZ." + WRITE (UNIT=iounit, FMT="(T2,A)") & + "Use PERIODIC NONE for molecular solvation diagnostics, or remove IMPLICIT_SOLVATION." + CPABORT("REFERENCE_CLI implicit solvation is incompatible with PERIODIC XYZ") + ELSE IF (n_periodic > 0) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI implicit solvation with finite periodicity is diagnostic only." + WRITE (UNIT=iounit, FMT="(T2,A,I0,A)") & + "The generated native tblite reference geometry has ", n_periodic, & + " periodic direction(s); continuum-solvation conventions are primarily molecular." + END IF + END IF unsupported_kpoints = .FALSE. IF (do_kpoints .AND. ASSOCIATED(kpoints)) THEN nkp = 0 @@ -4277,30 +4501,195 @@ CONTAINS END IF natom = SIZE(particle_set) method = tb_reference_method_name(dft_control%qs_control%xtb_control%tblite_method) + guess = tb_reference_guess_name(ref%guess) + solver = tb_reference_solver_name(dft_control%qs_control%xtb_control%tblite_mixer_solver) file_base = tb_join_path(ref%work_directory, ref%prefix) gen_file = TRIM(file_base)//".gen" grad_file = TRIM(file_base)//".grad" json_file = TRIM(file_base)//".json" log_file = TRIM(file_base)//".log" + post_processing_output_file = "" + IF (LEN_TRIM(ref%grad_file) > 0) grad_file = ref%grad_file + IF (LEN_TRIM(ref%json_file) > 0) json_file = ref%json_file + IF (LEN_TRIM(ref%post_processing_output_file) > 0) & + post_processing_output_file = ref%post_processing_output_file WRITE (charge_str, "(I0)") dft_control%charge spin = MAX(0, dft_control%multiplicity - 1) WRITE (spin_str, "(I0)") spin - WRITE (acc_str, "(ES16.8)") ref%accuracy - WRITE (iter_str, "(I0)") ref%iterations + WRITE (acc_str, "(ES16.8)") dft_control%qs_control%xtb_control%tblite_accuracy + reference_iterations = dft_control%qs_control%xtb_control%tblite_mixer_iterations + WRITE (iter_str, "(I0)") reference_iterations + efield_str = "" + IF (ref%efield_active) THEN + WRITE (efield_x_str, "(ES16.8)") ref%efield(1) + WRITE (efield_y_str, "(ES16.8)") ref%efield(2) + WRITE (efield_z_str, "(ES16.8)") ref%efield(3) + efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// & + TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str)) + END IF + solvation_str = "" + solvation_model_name = "" + IF (ref%solvation_active) THEN + SELECT CASE (ref%solvation_model) + CASE (tblite_cli_solvation_alpb) + solvation_model_name = "ALPB" + solvation_str = " --alpb " + CASE (tblite_cli_solvation_gbsa) + solvation_model_name = "GBSA" + solvation_str = " --gbsa " + CASE (tblite_cli_solvation_gbe) + solvation_model_name = "GBE" + solvation_str = " --gbe " + CASE (tblite_cli_solvation_gb) + solvation_model_name = "GB" + solvation_str = " --gb " + CASE (tblite_cli_solvation_cpcm) + solvation_model_name = "CPCM" + solvation_str = " --cpcm " + CASE DEFAULT + CPABORT("Unknown tblite reference CLI implicit-solvation model") + END SELECT + solvation_str = TRIM(solvation_str)//" "//TRIM(tb_shell_quote(ref%solvation_solvent)) + SELECT CASE (ref%solvation_born_kernel) + CASE (tblite_cli_born_kernel_auto) + CASE (tblite_cli_born_kernel_p16) + solvation_str = TRIM(solvation_str)//" --born-kernel p16" + CASE (tblite_cli_born_kernel_still) + solvation_str = TRIM(solvation_str)//" --born-kernel still" + CASE DEFAULT + CPABORT("Unknown tblite reference CLI Born kernel") + END SELECT + SELECT CASE (ref%solvation_state) + CASE (tblite_cli_solution_state_gsolv) + CASE (tblite_cli_solution_state_bar1mol) + solvation_str = TRIM(solvation_str)//" --solv-state bar1mol" + CASE (tblite_cli_solution_state_reference) + solvation_str = TRIM(solvation_str)//" --solv-state reference" + CASE DEFAULT + CPABORT("Unknown tblite reference CLI solution state") + END SELECT + END IF + IF (dft_control%qs_control%xtb_control%tblite_mixer_memory /= reference_iterations) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MEMORY." + WRITE (UNIT=iounit, FMT="(T2,A,I0,A,I0,A)") & + "The native reference run uses tblite's internal mixer memory tied to --iterations (", & + reference_iterations, "), while CP2K uses MEMORY ", & + dft_control%qs_control%xtb_control%tblite_mixer_memory, "." + IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MEMORY") + END IF + IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_damping - & + tblite_mixer_damping_default) > 10.0_dp*EPSILON(1.0_dp)) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/DAMPING." + WRITE (UNIT=iounit, FMT="(T2,A,F8.4,A,F8.4,A)") & + "The native reference run uses tblite's library default ", tblite_mixer_damping_default, & + ", while CP2K uses DAMPING ", dft_control%qs_control%xtb_control%tblite_mixer_damping, "." + IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/DAMPING") + END IF + IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_omega0 - & + tblite_mixer_omega0_default) > 10.0_dp*EPSILON(1.0_dp)) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/OMEGA0." + WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") & + "The native reference run uses tblite's library default ", tblite_mixer_omega0_default, & + ", while CP2K uses OMEGA0 ", dft_control%qs_control%xtb_control%tblite_mixer_omega0, "." + IF (ref%stop_on_error) CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/OMEGA0") + END IF + IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_min_weight - & + tblite_mixer_min_weight_default) > 10.0_dp*EPSILON(1.0_dp)) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MIN_WEIGHT." + WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") & + "The native reference run uses tblite's library default ", tblite_mixer_min_weight_default, & + ", while CP2K uses MIN_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_min_weight, "." + IF (ref%stop_on_error) & + CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MIN_WEIGHT") + END IF + IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_max_weight - & + tblite_mixer_max_weight_default) > 10.0_dp*EPSILON(1.0_dp)) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/MAX_WEIGHT." + WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") & + "The native reference run uses tblite's library default ", tblite_mixer_max_weight_default, & + ", while CP2K uses MAX_WEIGHT ", dft_control%qs_control%xtb_control%tblite_mixer_max_weight, "." + IF (ref%stop_on_error) & + CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/MAX_WEIGHT") + END IF + IF (ABS(dft_control%qs_control%xtb_control%tblite_mixer_weight_factor - & + tblite_mixer_weight_factor_default) > 10.0_dp*EPSILON(1.0_dp)) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI cannot reproduce XTB/TBLITE_MIXER/WEIGHT_FACTOR." + WRITE (UNIT=iounit, FMT="(T2,A,ES12.4,A,ES12.4,A)") & + "The native reference run uses tblite's library default ", tblite_mixer_weight_factor_default, & + ", while CP2K uses WEIGHT_FACTOR ", & + dft_control%qs_control%xtb_control%tblite_mixer_weight_factor, "." + IF (ref%stop_on_error) & + CPABORT("tblite reference CLI cannot reproduce TBLITE_MIXER/WEIGHT_FACTOR") + END IF + etemp = 300.0_dp + IF (ASSOCIATED(scf_control) .AND. ASSOCIATED(scf_control%smear)) THEN + IF (scf_control%smear%do_smear) THEN + etemp = cp_unit_from_cp2k(scf_control%smear%electronic_temperature, "K") + IF (scf_control%smear%method /= smear_fermi_dirac) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A,A,A)") & + "WARNING: tblite reference CLI cannot reproduce CP2K smearing method ", & + TRIM(tb_reference_smear_method_name(scf_control%smear%method)), "." + WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") & + "The native reference run uses Fermi-Dirac electronic temperature ", etemp, " K instead." + END IF + END IF + END IF + WRITE (etemp_str, "(ES16.8)") etemp + etemp_guess = 0.0_dp + etemp_guess_str = "" + IF (ref%electronic_temperature_guess > 0.0_dp) THEN + etemp_guess = cp_unit_from_cp2k(ref%electronic_temperature_guess, "K") + WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess + etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str)) + END IF + param_str = "" + IF (LEN_TRIM(dft_control%qs_control%xtb_control%tblite_param_file) > 0) & + param_str = " --param "//TRIM(tb_shell_quote(dft_control%qs_control%xtb_control%tblite_param_file)) spinpol_str = "" IF (dft_control%lsd) spinpol_str = " --spin-polarized" + post_processing_str = "" + IF (LEN_TRIM(ref%post_processing) > 0) & + post_processing_str = " --post-processing "//TRIM(tb_shell_quote(ref%post_processing)) + post_processing_output_str = "" + IF (LEN_TRIM(post_processing_output_file) > 0) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A)") & + "WARNING: tblite reference CLI POST_PROCESSING_OUTPUT was requested explicitly." + WRITE (UNIT=iounit, FMT="(T2,A)") & + "Some tblite 0.5.0 command-line builds document --post-processing-output but do not parse it." + post_processing_output_str = " --post-processing-output "// & + TRIM(tb_shell_quote(post_processing_output_file)) + END IF + restart_str = " --no-restart" + IF (LEN_TRIM(ref%restart_file) > 0) & + restart_str = " --restart "//TRIM(tb_shell_quote(ref%restart_file)) CALL tb_write_reference_gen(qs_env, TRIM(gen_file)) command = TRIM(tb_shell_quote(ref%program_name))//" run --method "//TRIM(method)// & + TRIM(param_str)// & TRIM(spinpol_str)// & " --charge "//TRIM(ADJUSTL(charge_str))// & " --spin "//TRIM(ADJUSTL(spin_str))// & " --acc "//TRIM(ADJUSTL(acc_str))// & + " --guess "//TRIM(guess)// & + " --solver "//TRIM(solver)// & " --iterations "//TRIM(ADJUSTL(iter_str))// & - " --etemp 0.0"// & - " --no-restart --input gen --grad "//TRIM(tb_shell_quote(grad_file))// & + " --etemp "//TRIM(ADJUSTL(etemp_str))// & + TRIM(etemp_guess_str)// & + TRIM(efield_str)// & + TRIM(solvation_str)// & + TRIM(post_processing_str)// & + TRIM(post_processing_output_str)// & + TRIM(restart_str)// & + TRIM(verbosity_str)//" --input "//TRIM(tb_shell_quote(ref%input_format))// & + " --grad "//TRIM(tb_shell_quote(grad_file))// & " --json "//TRIM(tb_shell_quote(json_file))//" "//TRIM(tb_shell_quote(gen_file))// & " > "//TRIM(tb_shell_quote(log_file))//" 2>&1" @@ -4312,7 +4701,7 @@ CONTAINS WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command) WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat IF (ref%stop_on_error) CPABORT("tblite reference CLI command failed") - CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file) + CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file) CALL timestop(handle) RETURN END IF @@ -4324,6 +4713,22 @@ CONTAINS WRITE (UNIT=iounit, FMT="(/,T2,A)") "tblite reference CLI check" WRITE (UNIT=iounit, FMT="(T2,A,A)") "Executable: ", TRIM(ref%program_name) WRITE (UNIT=iounit, FMT="(T2,A,A)") "Method: ", TRIM(method) + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Guess: ", TRIM(guess) + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Solver: ", TRIM(solver) + WRITE (UNIT=iounit, FMT="(T2,A,L1)") "Spin-pol.: ", dft_control%lsd + IF (ref%efield_active) & + WRITE (UNIT=iounit, FMT="(T2,A,3ES16.8,A)") "Efield: ", ref%efield, " V/Angstrom" + IF (ref%solvation_active) & + WRITE (UNIT=iounit, FMT="(T2,A,A,1X,A)") "Solvation: ", TRIM(solvation_model_name), & + TRIM(ref%solvation_solvent) + IF (LEN_TRIM(ref%post_processing) > 0) & + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Post proc.: ", TRIM(ref%post_processing) + IF (LEN_TRIM(post_processing_output_file) > 0) & + WRITE (UNIT=iounit, FMT="(T2,A,A)") "PP output: ", TRIM(post_processing_output_file) + IF (ref%electronic_temperature_guess > 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T2,A,F12.3,A)") "Guess etemp:", etemp_guess, " K" + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Grad file: ", TRIM(grad_file) + WRITE (UNIT=iounit, FMT="(T2,A,A)") "JSON file: ", TRIM(json_file) WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file) too_large = .FALSE. @@ -4372,7 +4777,9 @@ CONTAINS IF (ref%stop_on_error) CPABORT("tblite reference CLI deviation exceeded ERROR_LIMIT") END IF - CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file) + CALL tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit) + + CALL tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file) DEALLOCATE (cli_gradient, cli_virial) CALL timestop(handle) @@ -4401,6 +4808,234 @@ CONTAINS END FUNCTION tb_reference_method_name +! ************************************************************************************************** +!> \brief Map CP2K tblite reference CLI guess id to native tblite CLI guess name. +!> \param guess_id ... +!> \return ... +! ************************************************************************************************** + FUNCTION tb_reference_guess_name(guess_id) RESULT(guess) + INTEGER, INTENT(IN) :: guess_id + CHARACTER(LEN=8) :: guess + + SELECT CASE (guess_id) + CASE (tblite_guess_sad) + guess = "sad" + CASE (tblite_guess_eeq) + guess = "eeq" + CASE (tblite_guess_ceh) + guess = "ceh" + CASE DEFAULT + CPABORT("Unknown tblite reference CLI guess") + END SELECT + + END FUNCTION tb_reference_guess_name + +! ************************************************************************************************** +!> \brief Map CP2K tblite solver id to native tblite CLI solver name. +!> \param solver_id ... +!> \return ... +! ************************************************************************************************** + FUNCTION tb_reference_solver_name(solver_id) RESULT(solver) + INTEGER, INTENT(IN) :: solver_id + CHARACTER(LEN=8) :: solver + + SELECT CASE (solver_id) + CASE (tblite_solver_gvd) + solver = "gvd" + CASE (tblite_solver_gvr) + solver = "gvr" + CASE DEFAULT + CPABORT("Unknown tblite reference CLI solver") + END SELECT + + END FUNCTION tb_reference_solver_name + +! ************************************************************************************************** +!> \brief Map CP2K smearing method id to a diagnostic label. +!> \param method_id ... +!> \return ... +! ************************************************************************************************** + FUNCTION tb_reference_smear_method_name(method_id) RESULT(method) + INTEGER, INTENT(IN) :: method_id + CHARACTER(LEN=24) :: method + + SELECT CASE (method_id) + CASE (smear_fermi_dirac) + method = "FERMI_DIRAC" + CASE (smear_energy_window) + method = "ENERGY_WINDOW" + CASE (smear_list) + method = "LIST" + CASE (smear_gaussian) + method = "GAUSSIAN" + CASE (smear_mp) + method = "METHFESSEL_PAXTON" + CASE (smear_mv) + method = "MARZARI_VANDERBILT" + CASE DEFAULT + method = "UNKNOWN" + END SELECT + + END FUNCTION tb_reference_smear_method_name + +! ************************************************************************************************** +!> \brief Run optional native tblite auxiliary subcommands. +!> \param ref ... +!> \param dft_control ... +!> \param gen_file ... +!> \param file_base ... +!> \param verbosity_str ... +!> \param iounit ... +! ************************************************************************************************** + SUBROUTINE tb_reference_cli_aux_commands(ref, dft_control, gen_file, file_base, verbosity_str, iounit) + + TYPE(xtb_reference_cli_type), INTENT(IN) :: ref + TYPE(dft_control_type), INTENT(IN) :: dft_control + CHARACTER(LEN=*), INTENT(IN) :: gen_file, file_base, verbosity_str + INTEGER, INTENT(IN) :: iounit + + CHARACTER(LEN=32) :: charge_str, efield_x_str, efield_y_str, & + efield_z_str, etemp_guess_val_str, & + spin_str + CHARACTER(LEN=4*default_path_length+16) :: copy_str, dry_run_str, efield_str, & + etemp_guess_str, grad_str, json_str, & + method_str, output_str + CHARACTER(LEN=8*default_path_length) :: command + CHARACTER(LEN=default_path_length) :: guess_input, log_file + INTEGER :: spin + REAL(KIND=dp) :: etemp_guess + + WRITE (charge_str, "(I0)") dft_control%charge + spin = MAX(0, dft_control%multiplicity - 1) + WRITE (spin_str, "(I0)") spin + + IF (ref%guess_cli%enabled) THEN + guess_input = ref%guess_cli%input_file + IF (LEN_TRIM(guess_input) == 0) guess_input = gen_file + etemp_guess_str = "" + IF (ref%guess_cli%electronic_temperature_guess > 0.0_dp) THEN + etemp_guess = cp_unit_from_cp2k(ref%guess_cli%electronic_temperature_guess, "K") + WRITE (etemp_guess_val_str, "(ES16.8)") etemp_guess + etemp_guess_str = " --etemp-guess "//TRIM(ADJUSTL(etemp_guess_val_str)) + END IF + efield_str = "" + IF (ref%guess_cli%efield_active) THEN + WRITE (efield_x_str, "(ES16.8)") ref%guess_cli%efield(1) + WRITE (efield_y_str, "(ES16.8)") ref%guess_cli%efield(2) + WRITE (efield_z_str, "(ES16.8)") ref%guess_cli%efield(3) + efield_str = " --efield "//TRIM(ADJUSTL(efield_x_str))//","// & + TRIM(ADJUSTL(efield_y_str))//","//TRIM(ADJUSTL(efield_z_str)) + END IF + grad_str = "" + IF (ref%guess_cli%grad) grad_str = " --grad" + json_str = "" + IF (LEN_TRIM(ref%guess_cli%json_file) > 0) & + json_str = " --json "//TRIM(tb_shell_quote(ref%guess_cli%json_file)) + log_file = TRIM(file_base)//".guess.log" + command = TRIM(tb_shell_quote(ref%program_name))// & + " guess --charge "//TRIM(ADJUSTL(charge_str))// & + " --spin "//TRIM(ADJUSTL(spin_str))// & + " --method "//TRIM(tb_reference_guess_name(ref%guess_cli%method))// & + " --solver "//TRIM(tb_reference_solver_name(ref%guess_cli%solver))// & + TRIM(etemp_guess_str)// & + TRIM(efield_str)// & + TRIM(grad_str)// & + TRIM(json_str)// & + TRIM(verbosity_str)//" --input "//TRIM(tb_shell_quote(ref%guess_cli%input_format))// & + " "//TRIM(tb_shell_quote(guess_input))// & + " > "//TRIM(tb_shell_quote(log_file))//" 2>&1" + CALL tb_reference_cli_execute(ref, "guess", command, log_file, iounit) + IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%guess_cli%json_file) > 0) & + CALL tb_delete_file(ref%guess_cli%json_file) + END IF + + IF (ref%param_cli%enabled) THEN + method_str = "" + IF (ref%param_cli%method_explicit .OR. LEN_TRIM(ref%param_cli%input_file) == 0) & + method_str = " --method "// & + TRIM(tb_reference_method_name(MERGE(ref%param_cli%method, & + dft_control%qs_control%xtb_control%tblite_method, & + ref%param_cli%method_explicit))) + output_str = "" + IF (LEN_TRIM(ref%param_cli%output_file) > 0) & + output_str = " --output "//TRIM(tb_shell_quote(ref%param_cli%output_file)) + log_file = TRIM(file_base)//".param.log" + command = TRIM(tb_shell_quote(ref%program_name))//" param"// & + TRIM(method_str)// & + TRIM(output_str) + IF (LEN_TRIM(ref%param_cli%input_file) > 0) & + command = TRIM(command)//" "//TRIM(tb_shell_quote(ref%param_cli%input_file)) + command = TRIM(command)//" > "//TRIM(tb_shell_quote(log_file))//" 2>&1" + CALL tb_reference_cli_execute(ref, "param", command, log_file, iounit) + IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%param_cli%output_file) > 0) & + CALL tb_delete_file(ref%param_cli%output_file) + END IF + + IF (ref%fit_cli%enabled) THEN + dry_run_str = "" + IF (ref%fit_cli%dry_run) dry_run_str = " --dry-run" + copy_str = "" + IF (LEN_TRIM(ref%fit_cli%copy_file) > 0) & + copy_str = " --copy "//TRIM(tb_shell_quote(ref%fit_cli%copy_file)) + log_file = TRIM(file_base)//".fit.log" + command = TRIM(tb_shell_quote(ref%program_name))//" fit"// & + TRIM(dry_run_str)// & + TRIM(copy_str)// & + TRIM(verbosity_str)//" "//TRIM(tb_shell_quote(ref%fit_cli%param_file))// & + " "//TRIM(tb_shell_quote(ref%fit_cli%input_file))// & + " > "//TRIM(tb_shell_quote(log_file))//" 2>&1" + CALL tb_reference_cli_execute(ref, "fit", command, log_file, iounit) + IF (.NOT. ref%keep_files .AND. LEN_TRIM(ref%fit_cli%copy_file) > 0) & + CALL tb_delete_file(ref%fit_cli%copy_file) + END IF + + IF (ref%tagdiff_cli%enabled) THEN + method_str = "" + IF (ref%tagdiff_cli%fit) method_str = " --fit" + log_file = TRIM(file_base)//".tagdiff.log" + command = TRIM(tb_shell_quote(ref%program_name))//" tagdiff"// & + TRIM(method_str)//" "//TRIM(tb_shell_quote(ref%tagdiff_cli%actual_file))// & + " "//TRIM(tb_shell_quote(ref%tagdiff_cli%reference_file))// & + " > "//TRIM(tb_shell_quote(log_file))//" 2>&1" + CALL tb_reference_cli_execute(ref, "tagdiff", command, log_file, iounit) + END IF + + END SUBROUTINE tb_reference_cli_aux_commands + +! ************************************************************************************************** +!> \brief Execute one native tblite REFERENCE_CLI auxiliary command. +!> \param ref ... +!> \param label ... +!> \param command ... +!> \param log_file ... +!> \param iounit ... +! ************************************************************************************************** + SUBROUTINE tb_reference_cli_execute(ref, label, command, log_file, iounit) + + TYPE(xtb_reference_cli_type), INTENT(IN) :: ref + CHARACTER(LEN=*), INTENT(IN) :: label, command, log_file + INTEGER, INTENT(IN) :: iounit + + INTEGER :: cmdstat, exitstat + + cmdstat = 0 + exitstat = 0 + CALL execute_command_line(TRIM(command), exitstat=exitstat, cmdstat=cmdstat) + IF (cmdstat /= 0 .OR. exitstat /= 0) THEN + WRITE (UNIT=iounit, FMT="(/,T2,A,A)") "tblite reference CLI auxiliary command failed: ", & + TRIM(label) + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Command: ", TRIM(command) + WRITE (UNIT=iounit, FMT="(T2,A,I0,T32,A,I0)") "cmdstat:", cmdstat, "exitstat:", exitstat + IF (ref%stop_on_error) CPABORT("tblite reference CLI auxiliary command failed") + ELSE + WRITE (UNIT=iounit, FMT="(/,T2,A,A)") "tblite reference CLI auxiliary command completed: ", & + TRIM(label) + WRITE (UNIT=iounit, FMT="(T2,A,A)") "Log file: ", TRIM(log_file) + END IF + IF (.NOT. ref%keep_files) CALL tb_delete_file(log_file) + + END SUBROUTINE tb_reference_cli_execute + ! ************************************************************************************************** !> \brief Join directory and filename. !> \param directory ... @@ -4670,17 +5305,21 @@ CONTAINS !> \param grad_file ... !> \param json_file ... !> \param log_file ... +!> \param post_processing_output_file ... ! ************************************************************************************************** - SUBROUTINE tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file) + SUBROUTINE tb_reference_cleanup(ref, gen_file, grad_file, json_file, log_file, post_processing_output_file) TYPE(xtb_reference_cli_type), INTENT(IN) :: ref - CHARACTER(LEN=*), INTENT(IN) :: gen_file, grad_file, json_file, log_file + CHARACTER(LEN=*), INTENT(IN) :: gen_file, grad_file, json_file, & + log_file, post_processing_output_file IF (ref%keep_files) RETURN CALL tb_delete_file(gen_file) CALL tb_delete_file(grad_file) CALL tb_delete_file(json_file) CALL tb_delete_file(log_file) + IF (LEN_TRIM(post_processing_output_file) > 0) & + CALL tb_delete_file(post_processing_output_file) END SUBROUTINE tb_reference_cleanup diff --git a/src/tblite_scc_mixer.F b/src/tblite_scc_mixer.F new file mode 100644 index 0000000000..237f33576f --- /dev/null +++ b/src/tblite_scc_mixer.F @@ -0,0 +1,283 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief CP2K-side tblite-compatible SCC Broyden mixer. +! ************************************************************************************************** +MODULE tblite_scc_mixer +#if defined(__TBLITE) + USE mctc_env, ONLY: error_type, & + fatal_error, & + wp + USE tblite_lapack, ONLY: getrf, & + getrs + USE tblite_scf_mixer_type, ONLY: mixer_type +#endif + + IMPLICIT NONE + + PRIVATE + +#if defined(__TBLITE) + PUBLIC :: new_cp2k_tblite_mixer + + TYPE, EXTENDS(mixer_type) :: cp2k_tblite_broyden_mixer_type + INTEGER :: idif = 0 + INTEGER :: iget = 0 + INTEGER :: iset = 0 + INTEGER :: iter = 0 + INTEGER :: memory = 0 + INTEGER :: ndim = 0 + REAL(wp) :: damp = 0.0_wp + REAL(wp) :: max_weight = 0.0_wp + REAL(wp) :: min_weight = 0.0_wp + REAL(wp) :: omega0 = 0.0_wp + REAL(wp) :: weight_factor = 0.0_wp + REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: a, df, u + REAL(wp), ALLOCATABLE, DIMENSION(:) :: dq, dqlast, omega, q_in, qlast_in + CONTAINS + PROCEDURE :: diff_1d => cp2k_tblite_mixer_diff_1d + PROCEDURE :: get_1d => cp2k_tblite_mixer_get_1d + PROCEDURE :: get_error => cp2k_tblite_mixer_get_error + PROCEDURE :: next => cp2k_tblite_mixer_next + PROCEDURE :: set_1d => cp2k_tblite_mixer_set_1d + END TYPE cp2k_tblite_broyden_mixer_type +#endif + +CONTAINS + +#if defined(__TBLITE) +! ************************************************************************************************** +!> \brief Create a CP2K-side tblite-compatible Broyden mixer. +!> \param mixer mixer instance +!> \param memory Broyden history length +!> \param ndim number of SCC variables +!> \param damping first-step damping +!> \param omega0 Broyden regularization weight +!> \param min_weight minimum dynamic Broyden weight +!> \param max_weight maximum dynamic Broyden weight +!> \param weight_factor residual-to-weight scaling factor +! ************************************************************************************************** + SUBROUTINE new_cp2k_tblite_mixer(mixer, memory, ndim, damping, omega0, min_weight, max_weight, & + weight_factor) + CLASS(mixer_type), ALLOCATABLE, INTENT(OUT) :: mixer + INTEGER, INTENT(IN) :: memory, ndim + REAL(wp), INTENT(IN) :: damping, omega0, min_weight, max_weight, & + weight_factor + + TYPE(cp2k_tblite_broyden_mixer_type), ALLOCATABLE :: broyden + + ALLOCATE (broyden) + broyden%ndim = ndim + broyden%memory = memory + broyden%damp = damping + broyden%omega0 = omega0 + broyden%min_weight = min_weight + broyden%max_weight = max_weight + broyden%weight_factor = weight_factor + ALLOCATE (broyden%a(memory, memory)) + ALLOCATE (broyden%df(ndim, memory)) + ALLOCATE (broyden%u(ndim, memory)) + ALLOCATE (broyden%dq(ndim)) + ALLOCATE (broyden%dqlast(ndim)) + ALLOCATE (broyden%omega(memory)) + ALLOCATE (broyden%q_in(ndim)) + ALLOCATE (broyden%qlast_in(ndim)) + CALL MOVE_ALLOC(broyden, mixer) + + END SUBROUTINE new_cp2k_tblite_mixer + +! ************************************************************************************************** +!> \brief Set input SCC variables. +!> \param self mixer +!> \param qvec SCC-variable vector +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_mixer_set_1d(self, qvec) + CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self + REAL(wp), DIMENSION(:), INTENT(IN) :: qvec + + self%q_in(self%iset + 1:self%iset + SIZE(qvec)) = qvec + self%iset = self%iset + SIZE(qvec) + + END SUBROUTINE cp2k_tblite_mixer_set_1d + +! ************************************************************************************************** +!> \brief Set SCC residual. +!> \param self mixer +!> \param qvec output SCC-variable vector before mixing +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_mixer_diff_1d(self, qvec) + CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self + REAL(wp), DIMENSION(:), INTENT(IN) :: qvec + + self%dq(self%idif + 1:self%idif + SIZE(qvec)) = & + qvec - self%q_in(self%idif + 1:self%idif + SIZE(qvec)) + self%idif = self%idif + SIZE(qvec) + + END SUBROUTINE cp2k_tblite_mixer_diff_1d + +! ************************************************************************************************** +!> \brief Apply the next Broyden update. +!> \param self mixer +!> \param error error handling +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_mixer_next(self, error) + CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self + TYPE(error_type), ALLOCATABLE, INTENT(OUT) :: error + + INTEGER :: info + + self%iset = 0 + self%idif = 0 + self%iget = 0 + self%iter = self%iter + 1 + CALL cp2k_tblite_broyden(self%ndim, self%q_in, self%qlast_in, self%dq, self%dqlast, & + self%iter, self%memory, self%damp, self%omega0, & + self%min_weight, self%max_weight, self%weight_factor, & + self%omega, self%df, self%u, self%a, info) + IF (info /= 0) CALL fatal_error(error, "Broyden mixing failed to obtain next iteration") + + END SUBROUTINE cp2k_tblite_mixer_next + +! ************************************************************************************************** +!> \brief Get mixed SCC variables. +!> \param self mixer +!> \param qvec mixed SCC-variable vector +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_mixer_get_1d(self, qvec) + CLASS(cp2k_tblite_broyden_mixer_type), INTENT(INOUT) :: self + REAL(wp), DIMENSION(:), INTENT(OUT) :: qvec + + qvec(:) = self%q_in(self%iget + 1:self%iget + SIZE(qvec)) + self%iget = self%iget + SIZE(qvec) + + END SUBROUTINE cp2k_tblite_mixer_get_1d + +! ************************************************************************************************** +!> \brief Get RMS SCC residual. +!> \param self mixer +!> \return residual +! ************************************************************************************************** + PURE FUNCTION cp2k_tblite_mixer_get_error(self) RESULT(error) + CLASS(cp2k_tblite_broyden_mixer_type), INTENT(IN) :: self + REAL(wp) :: error + + error = SQRT(SUM(self%dq**2)/SIZE(self%dq)) + + END FUNCTION cp2k_tblite_mixer_get_error + +! ************************************************************************************************** +!> \brief Modified Broyden update matching tblite's default algorithm, with explicit constants. +!> \param n ... +!> \param q ... +!> \param qlast ... +!> \param dq ... +!> \param dqlast ... +!> \param iter ... +!> \param memory ... +!> \param alpha ... +!> \param omega0 ... +!> \param minw ... +!> \param maxw ... +!> \param wfac ... +!> \param omega ... +!> \param df ... +!> \param u ... +!> \param a ... +!> \param info ... +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_broyden(n, q, qlast, dq, dqlast, iter, memory, alpha, omega0, minw, & + maxw, wfac, omega, df, u, a, info) + INTEGER, INTENT(IN) :: n + REAL(wp), DIMENSION(n), INTENT(INOUT) :: q, qlast + REAL(wp), DIMENSION(n), INTENT(IN) :: dq + REAL(wp), DIMENSION(n), INTENT(INOUT) :: dqlast + INTEGER, INTENT(IN) :: iter, memory + REAL(wp), INTENT(IN) :: alpha, omega0, minw, maxw, wfac + REAL(wp), DIMENSION(memory), INTENT(INOUT) :: omega + REAL(wp), DIMENSION(n, memory), INTENT(INOUT) :: df, u + REAL(wp), DIMENSION(memory, memory), INTENT(INOUT) :: a + INTEGER, INTENT(OUT) :: info + + INTEGER :: i, it1, itn, j + REAL(wp) :: inv + REAL(wp), ALLOCATABLE, DIMENSION(:, :) :: beta, c + + info = 0 + itn = iter - 1 + it1 = MOD(itn - 1, memory) + 1 + + IF (iter == 1) THEN + dqlast(:) = dq + qlast(:) = q + q(:) = q + alpha*dq + RETURN + END IF + + ALLOCATE (beta(MIN(memory, itn), MIN(memory, itn)), c(MIN(memory, itn), 1)) + + omega(it1) = SQRT(DOT_PRODUCT(dq, dq)) + IF (omega(it1) > (wfac/maxw)) THEN + omega(it1) = wfac/omega(it1) + ELSE + omega(it1) = maxw + END IF + IF (omega(it1) < minw) omega(it1) = minw + + df(:, it1) = dq - dqlast + inv = MAX(SQRT(DOT_PRODUCT(df(:, it1), df(:, it1))), EPSILON(1.0_wp)) + inv = 1.0_wp/inv + df(:, it1) = inv*df(:, it1) + + DO j = MAX(1, itn - memory + 1), itn + i = MOD(j - 1, memory) + 1 + a(i, it1) = DOT_PRODUCT(df(:, i), df(:, it1)) + a(it1, i) = a(i, it1) + c(i, 1) = omega(i)*DOT_PRODUCT(df(:, i), dq) + END DO + + DO j = MAX(1, itn - memory + 1), itn + i = MOD(j - 1, memory) + 1 + beta(:it1, i) = omega(:it1)*omega(i)*a(:it1, i) + beta(i, i) = beta(i, i) + omega0*omega0 + END DO + + CALL cp2k_tblite_lineq(beta, c, info) + IF (info /= 0) RETURN + + u(:, it1) = alpha*df(:, it1) + inv*(q - qlast) + dqlast(:) = dq + qlast(:) = q + q(:) = q + alpha*dq + + DO j = MAX(1, itn - memory + 1), itn + i = MOD(j - 1, memory) + 1 + q(:) = q - omega(i)*c(i, 1)*u(:, i) + END DO + + END SUBROUTINE cp2k_tblite_broyden + +! ************************************************************************************************** +!> \brief Solve the small Broyden linear system. +!> \param a ... +!> \param c ... +!> \param info ... +! ************************************************************************************************** + SUBROUTINE cp2k_tblite_lineq(a, c, info) + REAL(wp), DIMENSION(:, :), INTENT(INOUT) :: a, c + INTEGER, INTENT(OUT) :: info + + INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv + + ALLOCATE (ipiv(SIZE(a, 1))) + CALL getrf(a, ipiv, info) + IF (info == 0) CALL getrs(a, c, ipiv, info, trans="t") + + END SUBROUTINE cp2k_tblite_lineq +#endif + +END MODULE tblite_scc_mixer diff --git a/src/tblite_types.F b/src/tblite_types.F index 7be9c26c1a..5fb1e505c7 100644 --- a/src/tblite_types.F +++ b/src/tblite_types.F @@ -20,10 +20,18 @@ MODULE tblite_types USE tblite_scf_mixer, ONLY: mixer_type USE tblite_wavefunction_type, ONLY: wavefunction_type USE tblite_scf_potential, ONLY: potential_type + USE tblite_param, ONLY: param_record #endif USE cp_dbcsr_api, ONLY: dbcsr_p_type USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set + USE input_constants, ONLY: tblite_mixer_damping_default, & + tblite_mixer_iterations_default, & + tblite_mixer_max_weight_default, & + tblite_mixer_min_weight_default, & + tblite_mixer_omega0_default, & + tblite_mixer_weight_factor_default, & + tblite_solver_gvd USE kinds, ONLY: dp #include "./base/base_uses.f90" @@ -62,6 +70,13 @@ MODULE tblite_types TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quadbra => Null(), & quadket => Null() TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp_ref => Null() + INTEGER :: mixer_memory = tblite_mixer_iterations_default + INTEGER :: mixer_solver = tblite_solver_gvd + REAL(KIND=dp) :: mixer_damping = tblite_mixer_damping_default + REAL(KIND=dp) :: mixer_max_weight = tblite_mixer_max_weight_default + REAL(KIND=dp) :: mixer_min_weight = tblite_mixer_min_weight_default + REAL(KIND=dp) :: mixer_omega0 = tblite_mixer_omega0_default + REAL(KIND=dp) :: mixer_weight_factor = tblite_mixer_weight_factor_default #if defined(__TBLITE) TYPE(structure_type) :: mol = structure_type() @@ -72,6 +87,7 @@ MODULE tblite_types TYPE(wavefunction_type) :: wfn = wavefunction_type() TYPE(potential_type) :: pot = potential_type() CLASS(mixer_type), ALLOCATABLE :: mixer + TYPE(param_record), ALLOCATABLE :: param #endif END TYPE tblite_type @@ -129,6 +145,7 @@ CONTAINS CALL dbcsr_deallocate_matrix_set(tb_tblite%rho_ao_kp_ref) #if defined(__TBLITE) IF (ALLOCATED(tb_tblite%mixer)) DEALLOCATE (tb_tblite%mixer) + IF (ALLOCATED(tb_tblite%param)) DEALLOCATE (tb_tblite%param) #endif DEALLOCATE (tb_tblite) diff --git a/src/xtb_ks_matrix.F b/src/xtb_ks_matrix.F index beb30d6f30..a139419cc7 100644 --- a/src/xtb_ks_matrix.F +++ b/src/xtb_ks_matrix.F @@ -357,8 +357,17 @@ CONTAINS CALL charge_mixing(scf_env%mixing_method, scf_env%mixing_store, & charges, para_env, scf_env%iter_count, & scc_mixer=dft_control%qs_control%xtb_control%tblite_scc_mixer, & + tblite_mixer_iterations= & + dft_control%qs_control%xtb_control%tblite_mixer_iterations, & tblite_mixer_damping=dft_control%qs_control%xtb_control%tblite_mixer_damping, & - tblite_mixer_memory=qs_env%scf_control%max_scf) + tblite_mixer_memory=dft_control%qs_control%xtb_control%tblite_mixer_memory, & + tblite_mixer_omega0=dft_control%qs_control%xtb_control%tblite_mixer_omega0, & + tblite_mixer_min_weight= & + dft_control%qs_control%xtb_control%tblite_mixer_min_weight, & + tblite_mixer_max_weight= & + dft_control%qs_control%xtb_control%tblite_mixer_max_weight, & + tblite_mixer_weight_factor= & + dft_control%qs_control%xtb_control%tblite_mixer_weight_factor) END IF DO iatom = 1, natom diff --git a/tests/DFTB/regtest-scc/ch2o_tblite_mixer.inp b/tests/DFTB/regtest-scc/ch2o_tblite_mixer.inp index 7800e96e3a..c934b03e11 100644 --- a/tests/DFTB/regtest-scc/ch2o_tblite_mixer.inp +++ b/tests/DFTB/regtest-scc/ch2o_tblite_mixer.inp @@ -23,11 +23,16 @@ ORTHOGONAL_BASIS F SCC_MIXER TBLITE SELF_CONSISTENT T - TBLITE_MIXER_DAMPING 0.005 &PARAMETER PARAM_FILE_NAME scc_parameter PARAM_FILE_PATH DFTB/scc &END PARAMETER + &TBLITE_MIXER + DAMPING 0.005 + ITERATIONS 500 + MEMORY 500 + SOLVER GVD + &END TBLITE_MIXER &END DFTB &END QS &SCF diff --git a/tests/xTB/regtest-3-spglib/TEST_FILES.toml b/tests/xTB/regtest-3-spglib/TEST_FILES.toml index d32ec2a963..07d43412a9 100644 --- a/tests/xTB/regtest-3-spglib/TEST_FILES.toml +++ b/tests/xTB/regtest-3-spglib/TEST_FILES.toml @@ -1,6 +1,6 @@ # runs are executed in the same order as in this file -"si_kp_kpsym.inp" = [{matcher="E_total", tol=1.0E-11, ref=-14.74208554207339}, +"si_kp_kpsym.inp" = [{matcher="E_total", tol=1.0E-11, ref=-14.74208554419625}, {matcher="N_special_kpoints", tol=0.0, ref=4}] -"si_kp_spglib_backend.inp" = [{matcher="E_total", tol=1.0E-11, ref=-14.74208554207337}, +"si_kp_spglib_backend.inp" = [{matcher="E_total", tol=1.0E-11, ref=-14.74208554419625}, {matcher="N_special_kpoints", tol=0.0, ref=4}] #EOF diff --git a/tests/xTB/regtest-3/TEST_FILES.toml b/tests/xTB/regtest-3/TEST_FILES.toml index 867ad404bb..dda2961d8b 100644 --- a/tests/xTB/regtest-3/TEST_FILES.toml +++ b/tests/xTB/regtest-3/TEST_FILES.toml @@ -6,8 +6,8 @@ "NdF3.inp" = [{matcher="E_total", tol=1.0E-12, ref=-16.30904020314352}] "h2o_rtp.inp" = [{matcher="E_total", tol=1.0E-12, ref=-5.76539102399686}] "h2o_emd.inp" = [{matcher="E_total", tol=1.0E-12, ref=-5.76569266115801}] -"si8_wan.inp" = [{matcher="E_total", tol=1.0E-12, ref=-14.36588818966667}] -"si_kp.inp" = [{matcher="E_total", tol=1.0E-12, ref=-14.73033123294306}] +"si8_wan.inp" = [{matcher="E_total", tol=1.0E-12, ref=-14.36588825040988}] +"si_kp.inp" = [{matcher="E_total", tol=1.0E-12, ref=-14.73033107961160}] "tmol.inp" = [{matcher="E_total", tol=1.0E-12, ref=-41.90845660778482}] "ch2o.inp" = [{matcher="E_total", tol=1.0E-12, ref=-7.84456570305608}] "ch2o_print.inp" = [{matcher="E_total", tol=1.0E-12, ref=-7.84456570305608}] diff --git a/tests/xTB/regtest-gfn1-d/TEST_FILES.toml b/tests/xTB/regtest-gfn1-d/TEST_FILES.toml index a687c06366..3076729742 100644 --- a/tests/xTB/regtest-gfn1-d/TEST_FILES.toml +++ b/tests/xTB/regtest-gfn1-d/TEST_FILES.toml @@ -26,7 +26,7 @@ {matcher="DEBUG_periodic_stress_sum", tol=1.0E-7, ref=0.0000000}, {matcher="DEBUG_force_sum", tol=1.0E-6, ref=0.0000000}, {matcher="N_special_kpoints", tol=0.0, ref=2}] -"si_kp_kpsym_debug_stress.inp" = [{matcher="DEBUG_stress_sum", tol=1.0E-7, ref=0.0000000}, +"si_kp_kpsym_debug_stress.inp" = [{matcher="DEBUG_stress_sum", tol=1.0E-6, ref=0.0000000}, {matcher="N_special_kpoints", tol=0.0, ref=32}] # #ch2o_t09.inp 0 diff --git a/tests/xTB/regtest-tblite-gfn1-grad/CH2O_gfn1_uks_force.inp b/tests/xTB/regtest-tblite-gfn1-grad/CH2O_gfn1_uks_force.inp index a93ac10fa5..1c4667c0e9 100644 --- a/tests/xTB/regtest-tblite-gfn1-grad/CH2O_gfn1_uks_force.inp +++ b/tests/xTB/regtest-tblite-gfn1-grad/CH2O_gfn1_uks_force.inp @@ -21,9 +21,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER TBLITE &TBLITE METHOD GFN1 - SCC_MIXER TBLITE &END TBLITE &END XTB &END QS diff --git a/tests/xTB/regtest-tblite-gfn1-grad/O2_gfn1_uks_force.inp b/tests/xTB/regtest-tblite-gfn1-grad/O2_gfn1_uks_force.inp index b64bcccc1a..8602d5d03e 100644 --- a/tests/xTB/regtest-tblite-gfn1-grad/O2_gfn1_uks_force.inp +++ b/tests/xTB/regtest-tblite-gfn1-grad/O2_gfn1_uks_force.inp @@ -20,9 +20,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER TBLITE &TBLITE METHOD GFN1 - SCC_MIXER TBLITE &END TBLITE &END XTB &END QS diff --git a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_force_cp2k_mixer.inp b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_force_cp2k_mixer.inp index ea5f547233..9fe90bf5a1 100644 --- a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_force_cp2k_mixer.inp +++ b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_force_cp2k_mixer.inp @@ -16,9 +16,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER CP2K &TBLITE METHOD GFN2 - SCC_MIXER CP2K &END TBLITE &END XTB &END QS diff --git a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_reference_cli_lone.inp b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_reference_cli_lone.inp index 94a2c6457a..d1eeaa896a 100644 --- a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_reference_cli_lone.inp +++ b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_reference_cli_lone.inp @@ -17,10 +17,24 @@ &REFERENCE_CLI CHECK_FORCES F CHECK_VIRIAL F + EFIELD 0.0 0.0 0.0 + GRAD tblite-reference-test.grad + JSON tblite-reference-test.json + POST_PROCESSING molmom + POST_PROCESSING_OUTPUT tblite-reference-test.npz PROGRAM_NAME /bin/false STOP_ON_ERROR F + &IMPLICIT_SOLVATION + BORN_KERNEL P16 + MODEL ALPB + SOLUTION_STATE BAR1MOL + SOLVENT water + &END IMPLICIT_SOLVATION &END REFERENCE_CLI &END TBLITE + &TBLITE_MIXER + SOLVER GVD + &END TBLITE_MIXER &END XTB &END QS &SCF diff --git a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_uks_force.inp b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_uks_force.inp index 18f7815e0e..1f6edf7a79 100644 --- a/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_uks_force.inp +++ b/tests/xTB/regtest-tblite-gfn2/CH2O_gfn2_uks_force.inp @@ -21,9 +21,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER TBLITE &TBLITE METHOD GFN2 - SCC_MIXER TBLITE &END TBLITE &END XTB &END QS diff --git a/tests/xTB/regtest-tblite-gfn2/O2_gfn2_uks_force.inp b/tests/xTB/regtest-tblite-gfn2/O2_gfn2_uks_force.inp index f2649b0056..57b7a93e34 100644 --- a/tests/xTB/regtest-tblite-gfn2/O2_gfn2_uks_force.inp +++ b/tests/xTB/regtest-tblite-gfn2/O2_gfn2_uks_force.inp @@ -20,9 +20,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER TBLITE &TBLITE METHOD GFN2 - SCC_MIXER TBLITE &END TBLITE &END XTB &END QS diff --git a/tests/xTB/regtest-tblite-gfn2/O2plus_gfn2_uks_cp2k_mixer.inp b/tests/xTB/regtest-tblite-gfn2/O2plus_gfn2_uks_cp2k_mixer.inp new file mode 100644 index 0000000000..5d352b6538 --- /dev/null +++ b/tests/xTB/regtest-tblite-gfn2/O2plus_gfn2_uks_cp2k_mixer.inp @@ -0,0 +1,49 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT O2plus_gfn2_uks_cp2k_mixer + RUN_TYPE ENERGY_FORCE +&END GLOBAL + +&FORCE_EVAL + &DFT + CHARGE 1 + LSD + MULTIPLICITY 2 + &QS + EPS_DEFAULT 1.00E-12 + METHOD xTB + &XTB + GFN_TYPE TBLITE + SCC_MIXER CP2K + &TBLITE + METHOD GFN2 + &END TBLITE + &END XTB + &END QS + &SCF + ADDED_MOS 1 2 + EPS_SCF 1.e-8 + MAX_SCF 80 + SCF_GUESS MOPAC + &MIXING + ALPHA 0.4 + METHOD BROYDEN_MIXING + NBUFFER 8 + &END MIXING + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 20.0 20.0 20.0 + PERIODIC NONE + &END CELL + &COORD + O 0.000000 0.000000 0.000000 + O 1.208000 0.000000 0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/xTB/regtest-tblite-gfn2/TEST_FILES.toml b/tests/xTB/regtest-tblite-gfn2/TEST_FILES.toml index a1dcd8ec1a..0500d2d438 100644 --- a/tests/xTB/regtest-tblite-gfn2/TEST_FILES.toml +++ b/tests/xTB/regtest-tblite-gfn2/TEST_FILES.toml @@ -10,6 +10,7 @@ "CH2O_gfn2_force.inp" = [{matcher="M082", tol=5.0E-6, ref=0.0000000}] "CH2O_gfn2_uks_force.inp" = [{matcher="M082", tol=5.0E-6, ref=0.0000000}] "O2_gfn2_uks_force.inp" = [{matcher="M082", tol=5.0E-6, ref=0.0000000}] +"O2plus_gfn2_uks_cp2k_mixer.inp" = [{matcher="E_total", tol=1.0E-8, ref=-7.2659684002}] "CH2O_gfn2_force_cp2k_mixer.inp" = [{matcher="M082", tol=5.0E-6, ref=0.0000000}] "CH2O_gfn2_stress.inp" = [{matcher="M042", tol=1.0E-7, ref=0.0000000}] "CH2O_gfn2_reference_cli_lone.inp" = [{matcher="E_total", tol=1.0E-8, ref=-7.1737299847}, diff --git a/tests/xTB/regtest-tblite-ipea1/O2_ipea1_uks_force.inp b/tests/xTB/regtest-tblite-ipea1/O2_ipea1_uks_force.inp index c9334ebb17..7b09ad7766 100644 --- a/tests/xTB/regtest-tblite-ipea1/O2_ipea1_uks_force.inp +++ b/tests/xTB/regtest-tblite-ipea1/O2_ipea1_uks_force.inp @@ -20,9 +20,9 @@ METHOD xTB &XTB GFN_TYPE TBLITE + SCC_MIXER TBLITE &TBLITE METHOD IPEA1 - SCC_MIXER TBLITE &END TBLITE &END XTB &END QS diff --git a/tools/doxify/fixcomments.pl b/tools/doxify/fixcomments.pl index 0beb3cb25e..287cfe0974 100755 --- a/tools/doxify/fixcomments.pl +++ b/tools/doxify/fixcomments.pl @@ -506,7 +506,8 @@ sub processSubroutineDefinition { # Split the subroutine or function definition by space or comma my @string = split(/([,\(\)\s+]+)/x, $functionLine); - while (my ($idx, $p) = each @string) { + for (my $idx = 0; $idx < scalar(@string); $idx++) { + my $p = $string[$idx]; defined $p or next; # continue with the next item if $p is undefined $p =~ s/^\s+|\s+$//gx;