diff --git a/src/input_constants.F b/src/input_constants.F index b349dc3a68..7bf4b498d2 100644 --- 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, & diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index 37bf74cfae..dd7e78a90f 100644 --- 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 ⟨ i | r | j ⟩ (valid for non-periodic molecular systems only)", & "Velocity form ⟨ i | d/dr | j ⟩", & - "Old velocity form ⟨ i | d/dr | j ⟩"), & + "Old velocity form ⟨ i | d/dr | j ⟩", & + "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) diff --git a/src/qs_tddfpt2_methods.F b/src/qs_tddfpt2_methods.F index 9fba9c3675..71684993fc 100644 --- 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) diff --git a/src/qs_tddfpt2_properties.F b/src/qs_tddfpt2_properties.F index 8f45fe2da7..c5432337c3 100644 --- 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 diff --git a/tests/QS/regtest-tddfpt/H2_kp_tddfpt_none_scf_moment.inp b/tests/QS/regtest-tddfpt/H2_kp_tddfpt_none_scf_moment.inp new file mode 100644 index 0000000000..18da028132 --- /dev/null +++ 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 diff --git a/tests/QS/regtest-tddfpt/TEST_FILES.toml b/tests/QS/regtest-tddfpt/TEST_FILES.toml index 9ca1f80012..06d479acb4 100644 --- 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