From fe04ab922dcd36428ef1171000fc68acc8d83b89 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Wed, 21 Oct 2020 14:51:22 +0200 Subject: [PATCH] xTB pairpotential (A. Hehn) (#1136) * xTB pairpotential (A. Hehn) * Adjust regtest * Adjust regtests * adjust test --- src/cp_control_types.F | 9 + src/cp_control_utils.F | 26 +- src/input_cp2k_dft.F | 9 + src/input_cp2k_tb.F | 46 ++- src/pair_potential.F | 3 +- src/qs_energy_types.F | 4 +- src/qs_environment_types.F | 6 +- src/qs_ks_types.F | 17 +- src/qs_neighbor_lists.F | 52 +++- src/qs_scf_output.F | 6 + src/xtb_matrices.F | 289 +++++++++++++++++- tests/xTB/regtest-debug/TEST_FILES | 3 + .../regtest-debug/ch3br_atprop_nonbond.inp | 71 +++++ tests/xTB/regtest-debug/ch3br_nonbond.inp | 56 ++++ tests/xTB/regtest-debug/ch3br_nonbond_2.inp | 70 +++++ .../regtest-debug/h2o-atprop_nonbonded.inp | 176 +++++++++++ 16 files changed, 812 insertions(+), 31 deletions(-) create mode 100644 tests/xTB/regtest-debug/ch3br_atprop_nonbond.inp create mode 100644 tests/xTB/regtest-debug/ch3br_nonbond.inp create mode 100644 tests/xTB/regtest-debug/ch3br_nonbond_2.inp create mode 100644 tests/xTB/regtest-debug/h2o-atprop_nonbonded.inp diff --git a/src/cp_control_types.F b/src/cp_control_types.F index e10e9bdb58..bba580fc67 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -16,6 +16,8 @@ MODULE cp_control_types USE kinds, ONLY: default_path_length,& default_string_length,& dp + USE pair_potential_types, ONLY: pair_potential_p_release,& + pair_potential_p_type USE qs_cdft_types, ONLY: cdft_control_create,& cdft_control_release,& cdft_control_type @@ -113,6 +115,7 @@ MODULE cp_control_types REAL(KIND=dp) :: kxr, kx2 ! LOGICAL :: xb_interaction + LOGICAL :: do_nonbonded LOGICAL :: coulomb_interaction LOGICAL :: tb3_interaction LOGICAL :: check_atomic_charges @@ -124,6 +127,8 @@ MODULE cp_control_types INTEGER, DIMENSION(:, :), POINTER :: kab_types INTEGER :: kab_nval REAL, DIMENSION(:), POINTER :: kab_vals + ! + TYPE(pair_potential_p_type), POINTER :: nonbonded END TYPE xtb_control_type ! ************************************************************************************************** @@ -1052,6 +1057,7 @@ CONTAINS NULLIFY (xtb_control%kab_param) NULLIFY (xtb_control%kab_vals) NULLIFY (xtb_control%kab_types) + NULLIFY (xtb_control%nonbonded) END SUBROUTINE xtb_control_create @@ -1072,6 +1078,9 @@ CONTAINS IF (ASSOCIATED(xtb_control%kab_types)) THEN DEALLOCATE (xtb_control%kab_types) END IF + IF (ASSOCIATED(xtb_control%nonbonded)) THEN + CALL pair_potential_p_release(xtb_control%nonbonded) + END IF DEALLOCATE (xtb_control) END IF END SUBROUTINE xtb_control_release diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 4910d98820..bfc5b2ec4c 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -26,6 +26,7 @@ MODULE cp_control_utils cp_print_key_unit_nr USE cp_units, ONLY: cp_unit_from_cp2k,& cp_unit_to_cp2k + USE force_fields_input, ONLY: read_gp_section USE input_constants, ONLY: & constant_env, custom_env, do_admm_basis_projection, do_admm_blocked_projection, & do_admm_blocking_purify_full, do_admm_charge_constrained_projection, & @@ -54,6 +55,7 @@ MODULE cp_control_utils USE kinds, ONLY: default_path_length,& default_string_length,& dp + USE pair_potential_types, ONLY: pair_potential_reallocate USE periodic_table, ONLY: get_ptable_info USE qs_cdft_utils, ONLY: read_cdft_control_section USE string_utilities, ONLY: uppercase @@ -613,21 +615,21 @@ CONTAINS CHARACTER(LEN=default_string_length), & DIMENSION(:), POINTER :: clist INTEGER :: handle, itmp, j, jj, k, n_rep, n_var, & - ngauss, nrep + ngauss, ngp, nrep INTEGER, DIMENSION(:), POINTER :: tmplist LOGICAL :: explicit, was_present REAL(dp) :: tmp, tmpsqrt, value REAL(dp), POINTER :: scal(:) TYPE(section_vals_type), POINTER :: cdft_control_section, ddapc_restraint_section, & - dftb_parameter, dftb_section, lri_optbas_section, mull_section, s2_restraint_section, & - se_section, xtb_parameter, xtb_section + dftb_parameter, dftb_section, genpot_section, lri_optbas_section, mull_section, & + nonbonded_section, s2_restraint_section, se_section, xtb_parameter, xtb_section 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) + cdft_control_section, genpot_section) mull_section => section_vals_get_subs_vals(qs_section, "MULLIKEN_RESTRAINT") ddapc_restraint_section => section_vals_get_subs_vals(qs_section, "DDAPC_RESTRAINT") @@ -639,6 +641,8 @@ CONTAINS xtb_parameter => section_vals_get_subs_vals(xtb_section, "PARAMETER") 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") ! Setup all defaults values and overwrite input parameters ! EPS_DEFAULT should set the target accuracy in the total energy (~per electron) or a closely related value @@ -1106,6 +1110,18 @@ CONTAINS CALL section_vals_val_get(xtb_parameter, "HALOGEN_BINDING", r_vals=scal) qs_control%xtb_control%kxr = scal(1) qs_control%xtb_control%kx2 = scal(2) + ! NONBONDED interactions + CALL section_vals_val_get(xtb_section, "DO_NONBONDED", & + l_val=qs_control%xtb_control%do_nonbonded) + CALL section_vals_get(nonbonded_section, explicit=explicit) + IF (explicit .AND. qs_control%xtb_control%do_nonbonded) THEN + CALL section_vals_get(genpot_section, explicit=explicit, n_repetition=ngp) + IF (explicit) THEN + CALL pair_potential_reallocate(qs_control%xtb_control%nonbonded, 1, ngp, gp=.TRUE.) + CALL read_gp_section(qs_control%xtb_control%nonbonded, genpot_section, 0) + END IF + END IF !nonbonded + ! XB_radius CALL section_vals_val_get(xtb_parameter, "XB_RADIUS", r_val=qs_control%xtb_control%xb_radius) ! Kab @@ -1557,6 +1573,8 @@ CONTAINS "xTB| Halogen interaction potential", xtb_control%xb_interaction WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") & "xTB| Halogen interaction potential cutoff radius", xtb_control%xb_radius + WRITE (UNIT=output_unit, FMT="(T2,A,T71,L10)") & + "xTB| Nonbonded interactions", xtb_control%do_nonbonded WRITE (UNIT=output_unit, FMT="(T2,A,T31,A50)") & "xTB| D3 Dispersion: Parameter file", ADJUSTR(TRIM(xtb_control%dispersion_parameter_file)) WRITE (UNIT=output_unit, FMT="(T2,A,T51,3F10.3)") & diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 45ceeda3ec..eeff2f924c 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -1140,6 +1140,15 @@ CONTAINS CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, & + name="sab_xtb_nonbond", & + description="Activates the printing of nonbonded interaction from (xTB) "// & + "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="soo_list", & description="Activates the printing of RI orbital-orbital "// & diff --git a/src/input_cp2k_tb.F b/src/input_cp2k_tb.F index ff31b045f5..8b815e0492 100644 --- a/src/input_cp2k_tb.F +++ b/src/input_cp2k_tb.F @@ -12,7 +12,6 @@ !> \author fawzi ! ************************************************************************************************** MODULE input_cp2k_tb -!qxtb citation USE bibliography, ONLY: Elstner1998,& Grimme2017,& Hu2007,& @@ -22,6 +21,7 @@ MODULE input_cp2k_tb USE input_constants, ONLY: dispersion_d3,& dispersion_uff,& slater + USE input_cp2k_mm, ONLY: create_GENPOT_section USE input_keyword_types, ONLY: keyword_create,& keyword_release,& keyword_type @@ -142,6 +142,10 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL create_xtb_nonbonded_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="DO_EWALD", & description="Use Ewald type method instead of direct sum for Coulomb interaction", & @@ -167,6 +171,13 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DO_NONBONDED", & + description="Controls the computation of real-space "// & + "(short-range) nonbonded interactions as correction to xTB.", & + usage="DO_NONBONDED 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="COULOMB_INTERACTION", & description="Use Coulomb interaction terms (electrostatics + TB3); for debug only", & usage="COULOMB_INTERACTION T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) @@ -395,7 +406,40 @@ CONTAINS CALL keyword_release(keyword) END SUBROUTINE create_xtb_parameter_section +! ************************************************************************************************** +!> \brief ... +!> \param section ... +! ************************************************************************************************** + SUBROUTINE create_xtb_nonbonded_section(section) + 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="NONBONDED", & + description="This section specifies the input parameters for NON-BONDED interactions.", & + n_keywords=1, n_subsections=0, repeats=.FALSE.) + NULLIFY (subsection) + + CALL create_GENPOT_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + CALL keyword_create(keyword, __LOCATION__, name="DX", & + description="Parameter used for computing the derivative with the Ridders' method.", & + usage="DX ", default_r_val=0.1_dp, unit_str="bohr") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ERROR_LIMIT", & + description="Checks that the error in computing the derivative is not larger than "// & + "the value set; in case error is larger a warning message is printed.", & + usage="ERROR_LIMIT ", default_r_val=1.0E-12_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_xtb_nonbonded_section ! ************************************************************************************************** !> \brief Creates the &ATOM_APRAMETER section !> \param section the section to create diff --git a/src/pair_potential.F b/src/pair_potential.F index 363ac0f144..f9cf0ed8f2 100644 --- a/src/pair_potential.F +++ b/src/pair_potential.F @@ -62,7 +62,8 @@ MODULE pair_potential INTEGER, PARAMETER, PRIVATE :: MAX_POINTS = 2000000 PUBLIC :: spline_nonbond_control, & - get_nonbond_storage + get_nonbond_storage, & + init_genpot CONTAINS diff --git a/src/qs_energy_types.F b/src/qs_energy_types.F index 7f48e39ec0..732c3aab7a 100644 --- a/src/qs_energy_types.F +++ b/src/qs_energy_types.F @@ -70,7 +70,9 @@ MODULE qs_energy_types kinetic, & !total kinetic energy [rk] s_square, & surf_dipole, & - embed_corr ! correction for embedding potential + embed_corr, & ! correction for embedding potential + xtb_xb_inter, & ! correction for halogen bonding within GFN1-xTB + xtb_nonbonded ! correction for nonbonded interactions within GFN1-xTB REAL(KIND=dp), DIMENSION(:), POINTER :: ddapc_restraint END TYPE qs_energy_type diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index ca00ee0bb4..87fcae1ca3 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -364,6 +364,7 @@ CONTAINS !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... +!> \param sab_xtb_nonbond ... !> \param sab_almo ... !> \param sab_kp ... !> \param particle_set ... @@ -506,7 +507,7 @@ CONTAINS SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, & dft_control, mos, mos_aux_fit, sab_orb, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_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_almo, sab_kp, particle_set, energy, force, & + sab_se, 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_ks_aux_fit, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, & @@ -547,7 +548,7 @@ CONTAINS 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_almo, sab_kp + sab_xtb_nonbond, sab_almo, sab_kp TYPE(particle_type), DIMENSION(:), OPTIONAL, & POINTER :: particle_set TYPE(qs_energy_type), OPTIONAL, POINTER :: energy @@ -830,6 +831,7 @@ CONTAINS sab_tbe=sab_tbe, & sab_core=sab_core, & sab_xb=sab_xb, & + sab_xtb_nonbond=sab_xtb_nonbond, & sab_almo=sab_almo, & sab_kp=sab_kp, & task_list=task_list, & diff --git a/src/qs_ks_types.F b/src/qs_ks_types.F index 6e2a48d328..f16034a1e3 100644 --- a/src/qs_ks_types.F +++ b/src/qs_ks_types.F @@ -116,6 +116,7 @@ MODULE qs_ks_types !> Ewald terms (DFTB) !> \param sab_core neighbor lists for the calculation of the core interactions !> \param sab_xb neighbor lists for the calculation of the XB interactions in xTB +!> \param sab_xtb_nonbond neighbor lists for the calculation of the nonbonded interactions in xTB !> \param sab_all neighbor lists for the calculation of the matrix element of !> non-symmetric operators !> \param sab_vdw: neighbor lists for the calculation of dispersion interactions @@ -194,6 +195,7 @@ MODULE qs_ks_types sab_tbe => Null(), & sab_core => Null(), & sab_xb => Null(), & + sab_xtb_nonbond => Null(), & sab_all => Null(), & sab_aux_fit => Null(), & sab_aux_fit_vs_orb => Null(), & @@ -304,6 +306,7 @@ CONTAINS !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... +!> \param sab_xtb_nonbond ... !> \param sab_vdw ... !> \param sab_scp ... !> \param sab_almo ... @@ -357,7 +360,7 @@ CONTAINS neighbor_list_id, & sab_orb, sab_aux_fit, sab_aux_fit_asymm, & sab_aux_fit_vs_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, & - sab_se, sab_tbe, sab_core, sab_xb, sab_vdw, sab_scp, sab_almo, sab_kp, & + sab_se, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & task_list, task_list_aux_fit, task_list_soft, & kpoints, do_kpoints, & atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, & @@ -385,7 +388,8 @@ CONTAINS INTEGER, OPTIONAL :: neighbor_list_id TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_aux_fit, & sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, & - sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, sab_vdw, sab_scp, sab_almo, sab_kp + sap_oce, sab_lrc, sab_se, 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_aux_fit, & task_list_soft TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints @@ -486,6 +490,7 @@ CONTAINS IF (PRESENT(sab_tbe)) sab_tbe => ks_env%sab_tbe 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 IF (PRESENT(sab_almo)) sab_almo => ks_env%sab_almo IF (PRESENT(sab_kp)) sab_kp => ks_env%sab_kp IF (PRESENT(dft_control)) dft_control => ks_env%dft_control @@ -583,6 +588,7 @@ CONTAINS !> \param sab_tbe ... !> \param sab_core ... !> \param sab_xb ... +!> \param sab_xtb_nonbond ... !> \param sab_vdw ... !> \param sab_scp ... !> \param sab_almo ... @@ -612,7 +618,7 @@ CONTAINS kpoints, & sab_orb, sab_aux_fit, sab_aux_fit_asymm, & sab_aux_fit_vs_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, & - sab_se, sab_tbe, sab_core, sab_xb, sab_vdw, sab_scp, sab_almo, sab_kp, & + sab_se, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, & task_list, task_list_aux_fit, task_list_soft, & subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, & para_env, blacs_env, admm_dm) @@ -635,7 +641,8 @@ CONTAINS TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_aux_fit, & sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, & - sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, sab_vdw, sab_scp, sab_almo, sab_kp + sap_oce, sab_lrc, sab_se, 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_aux_fit, & task_list_soft TYPE(qs_subsys_type), OPTIONAL, POINTER :: subsys @@ -713,6 +720,7 @@ CONTAINS IF (PRESENT(sab_tbe)) ks_env%sab_tbe => sab_tbe 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 IF (PRESENT(sab_almo)) ks_env%sab_almo => sab_almo IF (PRESENT(sab_kp)) ks_env%sab_kp => sab_kp IF (PRESENT(admm_dm)) ks_env%admm_dm => admm_dm @@ -860,6 +868,7 @@ CONTAINS CALL release_neighbor_list_sets(ks_env%sab_tbe) 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) CALL release_neighbor_list_sets(ks_env%sab_all) CALL release_neighbor_list_sets(ks_env%sab_aux_fit) CALL release_neighbor_list_sets(ks_env%sab_aux_fit_vs_orb) diff --git a/src/qs_neighbor_lists.F b/src/qs_neighbor_lists.F index dfc89f6633..0354ee5cf4 100644 --- a/src/qs_neighbor_lists.F +++ b/src/qs_neighbor_lists.F @@ -87,7 +87,8 @@ MODULE qs_neighbor_lists get_iterator_task, neighbor_list_iterate, neighbor_list_iterator_create, & neighbor_list_iterator_p_type, neighbor_list_iterator_release, neighbor_list_p_type, & neighbor_list_set_p_type, neighbor_list_set_type, release_neighbor_list_sets - USE string_utilities, ONLY: compress + USE string_utilities, ONLY: compress,& + uppercase USE subcell_types, ONLY: allocate_subcell,& deallocate_subcell,& give_ijk_subcell,& @@ -295,15 +296,16 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'build_qs_neighbor_lists' + CHARACTER(LEN=2) :: element_symbol, element_symbol2 CHARACTER(LEN=default_string_length) :: print_key_path - INTEGER :: handle, ikind, iw, maxatom, nkind, & - numshells, zat + INTEGER :: handle, ikind, ingp, iw, jkind, maxatom, & + ngp, nkind, numshells, zat LOGICAL :: all_potential_present, almo, dftb, do_hfx, dokp, explicit, gth_potential_present, & lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, paw_atom_present, rigpw, & sgp_potential_present, xtb LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, & - core_present, default_present, oce_present, orb_present, ppl_present, ppnl_present, & - ri_present, xb1_atom, xb2_atom + core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, orb_present, & + ppl_present, ppnl_present, ri_present, xb1_atom, xb2_atom REAL(dp) :: almo_rcov, almo_rvdw, rcut, roperator, & subcells REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, & @@ -325,8 +327,8 @@ CONTAINS TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, & sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_cn, sab_core, sab_gcp, sab_kp, & - sab_lrc, sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, sab_xb, sac_ae, sac_lri, sac_ppl, & - sap_oce, sap_ppnl, soa_list, soo_list + 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(paw_proj_set_type), POINTER :: paw_proj TYPE(qs_dftb_atom_type), POINTER :: dftb_atom @@ -356,6 +358,7 @@ CONTAINS NULLIFY (sab_tbe) NULLIFY (sab_core) NULLIFY (sab_xb) + NULLIFY (sab_xtb_nonbond) NULLIFY (sab_all) NULLIFY (sab_vdw) NULLIFY (sab_cn) @@ -401,6 +404,7 @@ CONTAINS sab_tbe=sab_tbe, & sab_core=sab_core, & sab_xb=sab_xb, & + sab_xtb_nonbond=sab_xtb_nonbond, & sab_scp=sab_scp, & sab_all=sab_all, & sab_almo=sab_almo, & @@ -801,6 +805,40 @@ CONTAINS CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, & symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xb") CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb) + CALL write_neighbor_lists(sab_xb, particle_set, cell, para_env, neighbor_list_section, & + "/SAB_XB", "sab_xb", "XB bonding") + + ! nonbonded interactions list + IF (dft_control%qs_control%xtb_control%do_nonbonded) THEN + ngp = SIZE(dft_control%qs_control%xtb_control%nonbonded%pot) + ALLOCATE (nonbond1_atom(nkind), nonbond2_atom(nkind)) + nonbond1_atom = .FALSE. + nonbond2_atom = .FALSE. + DO ingp = 1, ngp + DO ikind = 1, nkind + rcut = SQRT(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%rcutsq) + c_radius = rcut + CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=element_symbol) + CALL uppercase(element_symbol) + IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at1) == TRIM(element_symbol)) THEN + nonbond1_atom(ikind) = .TRUE. + DO jkind = 1, nkind + CALL get_atomic_kind(atomic_kind_set(jkind), element_symbol=element_symbol2) + CALL uppercase(element_symbol2) + IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at2) == TRIM(element_symbol2)) THEN + nonbond2_atom(jkind) = .TRUE. + ENDIF + ENDDO + END IF + END DO + CALL pair_radius_setup(nonbond1_atom, nonbond2_atom, c_radius, c_radius, pair_radius) + CALL build_neighbor_lists(sab_xtb_nonbond, particle_set, atom2d, cell, pair_radius, & + symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xtb_nonbond") + CALL set_ks_env(ks_env=ks_env, sab_xtb_nonbond=sab_xtb_nonbond) + CALL write_neighbor_lists(sab_xtb_nonbond, particle_set, cell, para_env, neighbor_list_section, & + "/SAB_XTB_NONBOND", "sab_xtb_nonbond", "XTB NONBONDED INTERACTIONS") + END DO + ENDIF END IF ! Build the neighbor lists for the vdW pair potential diff --git a/src/qs_scf_output.F b/src/qs_scf_output.F index 242fc63890..9d0b004825 100644 --- a/src/qs_scf_output.F +++ b/src/qs_scf_output.F @@ -372,6 +372,12 @@ CONTAINS "Electronic energy: ", energy%hartree, & "DFTB3 3rd order energy: ", energy%dftb3, & "Dispersion energy: ", energy%dispersion + IF (dft_control%qs_control%xtb_control%xb_interaction) & + WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & + "Correction for halogen bonding: ", energy%xtb_xb_inter + IF (dft_control%qs_control%xtb_control%do_nonbonded) & + WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & + "Correction for nonbonded interactions: ", energy%xtb_nonbonded IF (energy%efield /= 0.0_dp) & WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & "Electric field interaction energy: ", energy%efield diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F index 60f347f019..e4df89cf37 100644 --- a/src/xtb_matrices.F +++ b/src/xtb_matrices.F @@ -32,7 +32,8 @@ MODULE xtb_matrices USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_io_unit,& - cp_logger_type + cp_logger_type,& + cp_to_string USE cp_output_handling, ONLY: cp_p_file,& cp_print_key_finished_output,& cp_print_key_should_output,& @@ -42,16 +43,31 @@ MODULE xtb_matrices dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_finalize, dbcsr_get_block_p, & dbcsr_multiply, dbcsr_p_type, dbcsr_type USE efield_tb_methods, ONLY: efield_tb_matrix + USE ewald_environment_types, ONLY: ewald_env_get,& + ewald_environment_type + USE fparser, ONLY: evalfd,& + finalizef USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get - USE kinds, ONLY: dp + USE kinds, ONLY: default_string_length,& + dp USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type USE message_passing, ONLY: mp_minloc,& mp_sum USE mulliken, ONLY: ao_charges USE orbital_pointers, ONLY: ncoset + USE pair_potential, ONLY: init_genpot + USE pair_potential_types, ONLY: not_initialized,& + pair_potential_p_type,& + pair_potential_pp_create,& + pair_potential_pp_release,& + pair_potential_pp_type,& + pair_potential_single_clean,& + pair_potential_single_copy,& + pair_potential_single_type + USE pair_potential_util, ONLY: ener_pot USE particle_types, ONLY: particle_type USE qs_charge_mixing, ONLY: charge_mixing USE qs_condnum, ONLY: overlap_condnum @@ -83,6 +99,8 @@ MODULE xtb_matrices USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE qs_scf_types, ONLY: qs_scf_env_type + USE string_utilities, ONLY: compress,& + uppercase USE virial_methods, ONLY: virial_pair_force USE virial_types, ONLY: virial_type USE xtb_coulomb, ONLY: build_xtb_coulomb @@ -132,14 +150,13 @@ CONTAINS npgfb, nsgfa, nsgfb INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index - LOGICAL :: defined, diagblock, floating_a, found, & - ghost_a, norml1, norml2, omit_headers, & - use_arnoldi, use_virial, xb_inter + LOGICAL :: defined, diagblock, do_nonbonded, floating_a, found, ghost_a, norml1, norml2, & + omit_headers, use_arnoldi, use_virial, xb_inter LOGICAL, ALLOCATABLE, DIMENSION(:) :: floating, ghost - REAL(KIND=dp) :: alphaa, alphab, derepij, dfp, dhij, dr, drk, drx, ena, enb, erep, erepij, & - etaa, etab, exb, f0, f1, fen, fhua, fhub, fhud, foab, hij, k2sh, kab, kcnd, kcnp, kcns, & - kd, ken, kf, kg, kia, kjb, kp, ks, ksp, kx2, kxr, rcova, rcovab, rcovb, rrab, zneffa, & - zneffb + REAL(KIND=dp) :: alphaa, alphab, derepij, dfp, dhij, dr, drk, drx, ena, enb, enonbonded, & + erep, erepij, etaa, etab, exb, f0, f1, fen, fhua, fhub, fhud, foab, hij, k2sh, kab, kcnd, & + kcnp, kcns, kd, ken, kf, kg, kia, kjb, kp, ks, ksp, kx2, kxr, rcova, rcovab, rcovb, rrab, & + zneffa, zneffb REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers, kx REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, huckel, kcnlk, owork, rcab REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint @@ -168,7 +185,7 @@ CONTAINS TYPE(neighbor_list_iterator_p_type), & DIMENSION(:), POINTER :: nl_iterator TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_orb, sab_xb + POINTER :: sab_orb, sab_xb, sab_xtb_nonbond TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_dispersion_type), POINTER :: dispersion_env TYPE(qs_energy_type), POINTER :: energy @@ -180,6 +197,8 @@ CONTAINS TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b TYPE(xtb_control_type), POINTER :: xtb_control +!, na1 + CALL timeset(routineN, handle) NULLIFY (logger, virial, atprop) @@ -202,12 +221,16 @@ CONTAINS nkind = SIZE(atomic_kind_set) xtb_control => dft_control%qs_control%xtb_control xb_inter = xtb_control%xb_interaction + do_nonbonded = xtb_control%do_nonbonded nimg = dft_control%nimages nderivatives = 0 IF (calculate_forces) nderivatives = 1 IF (dft_control%tddfpt2_control%enabled) nderivatives = 1 maxder = ncoset(nderivatives) + IF (xb_inter) energy%xtb_xb_inter = 0.0_dp + IF (do_nonbonded) energy%xtb_nonbonded = 0.0_dp + ! global parameters (Table 2 from Ref.) ks = xtb_control%ks kp = xtb_control%kp @@ -842,8 +865,19 @@ CONTAINS CALL xb_interaction(exb, neighbor_atoms, atom_of_kind, kind_of, sab_xb, kx, kx2, rcab, & calculate_forces, use_virial, force, virial, atprop) END IF + + enonbonded = 0.0_dp + IF (do_nonbonded) THEN + ! nonbonded interactions + NULLIFY (sab_xtb_nonbond) + CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond) + CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, & + atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind) + ENDIF ! set repulsive energy - erep = erep + exb + erep = erep + exb + enonbonded + IF (xb_inter) energy%xtb_xb_inter = exb + IF (do_nonbonded) energy%xtb_nonbonded = enonbonded CALL mp_sum(erep, para_env%group) energy%repulsive = erep @@ -1148,6 +1182,12 @@ CONTAINS "Zeroth order Hamiltonian energy: ", energy%core, & "Charge fluctuation energy: ", energy%hartree, & "London dispersion energy: ", energy%dispersion + IF (dft_control%qs_control%xtb_control%xb_interaction) & + WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") & + "Correction for halogen bonding: ", energy%xtb_xb_inter + IF (dft_control%qs_control%xtb_control%do_nonbonded) & + WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") & + "Correction for nonbonded interactions: ", energy%xtb_nonbonded IF (ABS(energy%efield) > 1.e-12_dp) THEN WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") & "Electric field interaction energy: ", energy%efield @@ -1228,6 +1268,233 @@ CONTAINS END SUBROUTINE xb_neighbors ! ************************************************************************************************** +!> \brief Computes a correction for nonbonded interactions based on a generic potential +!> +!> \param enonbonded energy contribution +!> \param force ... +!> \param qs_env ... +!> \param xtb_control ... +!> \param sab_xtb_nonbond ... +!> \param atomic_kind_set ... +!> \param calculate_forces ... +!> \param use_virial ... +!> \param virial ... +!> \param atprop ... +!> \param atom_of_kind .. +!> \par History +!> 12.2018 JGH +!> \version 1.1 +! ************************************************************************************************** + SUBROUTINE nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, & + atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind) + REAL(dp), INTENT(INOUT) :: enonbonded + TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), & + POINTER :: force + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(xtb_control_type), POINTER :: xtb_control + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + INTENT(IN), POINTER :: sab_xtb_nonbond + TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), & + POINTER :: atomic_kind_set + LOGICAL, INTENT(IN) :: calculate_forces, use_virial + TYPE(virial_type), INTENT(IN), POINTER :: virial + TYPE(atprop_type), INTENT(IN), POINTER :: atprop + INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind + + CHARACTER(len=*), PARAMETER :: routineN = 'nonbonded_correction' + + CHARACTER(LEN=default_string_length) :: def_error, this_error + INTEGER :: atom_i, atom_j, handle, iatom, ikind, & + jatom, jkind, kk, ntype + INTEGER, DIMENSION(3) :: cell + LOGICAL :: do_ewald + REAL(KIND=dp) :: dedf, dr, dx, energy_cutoff, err, lerr, & + rcut + REAL(KIND=dp), DIMENSION(3) :: fij, rij + TYPE(ewald_environment_type), POINTER :: ewald_env + TYPE(neighbor_list_iterator_p_type), & + DIMENSION(:), POINTER :: nl_iterator + TYPE(pair_potential_p_type), POINTER :: nonbonded + TYPE(pair_potential_pp_type), POINTER :: potparm + TYPE(pair_potential_single_type), POINTER :: pot + TYPE(section_vals_type), POINTER :: nonbonded_section + +! + + CALL timeset(routineN, handle) + + NULLIFY (nonbonded) + NULLIFY (potparm) + NULLIFY (ewald_env) + nonbonded => xtb_control%nonbonded + do_ewald = xtb_control%do_ewald + CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env) + + ntype = SIZE(atomic_kind_set) + CALL pair_potential_pp_create(potparm, ntype) + !Assign input and potential info to potparm_nonbond + CALL force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald) + !Initialize genetic potential + CALL init_genpot(potparm, ntype) + + NULLIFY (pot) + enonbonded = 0._dp + energy_cutoff = 0._dp + + CALL neighbor_list_iterator_create(nl_iterator, sab_xtb_nonbond) + DO WHILE (neighbor_list_iterate(nl_iterator) == 0) + CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & + iatom=iatom, jatom=jatom, r=rij, cell=cell) + pot => potparm%pot(ikind, jkind)%pot + dr = SQRT(rij(1)**2 + rij(2)**2 + rij(3)**2) + rcut = SQRT(pot%rcutsq) + IF (dr <= rcut) THEN +! splines not implemented + enonbonded = enonbonded + ener_pot(pot, dr, energy_cutoff) + IF (atprop%energy) THEN + atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*ener_pot(pot, dr, energy_cutoff) + atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*ener_pot(pot, dr, energy_cutoff) + END IF + ENDIF + + IF (calculate_forces) THEN + + kk = SIZE(pot%type) + IF (kk /= 1) THEN + CALL cp_warn(__LOCATION__, "Generic potential with type > 1 not implemented.") + CPABORT("pot type") + ENDIF +! rmin and rmax and rcut + IF ((pot%set(kk)%rmin /= not_initialized) .AND. (dr < pot%set(kk)%rmin)) CYCLE + ! An upper boundary for the potential definition was defined + IF ((pot%set(kk)%rmax /= not_initialized) .AND. (dr >= pot%set(kk)%rmax)) CYCLE + ! If within limits let's compute the potential... + IF (dr <= rcut) THEN + + NULLIFY (nonbonded_section) + nonbonded_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%xTB%NONBONDED") + CALL section_vals_val_get(nonbonded_section, "DX", r_val=dx) + CALL section_vals_val_get(nonbonded_section, "ERROR_LIMIT", r_val=lerr) + + dedf = evalfd(pot%set(kk)%gp%myid, 1, pot%set(kk)%gp%values, dx, err) + IF (ABS(err) > lerr) THEN + WRITE (this_error, "(A,G12.6,A)") "(", err, ")" + WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")" + CALL compress(this_error, .TRUE.) + CALL compress(def_error, .TRUE.) + CALL cp_warn(__LOCATION__, & + 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// & + ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// & + TRIM(def_error)//' .') + END IF +! + atom_i = atom_of_kind(iatom) + atom_j = atom_of_kind(jatom) +! + fij(1:3) = dedf*rij(1:3)/pot%set(kk)%gp%values(1) + + force(ikind)%repulsive(:, atom_i) = force(ikind)%repulsive(:, atom_i) - fij(:) + force(jkind)%repulsive(:, atom_j) = force(jkind)%repulsive(:, atom_j) + fij(:) + + IF (use_virial) THEN + CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij) + IF (atprop%stress) THEN + CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij) + END IF + END IF + + END IF + END IF + NULLIFY (pot) + END DO + CALL neighbor_list_iterator_release(nl_iterator) + CALL finalizef() + + IF (ASSOCIATED(potparm)) THEN + CALL pair_potential_pp_release(potparm) + ENDIF + + CALL timestop(handle) +! + END SUBROUTINE nonbonded_correction +! ************************************************************************************************** +!> \brief ... +!> \param atomic_kind_set ... +!> \param nonbonded ... +!> \param potparm ... +!> \param ewald_env ... +!> \param do_ewald ... +! ************************************************************************************************** + SUBROUTINE force_field_pack_nonbond_pot_correction(atomic_kind_set, nonbonded, potparm, ewald_env, do_ewald) + +! routine based on force_field_pack_nonbond + TYPE(atomic_kind_type), DIMENSION(:), INTENT(IN), & + POINTER :: atomic_kind_set + TYPE(pair_potential_p_type), INTENT(IN), POINTER :: nonbonded + TYPE(pair_potential_pp_type), INTENT(INOUT), & + POINTER :: potparm + TYPE(ewald_environment_type), INTENT(IN), POINTER :: ewald_env + LOGICAL, INTENT(IN) :: do_ewald + + CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, & + name_atm_b, name_atm_b_local + INTEGER :: ikind, ingp, iw, jkind + LOGICAL :: found + REAL(KIND=dp) :: ewald_rcut + TYPE(atomic_kind_type), POINTER :: atomic_kind + TYPE(pair_potential_single_type), POINTER :: pot + + NULLIFY (pot) + DO ikind = 1, SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(ikind) + CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local) + DO jkind = ikind, SIZE(atomic_kind_set) + atomic_kind => atomic_kind_set(jkind) + CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local) + found = .FALSE. + name_atm_a = name_atm_a_local + name_atm_b = name_atm_b_local + CALL uppercase(name_atm_a) + CALL uppercase(name_atm_b) + pot => potparm%pot(ikind, jkind)%pot + IF (ASSOCIATED(nonbonded)) THEN + DO ingp = 1, SIZE(nonbonded%pot) + IF ((TRIM(nonbonded%pot(ingp)%pot%at1) == "*") .OR. & + (TRIM(nonbonded%pot(ingp)%pot%at2) == "*")) CYCLE + + IF (iw > 0) WRITE (iw, *) "TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & + " with ", TRIM(nonbonded%pot(ingp)%pot%at1), & + TRIM(nonbonded%pot(ingp)%pot%at2) + IF ((((name_atm_a) == (nonbonded%pot(ingp)%pot%at1)) .AND. & + ((name_atm_b) == (nonbonded%pot(ingp)%pot%at2))) .OR. & + (((name_atm_b) == (nonbonded%pot(ingp)%pot%at1)) .AND. & + ((name_atm_a) == (nonbonded%pot(ingp)%pot%at2)))) THEN + CALL pair_potential_single_copy(nonbonded%pot(ingp)%pot, pot) + ! multiple potential not implemented, simply overwriting + IF (found) & + CALL cp_warn(__LOCATION__, & + "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & + " and "//TRIM(name_atm_b)//" OVERWRITING! ") + END IF + IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) + found = .TRUE. + END DO + ENDIF + IF (.NOT. found) THEN + CALL pair_potential_single_clean(pot) + ENDIF + ENDDO !jkind + ENDDO !ikind + + ! Cutoff is defined always as the maximum between the FF and Ewald + IF (do_ewald) THEN + CALL ewald_env_get(ewald_env, rcut=ewald_rcut) + pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) + ENDIF + + END SUBROUTINE force_field_pack_nonbond_pot_correction +! ************************************************************************************************** !> \brief Computes the interaction term between Br/I/At and donor atoms !> !> \param exb ... diff --git a/tests/xTB/regtest-debug/TEST_FILES b/tests/xTB/regtest-debug/TEST_FILES index c1978147d7..1805a38c1c 100644 --- a/tests/xTB/regtest-debug/TEST_FILES +++ b/tests/xTB/regtest-debug/TEST_FILES @@ -9,4 +9,7 @@ ch2o_t03.inp 86 1.0E-12 0.11459198 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 +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.78101985886981 #EOF diff --git a/tests/xTB/regtest-debug/ch3br_atprop_nonbond.inp b/tests/xTB/regtest-debug/ch3br_atprop_nonbond.inp new file mode 100644 index 0000000000..dcb1f5fe57 --- /dev/null +++ b/tests/xTB/regtest-debug/ch3br_atprop_nonbond.inp @@ -0,0 +1,71 @@ +&FORCE_EVAL + STRESS_TENSOR ANALYTICAL + &PROPERTIES + &ATOMIC + ENERGY + PRESSURE + &END ATOMIC + &END PROPERTIES + METHOD Quickstep + &DFT + &QS + METHOD xTB + EPS_DEFAULT 1.0E-12 + &xTB + DO_EWALD T + CHECK_ATOMIC_CHARGES OFF + DO_NONBONDED T + &NONBONDED + &GENPOT + ATOMS Br Xe + FUNCTION Aparam*exp(-Bparam*r)-Cparam/r**8 + PARAMETERS Aparam Bparam Cparam + VALUES 49.507 1.441 6703.192 + VARIABLES r + &END GENPOT + DX 0.0001 + ERROR_LIMIT 0.0000001 + &END NONBONDED + &END XTB + &END QS + &SCF + SCF_GUESS RESTART + EPS_SCF 1.0E-8 + MAX_SCF 50 + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_S_INVERSE + &END + &OUTER_SCF + EPS_SCF 1.0E-8 + MAX_SCF 10 + &END + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 20.0 20.0 20.0 + &END CELL + &COORD + C -0.000001 -1.540826 -0.000000 + H 1.041396 -1.871752 -0.000000 + H -0.520681 -1.871772 0.901873 + H -0.520681 -1.871772 -0.901873 + Br -0.000001 0.424579 0.000000 + Xe -0.000000 3.0 -0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT CH3Br-Xe + PRINT_LEVEL LOW + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 3 + TIMESTEP 0.48 + TEMPERATURE 330.0 + &END MD +&END MOTION diff --git a/tests/xTB/regtest-debug/ch3br_nonbond.inp b/tests/xTB/regtest-debug/ch3br_nonbond.inp new file mode 100644 index 0000000000..85c69ab0bd --- /dev/null +++ b/tests/xTB/regtest-debug/ch3br_nonbond.inp @@ -0,0 +1,56 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD xTB + EPS_DEFAULT 1.0E-10 + EPS_PGF_ORB 1.0E-20 + &xTB + DO_EWALD F + CHECK_ATOMIC_CHARGES OFF + DO_NONBONDED T + &NONBONDED + &GENPOT + ATOMS Br Xe + FUNCTION Aparam*exp(-Bparam*r)-Cparam/r**8 + PARAMETERS Aparam Bparam Cparam + VALUES 49.507 1.441 6703.192 + VARIABLES r + &END GENPOT + &END NONBONDED + &END XTB + &END QS + &SCF + SCF_GUESS RESTART + EPS_SCF 1.0E-5 + MAX_SCF 100 + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_S_INVERSE + &END + &OUTER_SCF + EPS_SCF 1.0E-5 + MAX_SCF 900 + &END + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 10 10 10 + PERIODIC NONE + &END CELL + &COORD + C -0.000001 -1.540826 -0.000000 + H 1.041396 -1.871752 -0.000000 + H -0.520681 -1.871772 0.901873 + H -0.520681 -1.871772 -0.901873 + Br -0.000001 0.424579 0.000000 + Xe -0.000000 3.0 -0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT CH3Br-Xe + PRINT_LEVEL LOW + RUN_TYPE ENERGY +&END GLOBAL diff --git a/tests/xTB/regtest-debug/ch3br_nonbond_2.inp b/tests/xTB/regtest-debug/ch3br_nonbond_2.inp new file mode 100644 index 0000000000..f383834bbb --- /dev/null +++ b/tests/xTB/regtest-debug/ch3br_nonbond_2.inp @@ -0,0 +1,70 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD xTB + &xTB + DO_EWALD F + CHECK_ATOMIC_CHARGES OFF + DO_NONBONDED T + &NONBONDED + &GENPOT + ATOMS Br Xe + FUNCTION Aparam*exp(-Bparam*r)-Cparam/r**8 + PARAMETERS Aparam Bparam Cparam + VALUES 49.507 1.441 6703.192 + VARIABLES r + &END GENPOT + &END NONBONDED + &END XTB + &END QS + &SCF + SCF_GUESS RESTART + EPS_SCF 1.0E-5 + MAX_SCF 100 + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_S_INVERSE + &END + &OUTER_SCF + EPS_SCF 1.0E-5 + MAX_SCF 900 + &END + &END SCF + &EFIELD + &END + &PRINT + &MOMENTS ON + PERIODIC .FALSE. + REFERENCE COM + &END + &END + &END DFT + &SUBSYS + &CELL + ABC 10 10 10 + PERIODIC NONE + &END CELL + &COORD + C -0.000001 -1.540826 -0.000000 + H 1.041396 -1.871752 -0.000000 + H -0.520681 -1.871772 0.901873 + H -0.520681 -1.871772 -0.901873 + Br -0.000001 0.424579 0.000000 + Xe -0.000000 3.0 -0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT CH3Br-Xe + PRINT_LEVEL LOW + RUN_TYPE DEBUG +&END GLOBAL +&DEBUG + DEBUG_FORCES F + DEBUG_STRESS_TENSOR F + DEBUG_DIPOLE T + DEBUG_POLARIZABILITY F + STOP_ON_MISMATCH T + DE 0.0002 +&END diff --git a/tests/xTB/regtest-debug/h2o-atprop_nonbonded.inp b/tests/xTB/regtest-debug/h2o-atprop_nonbonded.inp new file mode 100644 index 0000000000..16f01bb594 --- /dev/null +++ b/tests/xTB/regtest-debug/h2o-atprop_nonbonded.inp @@ -0,0 +1,176 @@ +&FORCE_EVAL + STRESS_TENSOR ANALYTICAL + &PROPERTIES + &ATOMIC + ENERGY + PRESSURE + &END ATOMIC + &END PROPERTIES + &DFT + &QS + METHOD xTB + &xTB + DO_EWALD T + COULOMB_INTERACTION T + TB3_INTERACTION T + DO_NONBONDED T + &NONBONDED + &GENPOT + ATOMS O O + FUNCTION 4*epsilon*((sigma/(r+sigma))**12 - (sigma/(r+sigma))**6) + PARAMETERS epsilon sigma + VALUES 0.00380879917604 2.31204899506 + VARIABLES r + RCUT 4.0 + RMIN 1.0 + RMAX 4.0 + &END GENPOT + DX 0.001 + ERROR_LIMIT 0.0000000001 + &END NONBONDED + &END + &END QS + &SCF + SCF_GUESS MOPAC + EPS_SCF 1.E-5 + &MIXING + METHOD DIRECT_P_MIXING + ALPHA 0.75 + &END + MAX_SCF 20 + &END SCF + &POISSON + &EWALD + O_SPLINE 6 + &END EWALD + &END POISSON + &END DFT + &SUBSYS + &CELL + ABC 9.865 9.865 9.865 + &END CELL + &COORD + O -4.583 5.333 1.560 H2O + H -3.777 5.331 0.943 H2O + H -5.081 4.589 1.176 H2O + O -0.083 4.218 0.070 H2O + H -0.431 3.397 0.609 H2O + H 0.377 3.756 -0.688 H2O + O -1.488 2.692 4.125 H2O + H -2.465 2.433 3.916 H2O + H -1.268 2.145 4.952 H2O + O -2.461 -2.548 -6.136 H2O + H -1.892 -2.241 -6.921 H2O + H -1.970 -3.321 -5.773 H2O + O 4.032 0.161 2.183 H2O + H 4.272 -0.052 1.232 H2O + H 4.044 -0.760 2.641 H2O + O 2.950 -3.497 -1.006 H2O + H 2.599 -3.901 -0.129 H2O + H 3.193 -4.283 -1.533 H2O + O -2.890 -4.797 -2.735 H2O + H -2.810 -5.706 -2.297 H2O + H -2.437 -4.128 -2.039 H2O + O -0.553 0.922 -3.731 H2O + H -0.163 1.552 -3.085 H2O + H -1.376 0.544 -3.469 H2O + O 4.179 4.017 4.278 H2O + H 3.275 3.832 3.876 H2O + H 4.658 4.492 3.572 H2O + O 5.739 1.425 3.944 H2O + H 5.125 2.066 4.346 H2O + H 5.173 1.181 3.097 H2O + O 0.988 -0.760 -5.445 H2O + H 1.640 -1.372 -4.989 H2O + H 0.546 -0.220 -4.762 H2O + O -0.748 1.985 1.249 H2O + H -0.001 1.490 1.540 H2O + H -1.160 2.255 2.109 H2O + O 4.127 -0.234 -3.149 H2O + H 5.022 -0.436 -3.428 H2O + H 3.540 -0.918 -3.601 H2O + O -2.473 2.768 -1.395 H2O + H -1.533 2.719 -1.214 H2O + H -2.702 1.808 -1.479 H2O + O -0.124 -2.116 2.404 H2O + H 0.612 -2.593 2.010 H2O + H 0.265 -1.498 3.089 H2O + O 0.728 2.823 -2.190 H2O + H 0.646 3.694 -2.685 H2O + H 1.688 2.705 -1.947 H2O + O 4.256 -5.427 -2.644 H2O + H 5.222 -5.046 -2.479 H2O + H 4.174 -5.628 -3.593 H2O + O -3.178 -0.508 -4.227 H2O + H -2.762 -1.221 -4.818 H2O + H -3.603 0.073 -4.956 H2O + O -1.449 5.300 -4.805 H2O + H -1.397 4.470 -5.317 H2O + H -2.102 5.091 -4.067 H2O + O 3.354 2.192 -1.755 H2O + H 3.407 1.433 -2.405 H2O + H 3.971 2.958 -2.196 H2O + O 1.773 -4.018 1.769 H2O + H 1.121 -4.532 1.201 H2O + H 1.975 -4.529 2.618 H2O + O 1.526 1.384 2.712 H2O + H 2.317 1.070 2.251 H2O + H 1.353 0.657 3.364 H2O + O 2.711 -2.398 -4.253 H2O + H 2.202 -3.257 -4.120 H2O + H 3.305 -2.610 -5.099 H2O + O 6.933 0.093 -1.393 H2O + H 6.160 -0.137 -0.795 H2O + H 6.748 -0.394 -2.229 H2O + O -5.605 -2.549 3.151 H2O + H -4.756 -2.503 3.616 H2O + H -5.473 -3.187 2.378 H2O + O 0.821 -4.406 6.516 H2O + H 0.847 -3.675 7.225 H2O + H -0.014 -4.240 5.988 H2O + O 1.577 3.933 3.762 H2O + H 1.221 2.975 3.640 H2O + H 1.367 4.126 4.659 H2O + O -2.111 -3.741 -0.219 H2O + H -1.378 -4.425 -0.036 H2O + H -1.825 -2.775 0.003 H2O + O 0.926 -1.961 -2.063 H2O + H 0.149 -1.821 -1.402 H2O + H 1.725 -2.303 -1.536 H2O + O 4.531 -1.030 -0.547 H2O + H 4.290 -1.980 -0.581 H2O + H 4.292 -0.597 -1.390 H2O + O -0.740 -1.262 -0.029 H2O + H -1.272 -0.422 -0.099 H2O + H -0.403 -1.349 0.873 H2O + O 3.655 3.021 0.988 H2O + H 2.706 3.053 1.282 H2O + H 3.542 2.615 0.020 H2O + &END COORD + &END SUBSYS +&END FORCE_EVAL +&GLOBAL + PROJECT h2o + RUN_TYPE MD + PRINT_LEVEL LOW +&END GLOBAL +&MOTION + &MD + ENSEMBLE NPT_I + STEPS 3 + TIMESTEP 0.48 + TEMPERATURE 330.0 + &THERMOSTAT + &NOSE + LENGTH 3 + YOSHIDA 3 + TIMECON [wavenumber_t] 2000 + MTS 2 + &END NOSE + &END + &BAROSTAT + PRESSURE 1.0 + TIMECON 300 + &END BAROSTAT + &END MD +&END MOTION