TDDFPT transition dipole calculation: (#1975)

- new defaults for xTB
- correct velocity form for xTB

src/cp_control_utils.F

@@ -43,7 +43,7 @@ MODULE cp_control_utils
         real_time_propagation, sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, &
         sccs_derivative_cd7, sccs_derivative_fft, sccs_fattebert_gygi, sic_ad, sic_eo, &
         sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, slater, &
-        tddfpt_excitations, tddfpt_kernel_stda, use_mom_ref_user, xas_dip_len
+        tddfpt_dipole_length, tddfpt_excitations, tddfpt_kernel_stda, use_mom_ref_user
    USE input_cp2k_check,                ONLY: xc_functionals_expand
    USE input_cp2k_dft,                  ONLY: create_dft_section
    USE input_enumeration_types,         ONLY: enum_i2c,&
@@ -1328,14 +1328,19 @@ SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
 
       ! DIPOLE_MOMENTS subsection
       dipole_section => section_vals_get_subs_vals(t_section, "DIPOLE_MOMENTS")
-      CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", i_val=t_control%dipole_form)
+      CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", explicit=explicit)
+      IF (explicit) THEN
+         CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", i_val=t_control%dipole_form)
+      ELSE
+         t_control%dipole_form = 0
+      END IF
       CALL section_vals_val_get(dipole_section, "REFERENCE", i_val=t_control%dipole_reference)
       CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", explicit=explicit)
       IF (explicit) THEN
          CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", r_vals=t_control%dipole_ref_point)
       ELSE
          NULLIFY (t_control%dipole_ref_point)
-         IF (t_control%dipole_form == xas_dip_len .AND. t_control%dipole_reference == use_mom_ref_user) THEN
+         IF (t_control%dipole_form == tddfpt_dipole_length .AND. t_control%dipole_reference == use_mom_ref_user) THEN
             CPABORT("User-defined reference point should be given explicitly")
          END IF
       END IF

src/input_cp2k_properties_dft.F

@@ -1427,7 +1427,7 @@ SUBROUTINE create_tddfpt2_section(section)
 
       CALL keyword_create(keyword, __LOCATION__, name="DIPOLE_FORM", &
                           description="Form of dipole transition integrals.", &
-                          default_i_val=tddfpt_dipole_velocity, &
+                          type_of_var=integer_t, &
                           enum_c_vals=s2a("BERRY", "LENGTH", "VELOCITY"), &
                           enum_desc=s2a("Based on Berry phase formula (valid for fully periodic molecular systems only)", &
                                         "Length form ⟨ i | r | j ⟩ (valid for non-periodic molecular systems only)", &

src/qs_tddfpt2_properties.F

@@ -98,14 +98,17 @@ MODULE qs_tddfpt2_properties
    USE qs_environment_types,            ONLY: get_qs_env,&
                                               qs_environment_type
    USE qs_kind_types,                   ONLY: qs_kind_type
+   USE qs_ks_types,                     ONLY: qs_ks_env_type
    USE qs_mo_types,                     ONLY: allocate_mo_set,&
                                               deallocate_mo_set,&
                                               get_mo_set,&
                                               init_mo_set,&
                                               mo_set_p_type,&
                                               set_mo_set
    USE qs_moments,                      ONLY: build_berry_moment_matrix
+   USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
    USE qs_operators_ao,                 ONLY: rRc_xyz_ao
+   USE qs_overlap,                      ONLY: build_overlap_matrix
    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
                                               qs_subsys_type
    USE qs_tddfpt2_types,                ONLY: tddfpt_ground_state_mos
@@ -179,9 +182,13 @@ SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_
       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
       TYPE(cp_fm_type), POINTER                          :: ediff_inv, rRc_mos_occ, wfm_ao_ao, &
                                                             wfm_mo_virt_mo_occ
-      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: berry_cossin_xyz, matrix_s, rRc_xyz
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: berry_cossin_xyz, matrix_s, rRc_xyz, scrm
+      TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
+         POINTER                                         :: sab_orb
       TYPE(pw_env_type), POINTER                         :: pw_env
       TYPE(pw_poisson_type), POINTER                     :: poisson_env
+      TYPE(qs_ks_env_type), POINTER                      :: ks_env
 
       CALL timeset(routineN, handle)
 
@@ -222,6 +229,20 @@ SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_
       CALL pw_env_get(pw_env, poisson_env=poisson_env)
       ndim_periodic = COUNT(poisson_env%parameters%periodic == 1)
 
+      ! select default for dipole form
+      IF (tddfpt_control%dipole_form == 0) THEN
+         CALL get_qs_env(qs_env, dft_control=dft_control)
+         IF (dft_control%qs_control%xtb) THEN
+            IF (ndim_periodic == 0) THEN
+               tddfpt_control%dipole_form = tddfpt_dipole_length
+            ELSE
+               tddfpt_control%dipole_form = tddfpt_dipole_berry
+            END IF
+         ELSE
+            tddfpt_control%dipole_form = tddfpt_dipole_velocity
+         END IF
+      END IF
+
       SELECT CASE (tddfpt_control%dipole_form)
       CASE (tddfpt_dipole_berry)
          IF (ndim_periodic /= 3) THEN
@@ -434,9 +455,12 @@ SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_
          CALL dbcsr_deallocate_matrix_set(rRc_xyz)
 
       CASE (tddfpt_dipole_velocity)
-         IF (SIZE(matrix_s) < nderivs + 1) THEN
-            CPABORT("Where are the overlap matrix derivatives?")
-         END IF
+         ! generate overlap derivatives
+         CALL get_qs_env(qs_env, ks_env=ks_env, sab_orb=sab_orb)
+         NULLIFY (scrm)
+         CALL build_overlap_matrix(ks_env, matrix_s=scrm, nderivative=1, &
+                                   basis_type_a="ORB", basis_type_b="ORB", &
+                                   sab_nl=sab_orb)
 
          DO ispin = 1, nspins
             NULLIFY (fm_struct, ediff_inv, wfm_mo_virt_mo_occ)
@@ -466,7 +490,7 @@ SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_
 !$OMP       END PARALLEL DO
 
             DO ideriv = 1, nderivs
-               CALL cp_dbcsr_sm_fm_multiply(matrix_s(ideriv + 1)%matrix, &
+               CALL cp_dbcsr_sm_fm_multiply(scrm(ideriv + 1)%matrix, &
                                             gs_mos(ispin)%mos_occ, &
                                             dipole_op_mos_occ(ideriv, ispin)%matrix, &
                                             ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
@@ -497,6 +521,7 @@ SUBROUTINE tddfpt_dipole_operator(dipole_op_mos_occ, tddfpt_control, gs_mos, qs_
             CALL cp_fm_release(wfm_mo_virt_mo_occ)
             CALL cp_fm_release(ediff_inv)
          END DO
+         CALL dbcsr_deallocate_matrix_set(scrm)
 
       CASE DEFAULT
          CPABORT("Unimplemented form of the dipole operator")