Add SCF_MOMENT dipoles for k-point TDDFPT KERNEL NONE (#5199)

Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
2026-07-28 22:25:32 -04:00 · 2026-05-16 16:08:05 +02:00 · 2026-05-16 16:08:05 +02:00 · 479cecb1b0
commit 479cecb1b0
parent a91f8bf3b4
6 changed files with 154 additions and 37 deletions
--- a/src/input_constants.F
+++ b/src/input_constants.F
@ -712,7 +712,8 @@ MODULE input_constants
   INTEGER, PARAMETER, PUBLIC               :: tddfpt_dipole_berry = 1, &
                                               tddfpt_dipole_length = 2, &
                                               tddfpt_dipole_velocity = 3, &
-                                               tddfpt_dipole_velocity_old = 4
+                                               tddfpt_dipole_velocity_old = 4, &
+                                               tddfpt_dipole_scf_moment = 5

   ! XC Kernel derivative methods for forces
   INTEGER, PARAMETER, PUBLIC               :: xc_kernel_method_best = 100, &
--- a/src/input_cp2k_properties_dft.F
+++ b/src/input_cp2k_properties_dft.F
@ -40,10 +40,10 @@ MODULE input_cp2k_properties_dft
      oe_lb, oe_none, oe_saop, oe_shift, ot_precond_full_all, ot_precond_full_kinetic, &
      ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
      ot_precond_s_inverse, scan_x, scan_xy, scan_xyz, scan_xz, scan_y, scan_yz, scan_z, &
-      tddfpt_dipole_berry, tddfpt_dipole_length, tddfpt_dipole_velocity, tddfpt_dipole_velocity_old, &
-      tddfpt_kernel_full, &
-      tddfpt_kernel_none, tddfpt_kernel_stda, no_sf_tddfpt, tddfpt_sf_col, tddfpt_sf_noncol, &
-      use_mom_ref_coac, use_mom_ref_com, use_mom_ref_user, use_mom_ref_zero
+      tddfpt_dipole_berry, tddfpt_dipole_length, tddfpt_dipole_scf_moment, &
+      tddfpt_dipole_velocity, tddfpt_dipole_velocity_old, tddfpt_kernel_full, tddfpt_kernel_none, &
+      tddfpt_kernel_stda, no_sf_tddfpt, tddfpt_sf_col, tddfpt_sf_noncol, use_mom_ref_coac, &
+      use_mom_ref_com, use_mom_ref_user, use_mom_ref_zero
   USE input_cp2k_atprop, ONLY: create_atprop_section
   USE input_cp2k_dft, ONLY: create_interp_section, &
                             create_mgrid_section
@ -1848,13 +1848,16 @@ CONTAINS

      CALL keyword_create(keyword, __LOCATION__, name="DIPOLE_FORM", &
                          description="Form of dipole transition integrals.", &
-                          enum_c_vals=s2a("BERRY", "LENGTH", "VELOCITY", "VELOCITY_OLD"), &
+                          enum_c_vals=s2a("BERRY", "LENGTH", "VELOCITY", "VELOCITY_OLD", &
+                                          "SCF_MOMENT"), &
                          enum_desc=s2a("Based on Berry phase formula (valid for fully periodic molecular systems only)", &
                                        "Length form &lang; i | r | j &rang; (valid for non-periodic molecular systems only)", &
                                        "Velocity form &lang; i | d/dr | j &rang;", &
-                                        "Old velocity form &lang; i | d/dr | j &rang;"), &
+                                        "Old velocity form &lang; i | d/dr | j &rang;", &
+                                        "SCF molecular-orbital moment form for k-point TDDFPT"), &
                          enum_i_vals=[tddfpt_dipole_berry, tddfpt_dipole_length, &
-                                       tddfpt_dipole_velocity, tddfpt_dipole_velocity_old], &
+                                       tddfpt_dipole_velocity, tddfpt_dipole_velocity_old, &
+                                       tddfpt_dipole_scf_moment], &
                          default_i_val=tddfpt_dipole_velocity)
      CALL section_add_keyword(subsection, keyword)
      CALL keyword_release(keyword)
--- a/src/qs_tddfpt2_methods.F
+++ b/src/qs_tddfpt2_methods.F
@ -59,8 +59,8 @@ MODULE qs_tddfpt2_methods
   USE input_constants,                 ONLY: &
        do_admm_aux_exch_func_none, do_admm_basis_projection, do_admm_exch_scaling_none, &
        do_admm_purify_none, do_potential_truncated, no_sf_tddfpt, oe_none, &
