xTB short range Coulomb: use special neighbor list (#2287)

this solves the electronic energy dependence on EPS_DEFAULT
cp2k · Sep 16, 2022 · 70ce0ff · 70ce0ff
1 parent fd7edf4
commit 70ce0ff
Show file tree

Hide file tree

Showing 23 changed files with 144 additions and 94 deletions.
diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F
@@ -805,6 +805,15 @@ SUBROUTINE create_print_dft_section(section)
       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
@@ -420,13 +420,13 @@ SUBROUTINE create_xtb_parameter_section(section)
       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
@@ -992,11 +992,11 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
             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
@@ -361,6 +361,7 @@ MODULE qs_environment_types
 !> \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 @@ MODULE qs_environment_types
 !> \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 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
          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 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
                       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
@@ -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 @@ END SUBROUTINE qs_ks_env_create
 !> \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 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
                          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 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
       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 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
       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 @@ END SUBROUTINE get_ks_env
 !> \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 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
                          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 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
       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 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
       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 @@ SUBROUTINE qs_ks_release(ks_env)
             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 @@ SUBROUTINE qs_ks_part_release(ks_env)
             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
@@ -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 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
       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 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
       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 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
       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 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
                       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 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
                                       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
@@ -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 @@ SUBROUTINE stda_init_param(qs_env, stda_kernel, stda_control)
       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 @@ SUBROUTINE stda_init_param(qs_env, stda_kernel, stda_control)
 
       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 @@ SUBROUTINE stda_init_param(qs_env, stda_kernel, stda_control)
          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
@@ -180,8 +180,9 @@ SUBROUTINE setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim)
 
       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