Add D3_EXCLUDE_KIND_PAIR keyword (#956)

* Add D3_EXCLUDE_KIND_PAIR keyword

The D3_EXCLUDE_KIND_PAIR keyword allows to exclude specific
interactions from the D3 dispersion correction.

* Fix the description of the keyword

src/input_cp2k_xc.F

@@ -1038,6 +1038,14 @@ SUBROUTINE create_vdw_potential_section(section)
       CALL section_add_keyword(subsection, keyword)
       CALL keyword_release(keyword)
 
+      ! Ignore selected pair interactins
+      CALL keyword_create(keyword, __LOCATION__, name="D3_EXCLUDE_KIND_PAIR", &
+                          description="Specifies the atomic kinds for interactions excluded from the DFT-D3 calculation.", &
+                          usage="D3_EXCLUDE_KIND_PAIR kind1 kind2 ", repeats=.TRUE., &
+                          n_var=2, type_of_var=integer_t)
+      CALL section_add_keyword(subsection, keyword)
+      CALL keyword_release(keyword)
+
       ! Set coordination numbers by atom kinds
       CALL keyword_create(keyword, __LOCATION__, name="KIND_COORDINATION_NUMBERS", &
                           description="Specifies the coordination number for a kind for the C9 term in DFT-D3.", &

src/qs_dispersion_pairpot.F

@@ -328,6 +328,16 @@ SUBROUTINE qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_e
                      disp%defined = .FALSE.
                   END DO
                END IF
+               CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", n_rep_val=n_rep)
+               dispersion_env%nd3_exclude_pair = n_rep
+               IF (n_rep > 0) THEN
+                  ALLOCATE (dispersion_env%d3_exclude_pair(n_rep, 2))
+                  DO i = 1, n_rep
+                     CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", i_rep_val=i, &
+                                               i_vals=exlist)
+                     dispersion_env%d3_exclude_pair(i, :) = exlist
+                  END DO
+               END IF
             END IF
 
          END SELECT
@@ -1199,8 +1209,8 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomnumber, kind_of
       INTEGER, DIMENSION(3)                              :: cell_b, cell_c, ncell, periodic
       INTEGER, DIMENSION(:), POINTER                     :: atom_list
-      LOGICAL :: atenergy, atstress, debug, debugall, debugx, domol, floating_a, floating_b, &
-         floating_c, ghost_a, ghost_b, ghost_c, is000, use_virial
+      LOGICAL :: atenergy, atstress, debug, debugall, debugx, domol, exclude_pair, floating_a, &
+         floating_b, floating_c, ghost_a, ghost_b, ghost_c, is000, use_virial
       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: dodisp, floating, ghost
       REAL(KIND=dp) :: a1, a2, alp6, alp8, ang, c6, c8, c9, cc6ab, cc6bc, cc6ca, cnum, dc6a, dc6b, &
          dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, dd, de6, de8, de91, &
@@ -1550,6 +1560,10 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
                   atener(jatom) = atener(jatom) - 0.5_dp*xp*fdmp*fac
                END IF
             ELSE IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
+               IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
+                  IF (exclude_pair) CYCLE
+               END IF
                ! C6 coefficient
                CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
                           cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
@@ -1719,6 +1733,11 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
                               IF (r2ca >= rc2) CYCLE
                               r2abc = r2ab*r2bc*r2ca
                               IF (r2abc <= 0.001_dp) CYCLE
+                              IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                                 CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
+                                                           kkind, exclude_pair)
+                                 IF (exclude_pair) CYCLE
+                              END IF
                               ! this is an empirical scaling
                               IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
                                  kkind = kind_of(katom)
@@ -1962,6 +1981,10 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
                   CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
                   CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
                   IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
+                  IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                     CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
+                     IF (exclude_pair) CYCLE
+                  END IF
                   CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
                              cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
                   elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
@@ -1972,6 +1995,10 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
                         CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
                         CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
                         IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
+                        IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                           CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, kkind, exclude_pair)
+                           IF (exclude_pair) CYCLE
+                        END IF
                         CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
                                    cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
                         CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
@@ -3382,6 +3409,7 @@ SUBROUTINE d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating,
 
       INTEGER                                            :: handle, iatom, ikind, jatom, jkind, &
                                                             mepos, natom, ni, nj, num_pe, za, zb
+      LOGICAL                                            :: exclude_pair
       REAL(KIND=dp)                                      :: cnab, dcnab, eps_cn, rcc, rcova, rcovb
       REAL(KIND=dp), DIMENSION(3)                        :: rij
       TYPE(neighbor_list_iterator_p_type), &
@@ -3416,7 +3444,7 @@ SUBROUTINE d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating,
 !$OMP          PRIVATE( mepos                                         &
 !$OMP                 , ikind, jkind, iatom, jatom, rij, rcc          &
 !$OMP                 , za, zb, rcova, rcovb, cnab, dcnab             &
-!$OMP                 , ni, nj                                        &
+!$OMP                 , ni, nj, exclude_pair                          &
 !$OMP                 )
 
 !$OMP SINGLE