-        tddfpt_dipole_velocity, tddfpt_kernel_full, tddfpt_kernel_none, tddfpt_kernel_stda, &
-        tddfpt_sf_col, tddfpt_sf_noncol
+        tddfpt_dipole_scf_moment, tddfpt_dipole_velocity, tddfpt_kernel_full, tddfpt_kernel_none, &
+        tddfpt_kernel_stda, tddfpt_sf_col, tddfpt_sf_noncol
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
@ -692,8 +692,9 @@ CONTAINS
         IF (tddfpt_control%oe_corr /= oe_none) &
            CPABORT("Orbital-energy-corrected TDDFPT is not implemented for k-points")
         IF (tddfpt_control%dipole_form /= 0 .AND. &
-             tddfpt_control%dipole_form /= tddfpt_dipole_velocity) &
-            CPABORT("K-point TDDFPT supports only velocity-form transition dipoles")
+             tddfpt_control%dipole_form /= tddfpt_dipole_velocity .AND. &
+             tddfpt_control%dipole_form /= tddfpt_dipole_scf_moment) &
+            CPABORT("K-point TDDFPT supports only velocity-form or SCF_MOMENT transition dipoles")
      END IF

      IF (tddfpt_control%nstates <= 0) THEN
@ -995,7 +996,7 @@ CONTAINS
      INTEGER, DIMENSION(2)                              :: kp_range
      INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
      INTEGER, DIMENSION(maxspins)                       :: homo_spin, nao_spin, nmo_spin, nvirt_spin
-      LOGICAL                                            :: my_kpgrp
+      LOGICAL                                            :: my_kpgrp, use_scf_moment_dipoles
      REAL(kind=dp)                                      :: checksum, dipole_im, dipole_re, fsum, &
                                                            gap, oscillator_factor, spin_factor
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues_kp, evals, &
@ -1085,16 +1086,28 @@ CONTAINS
      transition_dipole_re = 0.0_dp
      transition_dipole_im = 0.0_dp
      oscillator_strength = 0.0_dp
+      use_scf_moment_dipoles = (tddfpt_control%dipole_form == tddfpt_dipole_scf_moment)
+      IF (use_scf_moment_dipoles) THEN
+         CALL cp_warn(__LOCATION__, "SCF_MOMENT k-point dipoles use direct SCF MO matrix "// &
+                      "elements; compare folded energy blocks, not individual degenerate states.")
+      END IF

      CALL build_overlap_matrix(ks_env, matrixkp_s=overlap_deriv, nderivative=1, &
                                basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb, &
                                ext_kpoints=kpoints)

      ALLOCATE (rmatrix, cmatrix)
-      CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
-                        matrix_type=dbcsr_type_symmetric)
-      CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
-                        matrix_type=dbcsr_type_antisymmetric)
+      IF (use_scf_moment_dipoles) THEN
+         CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
+                           matrix_type=dbcsr_type_antisymmetric)
+         CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
+                           matrix_type=dbcsr_type_symmetric)
+      ELSE
+         CALL dbcsr_create(rmatrix, template=overlap_deriv(1, 1)%matrix, &
+                           matrix_type=dbcsr_type_symmetric)
+         CALL dbcsr_create(cmatrix, template=overlap_deriv(1, 1)%matrix, &
+                           matrix_type=dbcsr_type_antisymmetric)
+      END IF
      CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_kp)
      CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_kp)

@ -1158,29 +1171,55 @@ CONTAINS
                                   ispin=ideriv + 1, xkp=kpoints%xkp(:, ikp), &
                                   cell_to_index=cell_to_index, sab_nl=sab_kp)

-               CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)
+               IF (use_scf_moment_dipoles) THEN
+                  CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)

-               CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
+                  CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)

-               CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
+                  CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)

-               CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
-               CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
-                                  -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
+                  CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
+               ELSE
+                  CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 0.0_dp, moment_re)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_im_global, fm_tmp, 0.0_dp, moment_im)
+
+                  CALL cp_dbcsr_sm_fm_multiply(rmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
+
+                  CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_re_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_im)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_re)
+
+                  CALL cp_dbcsr_sm_fm_multiply(cmatrix, mo_coeff_im_global, fm_tmp, nmo_spin(ispin))
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_re_global, fm_tmp, 1.0_dp, moment_re)
+                  CALL parallel_gemm("T", "N", nmo_spin(ispin), nmo_spin(ispin), nao, &
+                                     -1.0_dp, mo_coeff_im_global, fm_tmp, 1.0_dp, moment_im)
+               END IF

               DO iocc = 1, homo_spin(ispin)
                  DO ivirt = homo_spin(ispin) + 1, nmo_spin(ispin)
