Bug fix for CN forces and first part of XB interaction term

cp2k · Mar 19, 2019 · f489f15 · f489f15
1 parent 841ab0f
commit f489f15
Show file tree

Hide file tree

Showing 9 changed files with 137 additions and 42 deletions.
diff --git a/src/cp_control_types.F b/src/cp_control_types.F
@@ -109,6 +109,8 @@ MODULE cp_control_types
       REAL(KIND=dp)                        :: ken
       REAL(KIND=dp)                        :: kxr, kx2
       !
+      REAL(KIND=dp)                        :: xb_radius
+      !
       CHARACTER(LEN=default_string_length), &
          DIMENSION(:, :), POINTER          :: kab_param
    END TYPE xtb_control_type

diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F
@@ -1084,6 +1084,8 @@ SUBROUTINE read_qs_section(qs_control, qs_section)
          qs_control%xtb_control%kx2 = scal(2)
          ! Kab
          CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", n_rep_val=n_rep)
+         ! XB_radius
+         CALL section_vals_val_get(xtb_parameter, "XB_RADIUS", r_val=qs_control%xtb_control%xb_radius)
          IF (n_rep > 0) THEN
             ALLOCATE (qs_control%xtb_control%kab_param(3, n_rep))
             DO j = 1, n_rep

diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F
@@ -1066,6 +1066,15 @@ SUBROUTINE create_print_dft_section(section)
       CALL section_add_keyword(print_key, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword=keyword, &
+                          name="sab_xb", &
+                          description="Activates the printing of XB 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=keyword, &
                           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
@@ -367,6 +367,13 @@ SUBROUTINE create_xtb_parameter_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, name="XB_RADIUS", &
+                          description="Specifies the radius [Bohr] of the XB pair interaction in xTB.", &
+                          usage="XB_RADIUS 20.0 ", repeats=.FALSE., &
+                          n_var=1, default_r_val=20.0_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
    END SUBROUTINE create_xtb_parameter_section
 
 ! **************************************************************************************************

diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F
@@ -346,6 +346,7 @@ MODULE qs_environment_types
 !> \param sab_se ...
 !> \param sab_tbe ...
 !> \param sab_core ...
+!> \param sab_xb ...
 !> \param sab_almo ...
 !> \param sab_kp ...
 !> \param particle_set ...
@@ -483,7 +484,7 @@ MODULE qs_environment_types
    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_almo, sab_kp, particle_set, energy, force, &
+                         sab_se, sab_tbe, sab_core, sab_xb, 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, &
@@ -522,11 +523,9 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
          OPTIONAL, POINTER                               :: sab_orb, sab_aux_fit, sab_aux_fit_asymm, &
                                                             sab_aux_fit_vs_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_almo, &
-                                                            sab_kp
+      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
       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
          POINTER                                         :: particle_set
       TYPE(qs_energy_type), OPTIONAL, POINTER            :: energy
@@ -796,6 +795,7 @@ SUBROUTINE get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, ce
                       sab_lrc=sab_lrc, &
                       sab_tbe=sab_tbe, &
                       sab_core=sab_core, &
+                      sab_xb=sab_xb, &
                       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
@@ -113,6 +113,7 @@ MODULE qs_ks_types
 !> \param sac_tbe: neighbor lists for the calculation of the tight binding
 !>        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_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
@@ -190,6 +191,7 @@ MODULE qs_ks_types
                                                                sab_scp => Null(), &
                                                                sab_tbe => Null(), &
                                                                sab_core => Null(), &
+                                                               sab_xb => Null(), &
                                                                sab_all => Null(), &
                                                                sab_aux_fit => Null(), &
                                                                sab_aux_fit_vs_orb => Null(), &
@@ -302,6 +304,7 @@ END SUBROUTINE qs_ks_env_create
 !> \param sab_se ...
 !> \param sab_tbe ...
 !> \param sab_core ...
+!> \param sab_xb ...
 !> \param sab_vdw ...
 !> \param sab_scp ...
 !> \param sab_almo ...
@@ -355,7 +358,7 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
                          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_vdw, sab_scp, sab_almo, sab_kp, &
+                         sab_se, sab_tbe, sab_core, sab_xb, 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, &
@@ -383,7 +386,7 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
       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_vdw, sab_scp, sab_almo, sab_kp
+         sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, 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
@@ -483,6 +486,7 @@ SUBROUTINE get_ks_env(ks_env, v_hartree_rspace, &
       IF (PRESENT(sab_lrc)) sab_lrc => ks_env%sab_lrc
       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_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
@@ -579,6 +583,7 @@ END SUBROUTINE get_ks_env
 !> \param sab_se ...
 !> \param sab_tbe ...
 !> \param sab_core ...
+!> \param sab_xb ...
 !> \param sab_vdw ...
 !> \param sab_scp ...
 !> \param sab_almo ...
@@ -608,7 +613,7 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
                          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_vdw, sab_scp, sab_almo, sab_kp, &
+                         sab_se, sab_tbe, sab_core, sab_xb, 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)
@@ -631,7 +636,7 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
       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_vdw, sab_scp, sab_almo, sab_kp
+         sap_oce, sab_lrc, sab_se, sab_tbe, sab_core, sab_xb, 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
@@ -708,6 +713,7 @@ SUBROUTINE set_ks_env(ks_env, v_hartree_rspace, &
       IF (PRESENT(sab_lrc)) ks_env%sab_lrc => sab_lrc
       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_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
@@ -856,6 +862,7 @@ SUBROUTINE qs_ks_release(ks_env)
             CALL release_sab(ks_env%sab_scp)
             CALL release_sab(ks_env%sab_tbe)
             CALL release_sab(ks_env%sab_core)
+            CALL release_sab(ks_env%sab_xb)
             CALL release_sab(ks_env%sab_all)
             CALL release_sab(ks_env%sab_aux_fit)
             CALL release_sab(ks_env%sab_aux_fit_vs_orb)

diff --git a/src/qs_neighbor_lists.F b/src/qs_neighbor_lists.F
@@ -302,7 +302,7 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
          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
+         ri_present, xb1_atom, xb2_atom
       REAL(dp)                                           :: almo_rcov, almo_rvdw, alpha, roperator, &
                                                             subcells
       REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, &
@@ -324,8 +324,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_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, sab_cn, sab_core, sab_kp, sab_lrc, &
-         sab_orb, sab_scp, sab_se, sab_tbe, sab_vdw, sac_ae, sac_lri, sac_ppl, sap_oce, sap_ppnl, &
-         soa_list, soo_list
+         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
       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
       TYPE(paw_proj_set_type), POINTER                   :: paw_proj
       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_atom
@@ -353,6 +353,7 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
       NULLIFY (sab_lrc)
       NULLIFY (sab_tbe)
       NULLIFY (sab_core)
+      NULLIFY (sab_xb)
       NULLIFY (sab_all)
       NULLIFY (sab_vdw)
       NULLIFY (sab_cn)
@@ -397,6 +398,7 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
                       sab_lrc=sab_lrc, &
                       sab_tbe=sab_tbe, &
                       sab_core=sab_core, &
+                      sab_xb=sab_xb, &
                       sab_scp=sab_scp, &
                       sab_all=sab_all, &
                       sab_almo=sab_almo, &
@@ -766,7 +768,6 @@ SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_sectio
          END IF
       END IF
 
-!qxtb
       IF (xtb) THEN
          ! Build the neighbor lists for the DFTB Ewald methods
          IF (dft_control%qs_control%xtb_control%do_ewald) THEN
@@ -778,8 +779,27 @@ 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
+         ! XB list
+         ALLOCATE (xb1_atom(nkind), xb2_atom(nkind))
+         c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius
+         DO ikind = 1, nkind
+            CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
+            IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
+               xb1_atom(ikind) = .TRUE.
+            ELSE
+               xb1_atom(ikind) = .FALSE.
+            END IF
+            IF (zat == 7 .OR. zat == 8 .OR. zat == 15 .OR. zat == 16) THEN
+               xb2_atom(ikind) = .TRUE.
+            ELSE
+               xb2_atom(ikind) = .FALSE.
+            END IF
+         END DO
+         CALL pair_radius_setup(xb1_atom, xb2_atom, c_radius, c_radius, pair_radius)
+         CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, &
+                                   subcells=subcells, operator_type="PP", nlname="sab_xb")
+         CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb)
       END IF
-!qxtb
 
       ! Build the neighbor lists for the vdW pair potential
       CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)