@@ -3435,6 +3463,10 @@ SUBROUTINE d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating,
 
          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
          IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) CYCLE
+         IF (dispersion_env%nd3_exclude_pair > 0) THEN
+            CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, exclude=exclude_pair)
+            IF (exclude_pair) CYCLE
+         END IF
 
          rcc = SQRT(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
          IF (rcc > 1.e-6_dp) THEN
@@ -3502,4 +3534,47 @@ END SUBROUTINE d3_cnumber
 
 ! **************************************************************************************************
 
+! **************************************************************************************************
+!> \brief ...
+!> \param exclude_list List of kind pairs to exclude
+!> \param za first kind
+!> \param zb second kind
+!> \param zc third kind in case of the three-body term
+!> \param exclude whether to exclude the interaction or not
+! **************************************************************************************************
+   SUBROUTINE exclude_d3_kind_pair(exclude_list, za, zb, zc, exclude)
+
+      INTEGER, DIMENSION(:, :), INTENT(IN)               :: exclude_list
+      INTEGER, INTENT(in)                                :: za, zb
+      INTEGER, INTENT(in), OPTIONAL                      :: zc
+      LOGICAL, INTENT(out)                               :: exclude
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'exclude_d3_kind_pair', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle, i
+
+      CALL timeset(routineN, handle)
+      exclude = .FALSE.
+      DO i = 1, SIZE(exclude_list, 1)
+         IF (exclude_list(i, 1) == za .AND. exclude_list(i, 2) == zb .OR. &
+             exclude_list(i, 1) == zb .AND. exclude_list(i, 2) == za) THEN
+            exclude = .TRUE.
+            EXIT
+         END IF
+         IF (PRESENT(zc)) THEN
+            IF (exclude_list(i, 1) == za .AND. exclude_list(i, 2) == zc .OR. &
+                exclude_list(i, 1) == zc .AND. exclude_list(i, 2) == za .OR. &
+                exclude_list(i, 1) == zb .AND. exclude_list(i, 2) == zc .OR. &
+                exclude_list(i, 1) == zc .AND. exclude_list(i, 2) == zb) THEN
+               exclude = .TRUE.
+               EXIT
+            END IF
+         END IF
+      END DO
+
+      CALL timestop(handle)
+
+   END SUBROUTINE exclude_d3_kind_pair
+
 END MODULE qs_dispersion_pairpot

src/qs_dispersion_types.F

@@ -90,6 +90,8 @@ MODULE qs_dispersion_types
       !! kernel matrix at each of the q points.  Stored as
       !! d2phi_dk2(k_point, q1_value, q2_value)
       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: d2y_dx2 !! 2nd derivatives of q_mesh for interpolation
+      INTEGER, DIMENSION(:, :), POINTER  :: d3_exclude_pair
+      INTEGER  :: nd3_exclude_pair
    END TYPE qs_dispersion_type
 
    TYPE qs_atom_dispersion_type
@@ -162,6 +164,9 @@ SUBROUTINE qs_dispersion_release(dispersion_env)
          IF (ASSOCIATED(dispersion_env%d2y_dx2)) THEN
             DEALLOCATE (dispersion_env%d2y_dx2)
          END IF
+         IF (ASSOCIATED(dispersion_env%d3_exclude_pair)) THEN
+            DEALLOCATE (dispersion_env%d3_exclude_pair)
+         END IF
          ! neighborlists
          CALL release_neighbor_list_sets(dispersion_env%sab_vdw)
          CALL release_neighbor_list_sets(dispersion_env%sab_cn)

src/qs_dispersion_utils.F

@@ -75,8 +75,10 @@ SUBROUTINE qs_dispersion_env_set(dispersion_env, xc_section)
       dispersion_env%lrc = .FALSE.
       dispersion_env%srb = .FALSE.
       dispersion_env%verbose = .FALSE.
+      dispersion_env%nd3_exclude_pair = 0
       NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
-               dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist)
+               dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
+               dispersion_env%d3_exclude_pair)
       NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                dispersion_env%d2y_dx2)
       NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
@@ -246,6 +248,12 @@ SUBROUTINE qs_write_dispersion(qs_env, dispersion_env, ounit)
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'sr6 Scaling Factor:',T73,F8.4)") dispersion_env%sr6
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn
+               IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                  DO i = 1, dispersion_env%nd3_exclude_pair
+                     WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Excluded Pairs: ',T76,I2,' ',I2)") &
+                        dispersion_env%d3_exclude_pair(i, :)
+                  END DO
+               END IF
             ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D3 (Version 3.1)')")
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Potential Form: S. Grimme et al, JCP 132: 154104 (2010)')")
@@ -256,6 +264,12 @@ SUBROUTINE qs_write_dispersion(qs_env, dispersion_env, ounit)
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'a2 Damping Factor:',T73,F8.4)") dispersion_env%a2
                WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn
+               IF (dispersion_env%nd3_exclude_pair > 0) THEN
+                  DO i = 1, dispersion_env%nd3_exclude_pair
+                     WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Excluded Kind Pairs: ',T76,I2,' ',I2)") &
+                        dispersion_env%d3_exclude_pair(i, :)
+                  END DO
+               END IF
             END IF
          ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
             WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T61,'Non-local Functional')")

