Fix rtp partial occ density matrix (#4)

RTP partial occupation density matrix
cp2k · Jan 16, 2023 · d252d28 · d252d28
1 parent b578666
commit d252d28
Show file tree

Hide file tree

Showing 2 changed files with 58 additions and 15 deletions.
diff --git a/src/emd/rt_propagation_utils.F b/src/emd/rt_propagation_utils.F
@@ -408,25 +408,47 @@ SUBROUTINE calc_update_rho(qs_env)
       CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
       DO i = 1, SIZE(mos_new)/2
          re = 2*i - 1; im = 2*i
-         alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
          CALL dbcsr_set(rho_ao(i)%matrix, zero)
          CALL cp_fm_get_info(mos_new(re), ncol_global=ncol)
          CALL cp_fm_create(mos_occ, &
                            matrix_struct=mos(i)%mo_coeff%matrix_struct, &
                            name="mos_occ")
          CALL cp_fm_to_fm(mos_new(re), mos_occ)
-         CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
-         CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
-                                    matrix_v=mos_occ, &
-                                    ncol=ncol, &
-                                    alpha=alpha)
+         IF (mos(i)%uniform_occupation) THEN
+            alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
+            CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
+                                       matrix_v=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha)
+         ELSE
+            alpha = 1.0_dp
+            CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
+                                       matrix_v=mos_new(re), &
+                                       matrix_g=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha)
+         END IF
+
          ! It is actually complex conjugate but i*i=-1 therefore it must be added
          CALL cp_fm_to_fm(mos_new(im), mos_occ)
-         CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
-         CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
-                                    matrix_v=mos_occ, &
-                                    ncol=ncol, &
-                                    alpha=alpha)
+         IF (mos(i)%uniform_occupation) THEN
+            alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
+            CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
+                                       matrix_v=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha)
+         ELSE
+            alpha = 1.0_dp
+            CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
+                                       matrix_v=mos_new(im), &
+                                       matrix_g=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha)
+         END IF
          CALL cp_fm_release(mos_occ)
       END DO
       CALL qs_rho_update_rho(rho, qs_env)

diff --git a/src/emd/rt_propagator_init.F b/src/emd/rt_propagator_init.F
@@ -418,20 +418,41 @@ SUBROUTINE rt_initialize_rho_from_mos(rtp, mos)
 
       INTEGER                                            :: ispin, ncol, re
       REAL(KIND=dp)                                      :: alpha
+      TYPE(cp_fm_type)                                   :: mos_occ
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_new, rho_old
 
       CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
 
       DO ispin = 1, SIZE(mos)
          re = 2*ispin - 1
-         alpha = 3.0_dp - REAL(SIZE(mos), dp)
          CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
          CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
-         CALL cp_fm_column_scale(mos(ispin)%mo_coeff, mos(ispin)%occupation_numbers/alpha)
-         CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
-                                    matrix_v=mos(ispin)%mo_coeff, ncol=ncol, alpha=alpha, keep_sparsity=.FALSE.)
+
+         CALL cp_fm_create(mos_occ, &
+                           matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
+                           name="mos_occ")
+         CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
+         IF (mos(ispin)%uniform_occupation) THEN
+            alpha = 3.0_dp - REAL(SIZE(mos), dp)
+            CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
+                                       matrix_v=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha, keep_sparsity=.FALSE.)
+         ELSE
+            alpha = 1.0_dp
+            CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
+            CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
+                                       matrix_v=mos(ispin)%mo_coeff, &
+                                       matrix_g=mos_occ, &
+                                       ncol=ncol, &
+                                       alpha=alpha, keep_sparsity=.FALSE.)
+         END IF
+
          CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
          CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
+         CALL cp_fm_release(mos_occ)
+
       END DO
 
    END SUBROUTINE rt_initialize_rho_from_mos