From ada2ee8b55713fac823189a53727452fbcfe12c0 Mon Sep 17 00:00:00 2001 From: marcella <40494362+marci73@users.noreply.github.com> Date: Mon, 26 Jan 2026 18:56:45 +0100 Subject: [PATCH] BUG FIX: EMD + LSD + ADMM (#4739) --- src/rt_propagation_forces.F | 41 +++++++++++++++++++------------------ 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/src/rt_propagation_forces.F b/src/rt_propagation_forces.F index f8a7100a2b..0990ba870a 100644 --- a/src/rt_propagation_forces.F +++ b/src/rt_propagation_forces.F @@ -301,35 +301,35 @@ CONTAINS END DO ! First calculate H=KS_aux*C~, real part ends on work_aux_aux2, imaginary part ends at work_aux_aux3 - CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(re), tmp_aux_mo(re), nmo, 4.0_dp, 0.0_dp) - CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(im), tmp_aux_mo(im), nmo, 4.0_dp, 0.0_dp) - CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(im), tmp_aux_mo(re), nmo, -4.0_dp, 1.0_dp) - CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(re), tmp_aux_mo(im), nmo, 4.0_dp, 1.0_dp) + CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(re), tmp_aux_mo(1), nmo, 4.0_dp, 0.0_dp) + CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(im), tmp_aux_mo(2), nmo, 4.0_dp, 0.0_dp) + CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(im), tmp_aux_mo(1), nmo, -4.0_dp, 1.0_dp) + CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(re), tmp_aux_mo(2), nmo, 4.0_dp, 1.0_dp) ! Next step compute S-1*H - CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(re), 0.0_dp, tmp_aux_mo1(re)) - CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(im), 0.0_dp, tmp_aux_mo1(im)) + CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(1), 0.0_dp, tmp_aux_mo1(1)) + CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(2), 0.0_dp, tmp_aux_mo1(2)) ! Here we go on with Ws=S-1*H * C^H (take care of sign of the imaginary part!!!) - CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(re), mos(re), 0.0_dp, & - tmp_aux_nao(re)) + CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(1), mos(re), 0.0_dp, & + tmp_aux_nao(1)) CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(im), 1.0_dp, & - tmp_aux_nao(re)) - CALL parallel_gemm("N", "T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(re), mos(im), 0.0_dp, & - tmp_aux_nao(im)) - CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(re), 1.0_dp, & - tmp_aux_nao(im)) + tmp_aux_nao(1)) + CALL parallel_gemm("N", "T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(1), mos(im), 0.0_dp, & + tmp_aux_nao(2)) + CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(2), mos(re), 1.0_dp, & + tmp_aux_nao(2)) ! Let's do the final bit Wq=S-1*H * C^H * A^T - CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(re), admm_env%A, 0.0_dp, tmp_aux_aux(re)) - CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(im), admm_env%A, 0.0_dp, tmp_aux_aux(im)) + CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(1), admm_env%A, 0.0_dp, tmp_aux_aux(1)) + CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(2), admm_env%A, 0.0_dp, tmp_aux_aux(2)) ! *** copy to sparse matrix - CALL copy_fm_to_dbcsr(tmp_aux_nao(re), matrix_w_q, keep_sparsity=.TRUE.) + CALL copy_fm_to_dbcsr(tmp_aux_nao(1), matrix_w_q, keep_sparsity=.TRUE.) ! *** copy to sparse matrix - CALL copy_fm_to_dbcsr(tmp_aux_aux(re), matrix_w_s, keep_sparsity=.TRUE.) + CALL copy_fm_to_dbcsr(tmp_aux_aux(1), matrix_w_s, keep_sparsity=.TRUE.) DO jspin = 1, 2 CALL cp_fm_release(tmp_aux_mo(jspin)) @@ -353,11 +353,12 @@ CONTAINS CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set) DEALLOCATE (admm_force) - ! *** Deallocated weighted density matrices - CALL dbcsr_deallocate_matrix(matrix_w_s) - CALL dbcsr_deallocate_matrix(matrix_w_q) END DO + ! *** Deallocated weighted density matrices + CALL dbcsr_deallocate_matrix(matrix_w_s) + CALL dbcsr_deallocate_matrix(matrix_w_q) + DO jspin = 1, 2 CALL cp_fm_release(tmp_aux_aux(jspin)) CALL cp_fm_release(tmp_aux_nao(jspin))