src/qs_environment.F

@@ -1169,7 +1169,8 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          dispersion_env%srb = .FALSE.
          dispersion_env%verbose = .FALSE.
          NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
-                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist)
+                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
+                  dispersion_env%d3_exclude_pair)
          NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                   dispersion_env%d2y_dx2, dispersion_env%dftd_section)
          NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
@@ -1185,6 +1186,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
             dispersion_env%rc_disp = dftb_control%rcdisp
             dispersion_env%exp_pre = 0._dp
             dispersion_env%scaling = 0._dp
+            dispersion_env%nd3_exclude_pair = 0
             dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
             CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
          ELSE
@@ -1200,7 +1202,8 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          dispersion_env%srb = .FALSE.
          dispersion_env%verbose = .FALSE.
          NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
-                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist)
+                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
+                  dispersion_env%d3_exclude_pair)
          NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                   dispersion_env%d2y_dx2, dispersion_env%dftd_section)
          NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
@@ -1216,6 +1219,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          dispersion_env%rc_disp = xtb_control%rcdisp
          dispersion_env%exp_pre = 0._dp
          dispersion_env%scaling = 0._dp
+         dispersion_env%nd3_exclude_pair = 0
          dispersion_env%parameter_file_name = xtb_control%dispersion_parameter_file
          CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
          CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
@@ -1228,7 +1232,8 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
          dispersion_env%srb = .FALSE.
          dispersion_env%verbose = .FALSE.
          NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
-                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist)
+                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
+                  dispersion_env%d3_exclude_pair)
          NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                   dispersion_env%d2y_dx2, dispersion_env%dftd_section)
          NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
@@ -1244,6 +1249,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
             dispersion_env%rc_disp = se_control%rcdisp
             dispersion_env%exp_pre = 0._dp
             dispersion_env%scaling = 0._dp
+            dispersion_env%nd3_exclude_pair = 0
             dispersion_env%parameter_file_name = se_control%dispersion_parameter_file
             CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
          ELSE

tests/QS/regtest-dft-vdw-corr-2/TEST_FILES

@@ -10,6 +10,7 @@ dftd3_t11.inp                                          7      2e-12
 dftd3_t12.inp                                         33    1.0E-14              -0.00117296325816
 dftd3_t13.inp                                         33    1.0E-14              -0.00045037056834
 dftd3_t14.inp                                         33    1.0E-14              -0.00022349123009
+dftd3_t15.inp                                         33    1.0E-14              -0.00022349123009
 #
 argon-vdW-DF1.inp                                      1      3e-13             -85.04054534692069
 argon-vdW-DF2.inp                                      1      2e-13             -85.20254545254821

tests/QS/regtest-dft-vdw-corr-2/dftd3_t15.inp

@@ -0,0 +1,71 @@
+&GLOBAL
+  PROJECT dftd3_t12
+  RUN_TYPE ENERGY_FORCE
+  PRINT_LEVEL LOW
+&END GLOBAL
+&MOTION
+  &MD
+    ENSEMBLE NVE
+    STEPS 5
+    TIMESTEP 1.0
+    TEMPERATURE 300.0
+  &END
+&END
+&FORCE_EVAL
+  METHOD QS
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_MOLOPT
+    POTENTIAL_FILE_NAME POTENTIAL
+    &MGRID
+      CUTOFF 100
+    &END MGRID
+    &QS
+      METHOD GPW
+    &END QS
+    &SCF
+      SCF_GUESS ATOMIC
+      MAX_SCF 1
+      EPS_SCF 1.0e-5
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE 
+      &END XC_FUNCTIONAL
+      &vdW_POTENTIAL
+         DISPERSION_FUNCTIONAL PAIR_POTENTIAL
+         &PAIR_POTENTIAL
+            TYPE DFTD3
+            CALCULATE_C9_TERM .TRUE.
+            REFERENCE_C9_TERM .TRUE.
+            LONG_RANGE_CORRECTION .TRUE.
+            PARAMETER_FILE_NAME dftd3.dat
+            VERBOSE_OUTPUT
+            REFERENCE_FUNCTIONAL PBE
+            R_CUTOFF 8.
+            EPS_CN 0.01
+            D3_EXCLUDE_KIND_PAIR 2 2
+            D3_EXCLUDE_KIND_PAIR 1 2
+         &END PAIR_POTENTIAL
+      &END vdW_POTENTIAL
+    &END XC
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 8.30 8.30 8.30
+    &END
+    &COORD
+ SCALED
+ Ar    0.0 0.0 0.0
+ Ar    0.5 0.5 0.0
+ Ne    0.5 0.0 0.5
+ Ne    0.0 0.5 0.5
+    &END COORD
+    &KIND Ar
+      BASIS_SET SZV-MOLOPT-SR-GTH
+      POTENTIAL GTH-PBE-q8
+    &END
+    &KIND Ne
+      BASIS_SET SZV-MOLOPT-SR-GTH
+      POTENTIAL GTH-PBE-q8
+    &END
+  &END SUBSYS
+&END