comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

BUG FIX: EMD + LSD + ADMM (#4739)

Browse source
This commit is contained in:
marcella 2026-01-26 18:56:45 +01:00 • committed by GitHub
parent b049e7ca4f
commit ada2ee8b55
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
1 changed files with 21 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

41
src/rt_propagation_forces.F
View file

@ -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))