From 70ce0ff93c2c69218a8a8dad7f082561df0da4de Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Fri, 16 Sep 2022 09:41:52 +0200 Subject: [PATCH] xTB short range Coulomb: use special neighbor list (#2287) this solves the electronic energy dependence on EPS_DEFAULT --- src/input_cp2k_dft.F | 9 +++++++++ src/input_cp2k_tb.F | 4 ++-- src/qs_environment.F | 4 ++-- src/qs_environment_types.F | 12 +++++++----- src/qs_ks_types.F | 19 +++++++++++++------ src/qs_neighbor_lists.F | 18 +++++++++++++++++- src/qs_tddfpt2_stda_types.F | 20 +++++--------------- src/qs_tddfpt2_stda_utils.F | 3 ++- src/xtb_coulomb.F | 24 ++++++++++++++++-------- src/xtb_ehess.F | 24 ++++++++++++++++-------- src/xtb_ehess_force.F | 24 ++++++++++++++++-------- tests/QS/regtest-ri-rpa-grad/TEST_FILES | 2 +- tests/QS/regtest-spgr/quartz.inp | 1 + tests/QS/regtest-stda-force/TEST_FILES | 2 +- tests/QS/regtest-tddfpt-stda/TEST_FILES | 4 ++-- tests/xTB/regtest-1/TEST_FILES | 6 +++--- tests/xTB/regtest-2/TEST_FILES | 16 ++++++++-------- tests/xTB/regtest-3/TEST_FILES | 6 +++--- tests/xTB/regtest-4/TEST_FILES | 24 ++++++++++++------------ tests/xTB/regtest-5/TEST_FILES | 4 ++-- tests/xTB/regtest-debug/TEST_FILES | 6 +++--- tests/xTB/regtest-stda-force/TEST_FILES | 2 +- tests/xTB/regtest-stda/TEST_FILES | 4 ++-- 23 files changed, 144 insertions(+), 94 deletions(-) diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 1d72c1fd39..39508e0af7 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -805,6 +805,15 @@ CONTAINS CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, & + name="sab_xtbe", & + description="Activates the printing of the xTB sr-Coulomb "// & + "neighbor lists ", & + default_l_val=.FALSE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, & name="sab_core", & description="Activates the printing of core interaction "// & diff --git a/src/input_cp2k_tb.F b/src/input_cp2k_tb.F index 7c2d276813..e37427ba51 100644 --- a/src/input_cp2k_tb.F +++ b/src/input_cp2k_tb.F @@ -420,13 +420,13 @@ CONTAINS CALL keyword_create(keyword, __LOCATION__, name="COULOMB_SR_CUT", & description="Maximum range of short range part of Coulomb interaction.", & usage="COULOMB_SR_CUT 20.0 ", repeats=.FALSE., & - n_var=1, default_r_val=12.0_dp) + n_var=1, default_r_val=20.0_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="COULOMB_SR_EPS", & description="Cutoff for short range part of Coulomb interaction.", & - usage="COULOMB_SR_EPS 1.R-3 ", repeats=.FALSE., & + usage="COULOMB_SR_EPS 1.E-3 ", repeats=.FALSE., & n_var=1, default_r_val=1.0E-03_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/qs_environment.F b/src/qs_environment.F index 8d24a89d33..ec05f5dab1 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -992,11 +992,11 @@ CONTAINS qs_kind => qs_kind_set(ikind) IF (qs_kind%xtb_parameter%defined) THEN CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set) - CALL get_gto_basis_set(tmp_basis_set, kind_radius=rcut) IF (xtb_control%old_coulomb_damping) THEN + CALL get_gto_basis_set(tmp_basis_set, kind_radius=rcut) qs_kind%xtb_parameter%rcut = rcut ELSE - rcut = MIN(rcut, xtb_control%coulomb_sr_cut) + rcut = xtb_control%coulomb_sr_cut fxx = 2.0_dp*xtb_control%coulomb_sr_eps*qs_kind%xtb_parameter%eta**2 fxx = 0.80_dp*(1.0_dp/fxx)**0.3333_dp qs_kind%xtb_parameter%rcut = MIN(rcut, fxx) diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 8e794665eb..d72aa75e37 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -361,6 +361,7 @@ CONTAINS !> \param sap_oce ... !> \param sab_lrc ... !> \param sab_se ... +!> \param sab_xtbe ... !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... @@ -494,9 +495,9 @@ CONTAINS !> \version 1.0 ! ************************************************************************************************** SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, & - dft_control, mos, sab_orb, & - sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, & - sab_se, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, particle_set, energy, force, & + dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, & + sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, & + sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, particle_set, energy, force, & matrix_h, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, & matrix_h_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, & matrix_s, matrix_s_RI_aux, matrix_w, & @@ -534,8 +535,8 @@ CONTAINS OPTIONAL, POINTER :: sab_orb, sab_all LOGICAL, OPTIONAL :: qmmm, qmmm_periodic TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sac_ae, sac_ppl, sac_lri, & - sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, & - sab_xtb_nonbond, sab_almo, sab_kp + sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, & + sab_xb, sab_xtb_nonbond, sab_almo, sab_kp TYPE(particle_type), DIMENSION(:), OPTIONAL, & POINTER :: particle_set TYPE(qs_energy_type), OPTIONAL, POINTER :: energy @@ -800,6 +801,7 @@ CONTAINS sab_se=sab_se, & sab_lrc=sab_lrc, & sab_tbe=sab_tbe, & + sab_xtbe=sab_xtbe, & sab_core=sab_core, & sab_xb=sab_xb, & sab_xtb_nonbond=sab_xtb_nonbond, & diff --git a/src/qs_ks_types.F b/src/qs_ks_types.F index d8e7bec49d..f95eeb6355 100644 --- a/src/qs_ks_types.F +++ b/src/qs_ks_types.F @@ -183,6 +183,7 @@ MODULE qs_ks_types sab_vdw => Null(), & sab_scp => Null(), & sab_tbe => Null(), & + sab_xtbe => Null(), & sab_core => Null(), & sab_xb => Null(), & sab_xtb_nonbond => Null(), & @@ -275,6 +276,7 @@ CONTAINS !> \param sap_oce ... !> \param sab_lrc ... !> \param sab_se ... +!> \param sab_xtbe ... !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... @@ -328,7 +330,7 @@ CONTAINS vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, & neighbor_list_id, & sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, & - sab_se, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & + sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & task_list, task_list_soft, & kpoints, do_kpoints, & atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, & @@ -352,8 +354,8 @@ CONTAINS TYPE(pw_p_type), OPTIONAL, POINTER :: vppl, rho_core, rho_nlcc, rho_nlcc_g, vee INTEGER, OPTIONAL :: neighbor_list_id TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_all, sac_ae, & - sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, & - sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp + sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, & + sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_soft TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints LOGICAL, OPTIONAL :: do_kpoints @@ -438,6 +440,7 @@ CONTAINS IF (PRESENT(sab_se)) sab_se => ks_env%sab_se IF (PRESENT(sab_lrc)) sab_lrc => ks_env%sab_lrc IF (PRESENT(sab_tbe)) sab_tbe => ks_env%sab_tbe + IF (PRESENT(sab_xtbe)) sab_xtbe => ks_env%sab_xtbe IF (PRESENT(sab_core)) sab_core => ks_env%sab_core IF (PRESENT(sab_xb)) sab_xb => ks_env%sab_xb IF (PRESENT(sab_xtb_nonbond)) sab_xtb_nonbond => ks_env%sab_xtb_nonbond @@ -524,6 +527,7 @@ CONTAINS !> \param sap_oce ... !> \param sab_lrc ... !> \param sab_se ... +!> \param sab_xtbe ... !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... @@ -554,7 +558,7 @@ CONTAINS neighbor_list_id, & kpoints, & sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, & - sab_se, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & + sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & task_list, task_list_soft, & subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, & para_env, blacs_env) @@ -574,8 +578,8 @@ CONTAINS INTEGER, OPTIONAL :: neighbor_list_id TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_all, sac_ae, & - sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, & - sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp + sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, & + sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_soft TYPE(qs_subsys_type), OPTIONAL, POINTER :: subsys TYPE(dft_control_type), OPTIONAL, POINTER :: dft_control @@ -639,6 +643,7 @@ CONTAINS IF (PRESENT(sab_se)) ks_env%sab_se => sab_se IF (PRESENT(sab_lrc)) ks_env%sab_lrc => sab_lrc IF (PRESENT(sab_tbe)) ks_env%sab_tbe => sab_tbe + IF (PRESENT(sab_xtbe)) ks_env%sab_xtbe => sab_xtbe IF (PRESENT(sab_core)) ks_env%sab_core => sab_core IF (PRESENT(sab_xb)) ks_env%sab_xb => sab_xb IF (PRESENT(sab_xtb_nonbond)) ks_env%sab_xtb_nonbond => sab_xtb_nonbond @@ -767,6 +772,7 @@ CONTAINS CALL release_neighbor_list_sets(ks_env%sab_vdw) CALL release_neighbor_list_sets(ks_env%sab_scp) CALL release_neighbor_list_sets(ks_env%sab_tbe) + CALL release_neighbor_list_sets(ks_env%sab_xtbe) CALL release_neighbor_list_sets(ks_env%sab_core) CALL release_neighbor_list_sets(ks_env%sab_xb) CALL release_neighbor_list_sets(ks_env%sab_xtb_nonbond) @@ -856,6 +862,7 @@ CONTAINS CALL release_neighbor_list_sets(ks_env%sab_vdw) CALL release_neighbor_list_sets(ks_env%sab_scp) CALL release_neighbor_list_sets(ks_env%sab_tbe) + CALL release_neighbor_list_sets(ks_env%sab_xtbe) CALL release_neighbor_list_sets(ks_env%sab_core) CALL release_neighbor_list_sets(ks_env%sab_xb) CALL release_neighbor_list_sets(ks_env%sab_xtb_nonbond) diff --git a/src/qs_neighbor_lists.F b/src/qs_neighbor_lists.F index 097c98e98b..4fdc34aa92 100644 --- a/src/qs_neighbor_lists.F +++ b/src/qs_neighbor_lists.F @@ -95,6 +95,8 @@ MODULE qs_neighbor_lists subcell_type USE util, ONLY: locate,& sort + USE xtb_types, ONLY: get_xtb_atom_param,& + xtb_atom_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -327,7 +329,8 @@ CONTAINS TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, & sab_cn, sab_core, sab_gcp, sab_kp, sab_lrc, sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, & - sab_xb, sab_xtb_nonbond, sac_ae, sac_lri, sac_ppl, sap_oce, sap_ppnl, soa_list, soo_list + sab_xb, sab_xtb_nonbond, sab_xtbe, sac_ae, sac_lri, sac_ppl, sap_oce, sap_ppnl, soa_list, & + soo_list TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(paw_proj_set_type), POINTER :: paw_proj TYPE(qs_dftb_atom_type), POINTER :: dftb_atom @@ -337,6 +340,7 @@ CONTAINS TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section TYPE(sgp_potential_type), POINTER :: sgp_potential + TYPE(xtb_atom_type), POINTER :: xtb_atom CALL timeset(routineN, handle) NULLIFY (logger) @@ -355,6 +359,7 @@ CONTAINS NULLIFY (sab_se) NULLIFY (sab_lrc) NULLIFY (sab_tbe) + NULLIFY (sab_xtbe) NULLIFY (sab_core) NULLIFY (sab_xb) NULLIFY (sab_xtb_nonbond) @@ -396,6 +401,7 @@ CONTAINS sab_se=sab_se, & sab_lrc=sab_lrc, & sab_tbe=sab_tbe, & + sab_xtbe=sab_xtbe, & sab_core=sab_core, & sab_xb=sab_xb, & sab_xtb_nonbond=sab_xtb_nonbond, & @@ -757,6 +763,16 @@ CONTAINS subcells=subcells, nlname="sab_tbe") CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe) END IF + ! SR part of Coulomb interaction + DO ikind = 1, nkind + CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom) + CALL get_xtb_atom_param(xtb_parameter=xtb_atom, rcut=c_radius(ikind)) + END DO + default_present = .TRUE. + CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius) + CALL build_neighbor_lists(sab_xtbe, particle_set, atom2d, cell, pair_radius, & + subcells=subcells, nlname="sab_xtbe") + CALL set_ks_env(ks_env=ks_env, sab_xtbe=sab_xtbe) ! XB list ALLOCATE (xb1_atom(nkind), xb2_atom(nkind)) c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius diff --git a/src/qs_tddfpt2_stda_types.F b/src/qs_tddfpt2_stda_types.F index 8e6f56b916..ff12cd534b 100644 --- a/src/qs_tddfpt2_stda_types.F +++ b/src/qs_tddfpt2_stda_types.F @@ -11,8 +11,6 @@ MODULE qs_tddfpt2_stda_types USE atomic_kind_types, ONLY: atomic_kind_type,& get_atomic_kind - USE basis_set_types, ONLY: get_gto_basis_set,& - gto_basis_set_type USE cp_control_types, ONLY: dft_control_type,& stda_control_type USE cp_log_handling, ONLY: cp_get_default_logger,& @@ -25,8 +23,7 @@ MODULE qs_tddfpt2_stda_types USE physcon, ONLY: evolt USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type - USE qs_kind_types, ONLY: get_qs_kind,& - qs_kind_type + USE qs_kind_types, ONLY: qs_kind_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -130,9 +127,6 @@ CONTAINS TYPE(atomic_kind_type), POINTER :: atomic_kind TYPE(cp_logger_type), POINTER :: logger TYPE(dft_control_type), POINTER :: dft_control - TYPE(gto_basis_set_type), POINTER :: tmp_basis_set - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set - TYPE(qs_kind_type), POINTER :: qs_kind TYPE(section_vals_type), POINTER :: tddfpt_print_section TYPE(stda_kind_type), POINTER :: kind_param @@ -141,10 +135,9 @@ CONTAINS CPASSERT(ASSOCIATED(stda_kernel%kind_param_set)) - NULLIFY (atomic_kind_set, qs_kind_set) !get element symbol and atomic number - CALL get_qs_env(qs_env, dft_control=dft_control, & - atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) - nkind = SIZE(qs_kind_set) + NULLIFY (atomic_kind_set) + CALL get_qs_env(qs_env, dft_control=dft_control, atomic_kind_set=atomic_kind_set) + nkind = SIZE(atomic_kind_set) NULLIFY (tddfpt_print_section) tddfpt_print_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT%PRINT") @@ -197,13 +190,10 @@ CONTAINS stda_kernel%kind_param_set(ikind)%kind_param%hardness_param = & hardness(stda_kernel%kind_param_set(ikind)%kind_param%z)*2.0_dp/evolt ! rcut parameter - qs_kind => qs_kind_set(ikind) - CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set) - CALL get_gto_basis_set(tmp_basis_set, kind_radius=rcut) eta = stda_kernel%kind_param_set(ikind)%kind_param%hardness_param fxx = 2.0_dp*eta**2*stda_control%coulomb_sr_eps fxx = 0.5_dp*(1.0_dp/fxx)**0.33333_dp - rcut = MIN(rcut, stda_control%coulomb_sr_cut) + rcut = stda_control%coulomb_sr_cut stda_kernel%kind_param_set(ikind)%kind_param%rcut = MIN(rcut, fxx) END DO diff --git a/src/qs_tddfpt2_stda_utils.F b/src/qs_tddfpt2_stda_utils.F index 31ae9c11e8..692fac93f8 100644 --- a/src/qs_tddfpt2_stda_utils.F +++ b/src/qs_tddfpt2_stda_utils.F @@ -180,8 +180,9 @@ CONTAINS CALL get_qs_env(qs_env=qs_env, natom=natom) dbcsr_dist => sub_env%dbcsr_dist + ! Using the overlap list here can have a considerable effect on the number of + ! terms calculated. This makes gamma also dependent on EPS_DEFAULT -> Overlap n_list => sub_env%sab_orb - CALL get_qs_env(qs_env=qs_env, dbcsr_dist=dbcsr_dist, sab_orb=n_list) IF (PRESENT(ndim)) THEN nmat = ndim diff --git a/src/xtb_coulomb.F b/src/xtb_coulomb.F index f2d2d13a74..ba15ed88b0 100644 --- a/src/xtb_coulomb.F +++ b/src/xtb_coulomb.F @@ -22,7 +22,8 @@ MODULE xtb_coulomb USE cell_types, ONLY: cell_type,& get_cell,& pbc - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + xtb_control_type USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: dbcsr_add,& dbcsr_get_block_p,& @@ -153,6 +154,7 @@ CONTAINS TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int TYPE(virial_type), POINTER :: virial TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind + TYPE(xtb_control_type), POINTER :: xtb_control CALL timeset(routineN, handle) @@ -166,6 +168,8 @@ CONTAINS atprop=atprop, & dft_control=dft_control) + xtb_control => dft_control%qs_control%xtb_control + use_virial = .FALSE. IF (calculate_forces) THEN use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) @@ -195,10 +199,14 @@ CONTAINS gmcharge = 0._dp ! short range contribution (gamma) - ! loop over all atom pairs with a non-zero overlap (sab_orb) - kg = dft_control%qs_control%xtb_control%kg + ! loop over all atom pairs (sab_xtbe) + kg = xtb_control%kg NULLIFY (n_list) - CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + IF (xtb_control%old_coulomb_damping) THEN + CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + ELSE + CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list) + END IF CALL neighbor_list_iterator_create(nl_iterator, n_list) DO WHILE (neighbor_list_iterate(nl_iterator) == 0) CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & @@ -255,8 +263,8 @@ CONTAINS ! 1/R contribution - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN - do_ewald = dft_control%qs_control%xtb_control%do_ewald + IF (xtb_control%coulomb_lr) THEN + do_ewald = xtb_control%do_ewald IF (do_ewald) THEN ! Ewald sum NULLIFY (ewald_env, ewald_pw) @@ -312,7 +320,7 @@ CONTAINS CALL mp_sum(gmcharge(:, 1), para_env%group) CALL mp_sum(gchrg(:, :, 1), para_env%group) - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN + IF (xtb_control%coulomb_lr) THEN IF (do_ewald) THEN ! add self charge interaction and background charge contribution gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:) @@ -633,7 +641,7 @@ CONTAINS END IF END IF - IF (dft_control%qs_control%xtb_control%tb3_interaction) THEN + IF (xtb_control%tb3_interaction) THEN CALL get_qs_env(qs_env, nkind=nkind) ALLOCATE (zeffk(nkind), xgamma(nkind)) DO ikind = 1, nkind diff --git a/src/xtb_ehess.F b/src/xtb_ehess.F index 5b3d08bba4..d04c6cfe30 100644 --- a/src/xtb_ehess.F +++ b/src/xtb_ehess.F @@ -16,7 +16,8 @@ MODULE xtb_ehess USE cell_types, ONLY: cell_type,& get_cell,& pbc - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + xtb_control_type USE cp_para_types, ONLY: cp_para_env_type USE dbcsr_api, ONLY: dbcsr_get_block_p,& dbcsr_iterator_blocks_left,& @@ -115,6 +116,7 @@ CONTAINS TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(virial_type), POINTER :: virial TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind + TYPE(xtb_control_type), POINTER :: xtb_control CALL timeset(routineN, handle) @@ -125,6 +127,8 @@ CONTAINS cell=cell, & dft_control=dft_control) + xtb_control => dft_control%qs_control%xtb_control + IF (dft_control%nimages /= 1) THEN CPABORT("No kpoints allowed in xTB response calculation") END IF @@ -137,10 +141,14 @@ CONTAINS gmcharge = 0._dp ! short range contribution (gamma) - ! loop over all atom pairs with a non-zero overlap (sab_orb) - kg = dft_control%qs_control%xtb_control%kg + ! loop over all atom pairs (sab_xtbe) + kg = xtb_control%kg NULLIFY (n_list) - CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + IF (xtb_control%old_coulomb_damping) THEN + CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + ELSE + CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list) + END IF CALL neighbor_list_iterator_create(nl_iterator, n_list) DO WHILE (neighbor_list_iterate(nl_iterator) == 0) CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & @@ -171,8 +179,8 @@ CONTAINS ! 1/R contribution - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN - do_ewald = dft_control%qs_control%xtb_control%do_ewald + IF (xtb_control%coulomb_lr) THEN + do_ewald = xtb_control%do_ewald IF (do_ewald) THEN ! Ewald sum NULLIFY (ewald_env, ewald_pw) @@ -222,7 +230,7 @@ CONTAINS CALL mp_sum(gmcharge(:, 1), para_env%group) CALL mp_sum(gchrg(:, :, 1), para_env%group) - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN + IF (xtb_control%coulomb_lr) THEN IF (do_ewald) THEN ! add self charge interaction and background charge contribution gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:) @@ -272,7 +280,7 @@ CONTAINS END DO CALL dbcsr_iterator_stop(iter) - IF (dft_control%qs_control%xtb_control%tb3_interaction) THEN + IF (xtb_control%tb3_interaction) THEN CALL get_qs_env(qs_env, nkind=nkind) ALLOCATE (xgamma(nkind)) DO ikind = 1, nkind diff --git a/src/xtb_ehess_force.F b/src/xtb_ehess_force.F index de669697f4..8cc1859f73 100644 --- a/src/xtb_ehess_force.F +++ b/src/xtb_ehess_force.F @@ -16,7 +16,8 @@ MODULE xtb_ehess_force USE cell_types, ONLY: cell_type,& get_cell,& pbc - USE cp_control_types, ONLY: dft_control_type + USE cp_control_types, ONLY: dft_control_type,& + xtb_control_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_unit_nr,& cp_logger_type @@ -134,6 +135,7 @@ CONTAINS TYPE(qs_rho_type), POINTER :: rho TYPE(virial_type), POINTER :: virial TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind + TYPE(xtb_control_type), POINTER :: xtb_control CALL timeset(routineN, handle) @@ -156,6 +158,8 @@ CONTAINS atprop=atprop, & dft_control=dft_control) + xtb_control => dft_control%qs_control%xtb_control + calculate_forces = .TRUE. just_energy = .FALSE. use_virial = .FALSE. @@ -176,10 +180,14 @@ CONTAINS gmcharge1 = 0._dp ! short range contribution (gamma) - ! loop over all atom pairs with a non-zero overlap (sab_orb) - kg = dft_control%qs_control%xtb_control%kg + ! loop over all atom pairs (sab_xtbe) + kg = xtb_control%kg NULLIFY (n_list) - CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + IF (xtb_control%old_coulomb_damping) THEN + CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + ELSE + CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list) + END IF CALL neighbor_list_iterator_create(nl_iterator, n_list) DO WHILE (neighbor_list_iterate(nl_iterator) == 0) CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & @@ -227,8 +235,8 @@ CONTAINS ! 1/R contribution - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN - do_ewald = dft_control%qs_control%xtb_control%do_ewald + IF (xtb_control%coulomb_lr) THEN + do_ewald = xtb_control%do_ewald IF (do_ewald) THEN ! Ewald sum NULLIFY (ewald_env, ewald_pw) @@ -290,7 +298,7 @@ CONTAINS CALL mp_sum(gmcharge1(:, 1), para_env%group) CALL mp_sum(gchrg1(:, :, 1), para_env%group) - IF (dft_control%qs_control%xtb_control%coulomb_lr) THEN + IF (xtb_control%coulomb_lr) THEN IF (do_ewald) THEN ! add self charge interaction and background charge contribution gmcharge0(:, 1) = gmcharge0(:, 1) - 2._dp*alpha*oorootpi*mcharge0(:) @@ -411,7 +419,7 @@ CONTAINS IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H[P]*dS ", fodeb END IF - IF (dft_control%qs_control%xtb_control%tb3_interaction) THEN + IF (xtb_control%tb3_interaction) THEN CALL get_qs_env(qs_env, nkind=nkind) ALLOCATE (xgamma(nkind)) DO ikind = 1, nkind diff --git a/tests/QS/regtest-ri-rpa-grad/TEST_FILES b/tests/QS/regtest-ri-rpa-grad/TEST_FILES index 52c6aa68a4..a94f83e169 100644 --- a/tests/QS/regtest-ri-rpa-grad/TEST_FILES +++ b/tests/QS/regtest-ri-rpa-grad/TEST_FILES @@ -1,4 +1,4 @@ -RI_RPA_grad_MP2_H2O.inp 11 4e-08 -16.994626142112192 +RI_RPA_grad_MP2_H2O.inp 11 5e-08 -16.994626142112192 RI_RPA_grad_MP2_CH3.inp 11 7e-08 -7.515645339785070 RI_RPA_grad_MP2_H2O_HF_PBE.inp 11 6e-08 -17.057888508919117 RI_RPA_grad_MP2_H2O_HF_PBE_ADMM.inp 11 1e-07 -17.066276945129722 diff --git a/tests/QS/regtest-spgr/quartz.inp b/tests/QS/regtest-spgr/quartz.inp index 0435a75e4c..ba453d302d 100644 --- a/tests/QS/regtest-spgr/quartz.inp +++ b/tests/QS/regtest-spgr/quartz.inp @@ -36,6 +36,7 @@ METHOD XTB &XTB DO_EWALD T + CHECK_ATOMIC_CHARGES F &END XTB &END QS &POISSON diff --git a/tests/QS/regtest-stda-force/TEST_FILES b/tests/QS/regtest-stda-force/TEST_FILES index a872226637..1b558809bf 100644 --- a/tests/QS/regtest-stda-force/TEST_FILES +++ b/tests/QS/regtest-stda-force/TEST_FILES @@ -14,5 +14,5 @@ h2o_f07.inp 1 1.0E-12 h2o_f08a.inp 1 1.0E-12 -17.23033252463928 h2o_f08b.inp 1 1.0E-12 -17.23207360292983 CH2O_stda-xtb-s_virtual_shift.inp 37 1.0E-06 0.779158E-01 -water_xTB_virtual_shift.inp 37 1.0E-06 0.236248E+00 +water_xTB_virtual_shift.inp 37 1.0E-06 0.236164E+00 #EOF diff --git a/tests/QS/regtest-tddfpt-stda/TEST_FILES b/tests/QS/regtest-tddfpt-stda/TEST_FILES index 4ee261db1f..629c7abeb2 100644 --- a/tests/QS/regtest-tddfpt-stda/TEST_FILES +++ b/tests/QS/regtest-tddfpt-stda/TEST_FILES @@ -17,6 +17,6 @@ H2O_tddfpt_stda-pbe0-s.inp 37 1.0E-05 H2O_tddfpt_stda-pbe0-t.inp 37 1.0E-05 0.633489E+00 H2O_tddfpt_stda-s-1.inp 37 1.0E-05 0.173021E+00 NO_tddfpt_stda-s-1.inp 37 1.0E-05 0.375613E+00 -water_xTB.inp 37 2.0E-05 0.126001E+00 -water_xTBi_NTO.inp 37 2.0E-05 0.126001E+00 +water_xTB.inp 37 2.0E-05 0.125917E+00 +water_xTBi_NTO.inp 37 2.0E-05 0.125917E+00 #EOF diff --git a/tests/xTB/regtest-1/TEST_FILES b/tests/xTB/regtest-1/TEST_FILES index bbfcd96e40..ec1e6d4e72 100644 --- a/tests/xTB/regtest-1/TEST_FILES +++ b/tests/xTB/regtest-1/TEST_FILES @@ -10,9 +10,9 @@ ch2o_smear.inp 1 1.0E-12 -7.19650182 tmol.inp 1 1.0E-12 -41.90845660778482 h2.inp 1 1.0E-12 -1.03458223634251 h2_kab.inp 1 1.0E-12 -0.99816707611986 -h2o-md.inp 1 1.0E-12 -185.14802701162756 -h2o_str.inp 1 1.0E-12 -5.76544426551662 -h2o-atprop.inp 1 1.0E-10 -185.16704903931230 +h2o-md.inp 1 1.0E-12 -185.14385714384289 +h2o_str.inp 1 1.0E-12 -5.76524673249110 +h2o-atprop.inp 1 1.0E-10 -185.16262201741714 h2o-atprop0.inp 1 1.0E-12 -187.45499307089042 si_geo.inp 1 1.0E-12 -14.51194886943525 si_kp.inp 1 1.0E-12 -14.73191413260763 diff --git a/tests/xTB/regtest-2/TEST_FILES b/tests/xTB/regtest-2/TEST_FILES index e3c7da8672..6e1e0f28e6 100644 --- a/tests/xTB/regtest-2/TEST_FILES +++ b/tests/xTB/regtest-2/TEST_FILES @@ -4,18 +4,18 @@ # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest # test mulliken constraints -H2O-field-gopt.inp 1 3e-11 -5.76960964769063 -H2O-field-gopt-lsd.inp 1 4e-12 -5.76960970095199 +H2O-field-gopt.inp 1 3e-11 -5.76970795619747 +H2O-field-gopt-lsd.inp 1 4e-12 -5.76970797174834 H2O-field.inp 0 4e-13 -H2O-field-lsd.inp 1 4e-14 -5.76948986131183 -HF-field.inp 1 1e-12 -5.62817044142137 -HF-field-gopt.inp 1 5e-09 -5.66333179696722 +H2O-field-lsd.inp 1 4e-14 -5.76959233315201 +HF-field.inp 1 1e-12 -5.62813329327978 +HF-field-gopt.inp 1 5e-09 -5.66367831097363 HF-field-debug.inp 0 1e-12 HF-dfilter-debug.inp 0 1e-12 -HF-dfield-gopt.inp 1 1e-09 -5.66113975858277 -HF-dfield.inp 1 1e-12 -5.66099020966660 +HF-dfield-gopt.inp 1 1e-09 -5.66065888870620 +HF-dfield.inp 1 1e-12 -5.66052242442153 HF-dfield-debug.inp 0 1e-12 HF-loc-field.inp 1 1e-12 -5.65286303322644 -HF-loc-field-gopt.inp 1 5e-09 -5.67065466488875 +HF-loc-field-gopt.inp 1 5e-09 -5.67004654243380 HF-loc-field-debug.inp 0 1e-12 #EOF diff --git a/tests/xTB/regtest-3/TEST_FILES b/tests/xTB/regtest-3/TEST_FILES index f50e961683..009dbe05e4 100644 --- a/tests/xTB/regtest-3/TEST_FILES +++ b/tests/xTB/regtest-3/TEST_FILES @@ -4,8 +4,8 @@ # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest NdF3.inp 1 1.0E-12 -16.30904020314352 -h2o_rtp.inp 1 1.0E-12 -5.76549019124211 -h2o_emd.inp 1 1.0E-12 -5.76683221354591 +h2o_rtp.inp 1 1.0E-12 -5.76539152285971 +h2o_emd.inp 1 1.0E-12 -5.76621483846123 si8_wan.inp 1 1.0E-12 -14.36325382143829 si_kp.inp 1 1.0E-12 -14.73032998731771 tmol.inp 1 1.0E-12 -41.90845660778482 @@ -18,6 +18,6 @@ ch2o_mos.inp 1 1.0E-12 -7.84456570 si_print.inp 1 1.0E-12 -14.73032998731853 si_band.inp 1 1.0E-12 -14.73032998731853 H2O-geo-pdos.inp 1 1.0E-12 -5.76872484344198 -graphite-stm.inp 1 1.0E-12 -7.91352193418221 +graphite-stm.inp 1 1.0E-12 -7.91403228335400 si_dos.inp 1 1.0E-12 -14.73032998731853 #EOF diff --git a/tests/xTB/regtest-4/TEST_FILES b/tests/xTB/regtest-4/TEST_FILES index abcd52c78f..00b55b2ab0 100644 --- a/tests/xTB/regtest-4/TEST_FILES +++ b/tests/xTB/regtest-4/TEST_FILES @@ -3,18 +3,18 @@ # e.g. 0 means do not compare anything, running is enough # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest -H2O-32-xtb-ls-1.inp 11 1.0E-12 -46.141216327501148 -H2O-32-xtb-ls-2.inp 11 1.0E-12 -46.150675466617464 -H2O-32-xtb-ls-3.inp 11 1.0E-12 -92.315985244652225 -H2O-32-xtb-ls-4.inp 11 1.0E-12 -92.315984055047494 -H2O-32-xtb-ls-5.inp 11 1.0E-12 -92.315984900840121 -H2O-32-xtb-ls-6.inp 11 1.0E-12 -92.315984055047494 -H2O-32-xtb-ls-7.inp 11 1.0E-12 -92.315984880279260 -H2O-32-xtb-ls-8.inp 11 1.0E-12 -92.315984903564313 -H2O-xtb-tc2-1.inp 11 1.0E-12 -92.315984177511098 -H2O-xtb-tc2-2.inp 11 1.0E-12 -92.315984177511282 -H2O-32-xtb-trs4.inp 11 1.0E-12 -46.150675498649122 -H2O-32-xtb-trs4-dyn.inp 11 1.0E-12 -46.150675498649122 +H2O-32-xtb-ls-1.inp 11 1.0E-12 -46.140441340701784 +H2O-32-xtb-ls-2.inp 11 1.0E-12 -46.149897222623842 +H2O-32-xtb-ls-3.inp 11 1.0E-12 -92.314662149041851 +H2O-32-xtb-ls-4.inp 11 1.0E-12 -92.314660970750793 +H2O-32-xtb-ls-5.inp 11 1.0E-12 -92.314661815852020 +H2O-32-xtb-ls-6.inp 11 1.0E-12 -92.314660970750793 +H2O-32-xtb-ls-7.inp 11 1.0E-12 -92.314661795311508 +H2O-32-xtb-ls-8.inp 11 1.0E-12 -92.314661818553091 +H2O-xtb-tc2-1.inp 11 1.0E-12 -92.314661083181150 +H2O-xtb-tc2-2.inp 11 1.0E-12 -92.314661083180596 +H2O-32-xtb-trs4.inp 11 1.0E-12 -46.149897253047165 +H2O-32-xtb-trs4-dyn.inp 11 1.0E-12 -46.149897253047165 h2o_vib.inp 8 8.0E-06 1465.672 tmol.inp 11 1.0E-12 -41.908456607784821 #EOF diff --git a/tests/xTB/regtest-5/TEST_FILES b/tests/xTB/regtest-5/TEST_FILES index 63bd24cb11..85a02adb2f 100644 --- a/tests/xTB/regtest-5/TEST_FILES +++ b/tests/xTB/regtest-5/TEST_FILES @@ -5,8 +5,8 @@ # for details see cp2k/tools/do_regtest AdeThyvdW.inp 1 1.0E-12 -58.89501875071196 ice.inp 1 1.0E-10 -46.30708081122568 -ice2.inp 1 1.0E-10 -46.30638011626871 -Ru_geo.inp 1 1.0E-10 -7.01204333848255 +ice2.inp 1 1.0E-10 -46.30708080861037 +Ru_geo.inp 1 1.0E-10 -7.01245481257972 ef_stress1.inp 0 ef_stress2.inp 0 ef_stress3.inp 0 diff --git a/tests/xTB/regtest-debug/TEST_FILES b/tests/xTB/regtest-debug/TEST_FILES index fb780c334e..3eb36e9600 100644 --- a/tests/xTB/regtest-debug/TEST_FILES +++ b/tests/xTB/regtest-debug/TEST_FILES @@ -5,13 +5,13 @@ # for details see cp2k/tools/do_regtest ch2o_t01.inp 86 1.0E-12 0.108915750509E+01 ch2o_t02.inp 86 1.0E-12 0.108915750509E+01 -ch2o_t03.inp 86 1.0E-12 0.114591982142E+01 +ch2o_t03.inp 86 1.0E-12 0.114819181828E+01 ch2o_t04.inp 86 1.0E-12 0.109829662200E+01 ch2o_t05.inp 87 1.0E-12 0.226843652140E+02 -ch2o_t06.inp 87 1.0E-12 0.230117849452E+02 +ch2o_t06.inp 87 1.0E-12 0.229761888713E+02 ch2o_t07.inp 86 1.0E-12 0.954101834865E+00 ch2o_t08.inp 86 1.0E-12 0.963128388956E+00 ch3br_nonbond_2.inp 1 1.0E-12 -11.78062819229165 ch3br_nonbond.inp 1 1.0E-12 -11.78062802267188 -ch3br_atprop_nonbond.inp 1 1.0E-10 -11.78112872152760 +ch3br_atprop_nonbond.inp 1 1.0E-10 -11.78112909303997 #EOF diff --git a/tests/xTB/regtest-stda-force/TEST_FILES b/tests/xTB/regtest-stda-force/TEST_FILES index 083270c7fe..cea702fd45 100644 --- a/tests/xTB/regtest-stda-force/TEST_FILES +++ b/tests/xTB/regtest-stda-force/TEST_FILES @@ -10,6 +10,6 @@ h2o_f14.inp 1 1.0E-12 #h2o_f15.inp 1 1.0E-12 0.000000E+00 #h2o_f16.inp 1 1.0E-12 0.000000E+00 #h2o_f17.inp 1 1.0E-12 0.000000E+00 -h2o_f18.inp 1 1.0E-12 -5.77037437397448 +h2o_f18.inp 1 1.0E-12 -5.76959310723082 #h2o_f19.inp 1 1.0E-12 0.000000E+00 #EOF diff --git a/tests/xTB/regtest-stda/TEST_FILES b/tests/xTB/regtest-stda/TEST_FILES index 2126b21768..6ee9801b4d 100644 --- a/tests/xTB/regtest-stda/TEST_FILES +++ b/tests/xTB/regtest-stda/TEST_FILES @@ -4,6 +4,6 @@ # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest # -water_xTB.inp 37 2.0E-05 0.126001E+00 -water_xTB_2.inp 37 2.0E-05 0.186902E+00 +water_xTB.inp 37 2.0E-05 0.125917E+00 +water_xTB_2.inp 37 2.0E-05 0.186821E+00 #EOF