--- a/src/qs_tddfpt2_properties.F
+++ b/src/qs_tddfpt2_properties.F
@ -60,6 +60,7 @@ MODULE qs_tddfpt2_properties
   USE input_constants,                 ONLY: no_sf_tddfpt,&
                                              tddfpt_dipole_berry,&
                                              tddfpt_dipole_length,&
+                                              tddfpt_dipole_scf_moment,&
                                              tddfpt_dipole_velocity,&
                                              tddfpt_dipole_velocity_old
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
@ -572,6 +573,8 @@ CONTAINS
            WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using velocity formulation"
         CASE (tddfpt_dipole_velocity_old)
            WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using old velocity formulation"
+         CASE (tddfpt_dipole_scf_moment)
+            WRITE (log_unit, '(1X,A,/)') "Transition dipoles calculated using SCF-MO moment formulation"
         CASE DEFAULT
            CPABORT("Unimplemented form of the dipole operator")
         END SELECT
--- a/tests/QS/regtest-tddfpt/H2_kp_tddfpt_none_scf_moment.inp
+++ b/tests/QS/regtest-tddfpt/H2_kp_tddfpt_none_scf_moment.inp
@ -0,0 +1,69 @@
+&GLOBAL
+  PRINT_LEVEL LOW
+  PROJECT H2_kp_tddfpt_none_scf_moment
+  RUN_TYPE ENERGY
+&END GLOBAL
+
+&FORCE_EVAL
+  METHOD Quickstep
+  &DFT
+    BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH
+    POTENTIAL_FILE_NAME POTENTIAL_UZH
+    &KPOINTS
+      FULL_GRID ON
+      PARALLEL_GROUP_SIZE -1
+      SCHEME MONKHORST-PACK 2 1 1
+      WAVEFUNCTIONS COMPLEX
+    &END KPOINTS
+    &MGRID
+      CUTOFF 100
+      REL_CUTOFF 30
+    &END MGRID
+    &QS
+      EPS_DEFAULT 1.0E-10
+      METHOD GPW
+    &END QS
+    &SCF
+      ADDED_MOS 4
+      CHOLESKY OFF
+      EPS_EIGVAL 1.0E-8
+      EPS_SCF 1.0E-8
+      MAX_SCF 50
+      SCF_GUESS ATOMIC
+      &MIXING
+        ALPHA 0.35
+        METHOD BROYDEN_MIXING
+      &END MIXING
+      &PRINT
+        &RESTART OFF
+        &END RESTART
+      &END PRINT
+    &END SCF
+    &XC
+      &XC_FUNCTIONAL PBE
+      &END XC_FUNCTIONAL
+    &END XC
+  &END DFT
+  &PROPERTIES
+    &TDDFPT
+      KERNEL NONE
+      NSTATES 2
+      &DIPOLE_MOMENTS
+        DIPOLE_FORM SCF_MOMENT
+      &END DIPOLE_MOMENTS
+    &END TDDFPT
+  &END PROPERTIES
+  &SUBSYS
+    &CELL
+      ABC 6.0 6.0 6.0
+    &END CELL
+    &COORD
+      H 2.65 3.00 3.00
+      H 3.35 3.00 3.00
+    &END COORD
+    &KIND H
+      BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q1
+      POTENTIAL GTH-GGA-q1
+    &END KIND
+  &END SUBSYS
+&END FORCE_EVAL
--- a/tests/QS/regtest-tddfpt/TEST_FILES.toml
+++ b/tests/QS/regtest-tddfpt/TEST_FILES.toml
@ -39,5 +39,7 @@
                                           {matcher="TDDFPT_Check_Osc_Strength", tol=4.0E-06, ref=0.36697616E+00}]
 "H2_kp_tddfpt_none_sym.inp"             = [{matcher="TDDFPT_Check_Energy", tol=4.0E-06, ref=0.465466E+00},
                                           {matcher="TDDFPT_Check_Osc_Strength", tol=4.0E-06, ref=0.51898266E+00}]
+"H2_kp_tddfpt_none_scf_moment.inp"      = [{matcher="TDDFPT_Check_Energy", tol=4.0E-06, ref=0.658268E+00},
+                                           {matcher="TDDFPT_Check_Osc_Strength", tol=4.0E-06, ref=0.42166952E+00}]
 #
 #EOF