From 8817145e271dba2fc8cd57e0e6bba2cd65c5ad93 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 13 Feb 2023 20:12:33 +0100 Subject: [PATCH 1/2] minor modifs --- external/qp2-dependencies | 2 +- src/dft_utils_func/rho_ab_to_rho_tot.irp.f | 35 +++++++++++----------- src/dft_utils_in_r/dm_in_r.irp.f | 3 ++ 3 files changed, 21 insertions(+), 19 deletions(-) diff --git a/external/qp2-dependencies b/external/qp2-dependencies index b8cd5815..f40bde09 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit b8cd5815bce14c9b880e3c5ea3d5fc2652f5955c +Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 diff --git a/src/dft_utils_func/rho_ab_to_rho_tot.irp.f b/src/dft_utils_func/rho_ab_to_rho_tot.irp.f index 919543fe..d4d5a536 100644 --- a/src/dft_utils_func/rho_ab_to_rho_tot.irp.f +++ b/src/dft_utils_func/rho_ab_to_rho_tot.irp.f @@ -66,10 +66,27 @@ subroutine v_rho_oc_to_v_rho_ab(v_rho_o,v_rho_c,v_rho_a,v_rho_b) END_DOC double precision, intent(in) :: v_rho_o,v_rho_c double precision, intent(out) :: v_rho_a,v_rho_b +! print*,'in v_rho_oc_to_v_rho_ab' +! print*, v_rho_c , v_rho_o v_rho_a = v_rho_c + v_rho_o v_rho_b = v_rho_c - v_rho_o end +subroutine v_grad_rho_ab_to_v_grad_rho_oc(v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b,v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c) + implicit none + double precision, intent(in) :: v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b + double precision, intent(out) :: v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c + BEGIN_DOC +! convert (v_grad_rho_a_2, v_grad_rho_b_2, v_grad_rho_a.grad_rho_b) +! +! to (v_grad_rho_c_2, v_grad_rho_o_2, v_grad_rho_o.grad_rho_c) +! +! rho_c = total density, rho_o spin density + END_DOC + v_grad_rho_c_2 = 0.25d0 * (v_grad_rho_a_2 + v_grad_rho_b_2 + v_grad_rho_a_b) + v_grad_rho_o_2 = 0.25d0 * (v_grad_rho_a_2 + v_grad_rho_b_2 - v_grad_rho_a_b) + v_grad_rho_o_c = 0.25d0 * (2d0 * v_grad_rho_a_2 - 2d0 * v_grad_rho_b_2 ) +end subroutine v_grad_rho_oc_to_v_grad_rho_ab(v_grad_rho_o_2,v_grad_rho_c_2,v_grad_rho_o_c,v_grad_rho_a_2,v_grad_rho_b_2,v_grad_rho_a_b) @@ -88,21 +105,3 @@ subroutine v_grad_rho_oc_to_v_grad_rho_ab(v_grad_rho_o_2,v_grad_rho_c_2,v_grad_r v_grad_rho_a_b = -2d0 * v_grad_rho_o_2 + 2d0 * v_grad_rho_c_2 end - - - - - - - - - - - - - - - - - - diff --git a/src/dft_utils_in_r/dm_in_r.irp.f b/src/dft_utils_in_r/dm_in_r.irp.f index 53e15b06..feb174fd 100644 --- a/src/dft_utils_in_r/dm_in_r.irp.f +++ b/src/dft_utils_in_r/dm_in_r.irp.f @@ -45,6 +45,8 @@ call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array) ! alpha/beta density + dm_a(istate) = max(dm_a(istate),1.d-12) + dm_b(istate) = max(dm_b(istate),1.d-12) one_e_dm_and_grad_alpha_in_r(4,i,istate) = dm_a(istate) one_e_dm_and_grad_beta_in_r(4,i,istate) = dm_b(istate) @@ -80,6 +82,7 @@ enddo enddo !$OMP END PARALLEL DO + print*,'density and gradients provided' END_PROVIDER From 69f014bc9cdde7c01c9481ecff86935ef522907e Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 14 Feb 2023 23:45:15 +0100 Subject: [PATCH 2/2] added nuclear repulsion in the diagonal TC matrix element bi_ortho --- src/tc_bi_ortho/slater_tc_opt_diag.irp.f | 2 +- src/two_body_rdm/example.irp.f | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f index 68f647dd..7524e11a 100644 --- a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -16,7 +16,7 @@ else ref_tc_energy_3e = 0.d0 endif - ref_tc_energy_tot = ref_tc_energy_1e + ref_tc_energy_2e + ref_tc_energy_3e + ref_tc_energy_tot = ref_tc_energy_1e + ref_tc_energy_2e + ref_tc_energy_3e + nuclear_repulsion END_PROVIDER subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, htot) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 6cf0209e..de3d97b9 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -17,6 +17,7 @@ subroutine routine_active_only double precision :: wee_ab_st_av, rdm_ab_st_av double precision :: wee_tot_st_av, rdm_tot_st_av,spin_trace double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 + double precision :: accu_aa, accu_bb, accu_ab, accu_tot wee_ab = 0.d0 wee_bb = 0.d0 @@ -64,14 +65,23 @@ subroutine routine_active_only do istate = 1, N_states !! PURE ACTIVE PART !! + accu_aa = 0.d0 + accu_bb = 0.d0 + accu_ab = 0.d0 + accu_tot = 0.d0 do i = 1, n_act_orb iorb = list_act(i) do j = 1, n_act_orb jorb = list_act(j) + accu_bb += act_2_rdm_bb_mo(j,i,j,i,1) + accu_aa += act_2_rdm_aa_mo(j,i,j,i,1) + accu_ab += act_2_rdm_ab_mo(j,i,j,i,1) + accu_tot += act_2_rdm_spin_trace_mo(j,i,j,i,1) do k = 1, n_act_orb korb = list_act(k) do l = 1, n_act_orb lorb = list_act(l) + ! 1 2 1 2 2 1 2 1 if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)).gt.1.d-10)then print*,'Error in act_2_rdm_spin_trace_mo' print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l) - act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10" @@ -79,6 +89,7 @@ subroutine routine_active_only print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(j,i,l,k,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)) endif + ! 1 2 1 2 1 2 1 2 if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)).gt.1.d-10)then print*,'Error in act_2_rdm_spin_trace_mo' print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate),istate).gt.1.d-10" @@ -131,6 +142,15 @@ subroutine routine_active_only print*,'wee_tot = ',wee_tot(istate) print*,'Full energy ' print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) + print*,'--------------------------' + print*,'accu_aa = ',accu_aa + print*,'N_a (N_a-1)/2 = ', elec_alpha_num*(elec_alpha_num-1)*0.5 + print*,'accu_bb = ',accu_bb + print*,'N_b (N_b-1)/2 = ', elec_beta_num*(elec_beta_num-1)*0.5 + print*,'accu_ab = ',accu_ab + print*,'N_a N_b = ', elec_beta_num*elec_alpha_num + print*,'accu_tot = ',accu_tot + print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num * 0.5 enddo wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0