From 4e0c71df10ccf6653fbc802692e6bba4f663965f Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 6 Feb 2017 21:28:01 +0100 Subject: [PATCH] density based mrpt2 --- plugins/MRPT/print_1h2p.irp.f | 7 +- plugins/MRPT_Utils/density_matrix_based.irp.f | 193 +++++++++++++++++ plugins/MRPT_Utils/energies_cas.irp.f | 28 ++- plugins/MRPT_Utils/excitations_cas.irp.f | 9 +- plugins/MRPT_Utils/mrpt_dress.irp.f | 17 +- plugins/MRPT_Utils/mrpt_utils.irp.f | 64 +++--- plugins/MRPT_Utils/new_way.irp.f | 5 +- plugins/MRPT_Utils/psi_active_prov.irp.f | 1 + src/Determinants/density_matrix.irp.f | 51 +++++ src/Determinants/two_body_dm_map.irp.f | 199 ++++-------------- 10 files changed, 367 insertions(+), 207 deletions(-) create mode 100644 plugins/MRPT_Utils/density_matrix_based.irp.f diff --git a/plugins/MRPT/print_1h2p.irp.f b/plugins/MRPT/print_1h2p.irp.f index 85ddcda8..b9f6575b 100644 --- a/plugins/MRPT/print_1h2p.irp.f +++ b/plugins/MRPT/print_1h2p.irp.f @@ -2,7 +2,7 @@ program print_1h2p implicit none read_wf = .True. touch read_wf - call routine_2 + call routine end subroutine routine_2 @@ -35,7 +35,8 @@ subroutine routine integer :: i,j,istate accu = 0.d0 matrix_1h2p = 0.d0 - call H_apply_mrpt_2p(matrix_1h2p,N_det_ref) +!call H_apply_mrpt_1h2p(matrix_1h2p,N_det_ref) + call give_1h2p_contrib(matrix_1h2p) do istate = 1, N_states do i = 1, N_det do j = 1, N_det @@ -44,6 +45,8 @@ subroutine routine enddo print*,accu(istate) enddo + call contrib_1h2p_dm_based(accu) + print*,accu(:) deallocate (matrix_1h2p) end diff --git a/plugins/MRPT_Utils/density_matrix_based.irp.f b/plugins/MRPT_Utils/density_matrix_based.irp.f new file mode 100644 index 00000000..ac135807 --- /dev/null +++ b/plugins/MRPT_Utils/density_matrix_based.irp.f @@ -0,0 +1,193 @@ +subroutine contrib_1h2p_dm_based(accu) + implicit none + integer :: i_i,i_r,i_v,i_a,i_b + integer :: i,r,v,a,b + integer :: ispin,jspin + integer :: istate + double precision, intent(out) :: accu(N_states) + double precision :: active_int(n_act_orb,2) + double precision :: delta_e(n_act_orb,2,N_states) + double precision :: get_mo_bielec_integral + accu = 0.d0 +!do i_i = 1, 1 + do i_i = 1, n_inact_orb + i = list_inact(i_i) +! do i_r = 1, 1 + do i_r = 1, n_virt_orb + r = list_virt(i_r) +! do i_v = 1, 1 + do i_v = 1, n_virt_orb + v = list_virt(i_v) + do i_a = 1, n_act_orb + a = list_act(i_a) + active_int(i_a,1) = get_mo_bielec_integral(i,a,r,v,mo_integrals_map) ! direct + active_int(i_a,2) = get_mo_bielec_integral(i,a,v,r,mo_integrals_map) ! exchange + do istate = 1, N_states + do jspin=1, 2 + delta_e(i_a,jspin,istate) = one_anhil(i_a,jspin,istate) & + - fock_virt_total_spin_trace(r,istate) & + - fock_virt_total_spin_trace(v,istate) & + + fock_core_inactive_total_spin_trace(i,istate) + delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate) + enddo + enddo + enddo + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb +! do i_b = i_a, i_a + b = list_act(i_b) + do ispin = 1, 2 ! spin of (i --> r) + do jspin = 1, 2 ! spin of (a --> v) + if(ispin == jspin .and. r.le.v)cycle ! condition not to double count + do istate = 1, N_states + if(ispin == jspin)then + accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_mo_spin_index(a,b,istate,ispin) & + * (active_int(i_b,1) - active_int(i_b,2)) & + * delta_e(i_a,jspin,istate) + else + accu(istate) += active_int(i_a,1) * one_body_dm_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & + * active_int(i_b,1) + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + +end + +subroutine contrib_2h1p_dm_based(accu) + implicit none + integer :: i_i,i_j,i_v,i_a,i_b + integer :: i,j,v,a,b + integer :: ispin,jspin + integer :: istate + double precision, intent(out) :: accu(N_states) + double precision :: active_int(n_act_orb,2) + double precision :: delta_e(n_act_orb,2,N_states) + double precision :: get_mo_bielec_integral + accu = 0.d0 + do i_i = 1, n_inact_orb + i = list_inact(i_i) + do i_j = 1, n_inact_orb + j = list_inact(i_j) + do i_v = 1, n_virt_orb + v = list_virt(i_v) + do i_a = 1, n_act_orb + a = list_act(i_a) + active_int(i_a,1) = get_mo_bielec_integral(i,j,v,a,mo_integrals_map) ! direct + active_int(i_a,2) = get_mo_bielec_integral(i,j,a,v,mo_integrals_map) ! exchange + do istate = 1, N_states + do jspin=1, 2 +! delta_e(i_a,jspin,istate) = +! + delta_e(i_a,jspin,istate) = one_creat(i_a,jspin,istate) - fock_virt_total_spin_trace(v,istate) & + + fock_core_inactive_total_spin_trace(i,istate) & + + fock_core_inactive_total_spin_trace(j,istate) + delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate) + enddo + enddo + enddo + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb +! do i_b = i_a, i_a + b = list_act(i_b) + do ispin = 1, 2 ! spin of (i --> v) + do jspin = 1, 2 ! spin of (j --> a) + if(ispin == jspin .and. i.le.j)cycle ! condition not to double count + do istate = 1, N_states + if(ispin == jspin)then + accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) & + * (active_int(i_b,1) - active_int(i_b,2)) & + * delta_e(i_a,jspin,istate) + else + accu(istate) += active_int(i_a,1) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & + * active_int(i_b,1) + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + +end + + +subroutine contrib_2p_dm_based(accu) + implicit none + integer :: i_r,i_v,i_a,i_b,i_c,i_d + integer :: r,v,a,b,c,d + integer :: ispin,jspin + integer :: istate + double precision, intent(out) :: accu(N_states) + double precision :: active_int(n_act_orb,n_act_orb,2) + double precision :: delta_e(n_act_orb,n_act_orb,2,2,N_states) + double precision :: get_mo_bielec_integral + accu = 0.d0 + do i_r = 1, n_virt_orb + r = list_virt(i_r) + do i_v = 1, n_virt_orb + v = list_virt(i_v) + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb + b = list_act(i_b) + active_int(i_a,i_b,1) = get_mo_bielec_integral(a,b,r,v,mo_integrals_map) ! direct + active_int(i_a,i_b,2) = get_mo_bielec_integral(a,b,v,r,mo_integrals_map) ! direct + do istate = 1, N_states + do jspin=1, 2 ! spin of i_a + do ispin = 1, 2 ! spin of i_b + delta_e(i_a,i_b,jspin,ispin,istate) = two_anhil(i_a,i_b,jspin,ispin,istate) & + - fock_virt_total_spin_trace(r,istate) & + - fock_virt_total_spin_trace(v,istate) + delta_e(i_a,i_b,jspin,ispin,istate) = 1.d0/delta_e(i_a,i_b,jspin,ispin,istate) + enddo + enddo + enddo + enddo + enddo + ! diagonal terms + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb + b = list_act(i_b) + do ispin = 1, 2 ! spin of (a --> r) + do jspin = 1, 2 ! spin of (b --> v) + if(ispin == jspin .and. r.le.v)cycle ! condition not to double count + if(ispin == jspin .and. a.le.b)cycle ! condition not to double count + do istate = 1, N_states + if(ispin == jspin)then + double precision :: contrib_spin + if(ispin == 1)then + contrib_spin = two_body_dm_aa_diag_act(i_a,i_b) + else + contrib_spin = two_body_dm_bb_diag_act(i_a,i_b) + endif + accu(istate) += (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) * contrib_spin & + * (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) & + * delta_e(i_a,i_b,ispin,jspin,istate) + else + accu(istate) += 0.5d0 * active_int(i_a,i_b,1) * two_body_dm_ab_diag_act(i_a,i_b) * delta_e(i_a,i_b,ispin,jspin,istate) & + * active_int(i_a,i_b,1) + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + +end + diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index f8782bec..563b4bdf 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -648,10 +648,11 @@ END_PROVIDER double precision :: coef,contrib coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j) psi_in_out_coef(i,j) = coef * hij -! if(vorb == 1.and. iorb == 1)then -! if(vorb == 1.and. iorb == 3)then -! print*, i,hij,coef -! endif + if(orb_i == 5 .and. orb_v == 20)then +! if(orb_i == 2 .and. orb_v == 6 )then + print*, i, ispin + print*, coef * hij,coef,hij + endif norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo enddo @@ -664,6 +665,10 @@ END_PROVIDER norm_no_inv(j,ispin) = norm(j,ispin) one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin) norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) + if(orb_i == 5 .and. orb_v == 20)then +! if(orb_i == 2 .and. orb_v == 6 )then + print*,ispin ,norm(j,ispin) + endif endif enddo do i = 1, N_det_ref @@ -679,21 +684,32 @@ END_PROVIDER do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then -! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) + if(orb_i == 5 .and. orb_v == 20)then +! if(orb_i == 2 .and. orb_v == 6 )then + print*, ispin, energy_cas_dyall_no_exchange(1) , energies_alpha_beta(state_target, ispin) + print*, ispin, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target, ispin) + endif endif enddo enddo ! ispin do state_target = 1, N_states if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then -! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - & one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - & ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) else one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0 endif + if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then + print*, orb_i,orb_v + print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) / norm_bis(state_target,1) + print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) / norm_bis(state_target,2) + print*, fock_core_inactive_total_spin_trace(orb_i,1) + print*, fock_virt_total_spin_trace(orb_v,1) + print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target) + endif enddo enddo enddo diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f index 6028d4fa..768abe8c 100644 --- a/plugins/MRPT_Utils/excitations_cas.irp.f +++ b/plugins/MRPT_Utils/excitations_cas.irp.f @@ -774,6 +774,7 @@ subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij) case (0) hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint) +! hij = 0.d0 end select end @@ -798,7 +799,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint) ! alpha - alpha do i = 1, elec_num_tab_local(1) iorb = occ(i,1) - diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) + diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb) do j = i+1, elec_num_tab_local(1) jorb = occ(j,1) diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb) @@ -808,7 +809,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint) ! beta - beta do i = 1, elec_num_tab_local(2) iorb = occ(i,2) - diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) + diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb) !+ fock_operator_active_from_core_inact(iorb,iorb) do j = i+1, elec_num_tab_local(2) jorb = occ(j,2) diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb) @@ -826,6 +827,8 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint) enddo +! return + ! alpha - core-act do i = 1, elec_num_tab_local(1) iorb = occ(i,1) @@ -833,6 +836,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint) jorb = list_core_inact(j) diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb) core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb) +! diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(1) enddo enddo @@ -843,6 +847,7 @@ double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint) jorb = list_core_inact(j) diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb) core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb) +! diag_H_mat_elem_no_elec_check_no_exchange += core_act_exchange(2) enddo enddo diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 1fd8cb03..91bb7a54 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -113,19 +113,20 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip integer :: degree_scalar call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,index_i),degree_scalar,N_int) +! if(degree_scalar == 2)then +! hialpha = 0.d0 +! endif if(dabs(hialpha).le.1.d-20)then do i_state = 1, N_states delta_e(i_state) = 1.d+20 enddo else call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) - - ! !!!!!!!!!!!!! SHIFTED BK - ! double precision :: hjj - ! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj) - ! delta_e(1) = CI_electronic_energy(1) - hjj - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - +! !!!!!!!!!!!!! SHIFTED BK +! double precision :: hjj +! call i_h_j(tq(1,1,i_alpha),tq(1,1,i_alpha),Nint,hjj) +! delta_e(1) = CI_electronic_energy(1) - hjj +! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! endif hij_array(index_i) = hialpha do i_state = 1,N_states @@ -139,6 +140,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call omp_set_lock( psi_ref_bis_lock(index_i) ) do j = 1, idx_alpha(0) index_j = idx_alpha(j) + !!!!!!!!!!!!!!!!!! WARNING TEST + if(index_j .ne. index_i)cycle do i_state=1,N_states ! standard dressing first order delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state) diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 35940404..cc62295f 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -104,23 +104,23 @@ !print*, 'ionic = ',ionic print*, '1h1p = ',accu - ! 1h1p third order - if(do_third_order_1h1p)then - delta_ij_tmp = 0.d0 - call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state) - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo - enddo - second_order_pt_new_1h1p(i_state) = accu(i_state) - enddo - print*, '1h1p(3)',accu - endif + !! 1h1p third order + !if(do_third_order_1h1p)then + ! delta_ij_tmp = 0.d0 + ! call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) + ! accu = 0.d0 + ! do i_state = 1, N_states + ! do i = 1, N_det_ref + ! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state) + ! do j = 1, N_det_ref + ! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) + ! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + ! enddo + ! enddo + ! second_order_pt_new_1h1p(i_state) = accu(i_state) + ! enddo + ! print*, '1h1p(3)',accu + !endif ! 2h delta_ij_tmp = 0.d0 @@ -200,21 +200,21 @@ enddo print*, '2h2p = ',contrib_2h2p(:) - !! 2h2p old fashion - !delta_ij_tmp = 0.d0 - !call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref) - !accu = 0.d0 - !do i_state = 1, N_states - !do i = 1, N_det_ref - ! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state) - ! do j = 1, N_det_ref - ! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) - ! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - ! enddo - !enddo - !second_order_pt_new_2h2p(i_state) = accu(i_state) - !enddo - !print*, '2h2p = ',accu +! ! 2h2p old fashion +! delta_ij_tmp = 0.d0 +! call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref) +! accu = 0.d0 +! do i_state = 1, N_states +! do i = 1, N_det_ref +! write(*,'(1000(F16.10,x))')delta_ij_tmp(i,:,i_state) +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +! enddo +! second_order_pt_new_2h2p(i_state) = accu(i_state) +! enddo +! print*, '2h2p = ',accu ! total diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index 3cfa7154..a007e761 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -129,6 +129,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) integer :: kspin do jdet = 1, idx(0) if(idx(jdet).ne.idet)then +! cycle ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator ! are connected by the presence of the perturbers determinants |det_tmp> aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb} @@ -213,16 +214,18 @@ subroutine give_1h2p_contrib(matrix_1h2p) double precision :: active_int(n_act_orb,2) double precision :: hij,phase !matrix_1h2p = 0.d0 - elec_num_tab_local = 0 do inint = 1, N_int elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1)) elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1)) enddo +!do i = 1, 1 ! First inactive do i = 1, n_inact_orb ! First inactive iorb = list_inact(i) +! do v = 1, 1 do v = 1, n_virt_orb ! First virtual vorb = list_virt(v) +! do r = 1, 1 do r = 1, n_virt_orb ! Second virtual rorb = list_virt(r) ! take all the integral you will need for i,j,r fixed diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 95f7479e..9d0d1fc6 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -14,6 +14,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] psi_active(j,1,i) = iand(psi_ref(j,1,i),cas_bitmask(j,1,1)) psi_active(j,2,i) = iand(psi_ref(j,2,i),cas_bitmask(j,1,1)) enddo + call debug_det(psi_active(1,1,i),N_int) enddo END_PROVIDER diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 118bbdf7..6bafa287 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -15,6 +15,57 @@ enddo END_PROVIDER + BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] + implicit none + integer :: i,j,ispin,istate + ispin = 1 + do istate = 1, N_states + do j = 1, mo_tot_num + do i = 1, mo_tot_num + one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_alpha(i,j,istate) + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do j = 1, mo_tot_num + do i = 1, mo_tot_num + one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_beta(i,j,istate) + enddo + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [ double precision, one_body_dm_dagger_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] + implicit none + integer :: i,j,ispin,istate + ispin = 1 + do istate = 1, N_states + do j = 1, mo_tot_num + one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_alpha(j,j,istate) + do i = j+1, mo_tot_num + one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate) + one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate) + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do j = 1, mo_tot_num + one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_beta(j,j,istate) + do i = j+1, mo_tot_num + one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_beta(i,j,istate) + one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_beta(i,j,istate) + enddo + enddo + enddo + + END_PROVIDER + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num,N_states) ] &BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num,N_states) ] implicit none diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f index aa8f630b..bb1a341e 100644 --- a/src/Determinants/two_body_dm_map.irp.f +++ b/src/Determinants/two_body_dm_map.irp.f @@ -194,6 +194,8 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl) end BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_act, (n_act_orb, n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_aa_diag_act, (n_act_orb, n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_bb_diag_act, (n_act_orb, n_act_orb)] &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_inact, (n_inact_orb_allocate, n_inact_orb_allocate)] &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_core, (n_core_orb_allocate, n_core_orb_allocate)] &BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_all, (mo_tot_num, mo_tot_num)] @@ -234,6 +236,8 @@ end two_body_dm_ab_diag_all = 0.d0 two_body_dm_ab_diag_act = 0.d0 + two_body_dm_aa_diag_act = 0.d0 + two_body_dm_bb_diag_act = 0.d0 two_body_dm_ab_diag_core = 0.d0 two_body_dm_ab_diag_inact = 0.d0 two_body_dm_diag_core_a_act_b = 0.d0 @@ -269,8 +273,20 @@ end two_body_dm_ab_diag_act(k,m) += 0.5d0 * contrib two_body_dm_ab_diag_act(m,k) += 0.5d0 * contrib enddo + do l = 1, n_occ_ab_act(2) + m = list_act_reverse(occ_act(l,2)) + two_body_dm_bb_diag_act(k,m) += 0.5d0 * contrib + two_body_dm_bb_diag_act(m,k) += 0.5d0 * contrib + enddo + enddo + do j = 1,n_occ_ab_act(1) + k = list_act_reverse(occ_act(j,1)) + do l = 1, n_occ_ab_act(1) + m = list_act_reverse(occ_act(l,1)) + two_body_dm_aa_diag_act(k,m) += 0.5d0 * contrib + two_body_dm_aa_diag_act(m,k) += 0.5d0 * contrib + enddo enddo - ! CORE PART of the diagonal part of the two body dm do j = 1, N_int key_tmp_core(j,1) = psi_det(j,1,i) @@ -325,6 +341,8 @@ end END_PROVIDER BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_aa_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_bb_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] &BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_core_act, (n_core_orb_allocate,n_act_orb,n_act_orb)] implicit none use bitmasks @@ -394,14 +412,22 @@ END_PROVIDER call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ - if(s1==s2)cycle ! Only the alpha/beta two body density matrix ! * c_I * c_J h1 = list_act_reverse(h1) h2 = list_act_reverse(h2) p1 = list_act_reverse(p1) p2 = list_act_reverse(p2) - call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) - + if(s1==s2)then + if(s1==1)then + call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) +! call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,p2,h2,p1) + else + call insert_into_two_body_dm_big_array( two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) +! call insert_into_two_body_dm_big_array( two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,p2,h2,p1) + endif + else ! alpha/beta two body density matrix + call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) + endif else if(degree==1)then! case of the SINGLE EXCITATIONS *************************************************** print*,'h1 = ',h1 h1 = list_act_reverse(h1) @@ -417,6 +443,12 @@ END_PROVIDER ! * c_I * c_J call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) enddo + do k = 1, n_occ_ab(1) + m = list_act_reverse(occ(k,1)) + ! * c_I * c_J + call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) +! call insert_into_two_body_dm_big_array( two_body_dm_aa_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,m,p1,m) + enddo ! core <-> active part of the extra diagonal two body dm do k = 1, n_occ_ab_core(2) @@ -432,6 +464,12 @@ END_PROVIDER ! * c_I * c_J call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) enddo + do k = 1, n_occ_ab(2) + m = list_act_reverse(occ(k,2)) + ! * c_I * c_J + call insert_into_two_body_dm_big_array(two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) +! call insert_into_two_body_dm_big_array(two_body_dm_bb_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,-contrib,h1,m,p1,m) + enddo ! core <-> active part of the extra diagonal two body dm do k = 1, n_occ_ab_core(1) @@ -464,156 +502,3 @@ subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contr end - -double precision function compute_extra_diag_two_body_dm_ab(r1,r2) - implicit none - BEGIN_DOC -! compute the extra diagonal contribution to the alpha/bet two body density at r1, r2 - END_DOC - double precision :: r1(3), r2(3) - double precision :: compute_extra_diag_two_body_dm_ab_act,compute_extra_diag_two_body_dm_ab_core_act - compute_extra_diag_two_body_dm_ab = compute_extra_diag_two_body_dm_ab_act(r1,r2)+compute_extra_diag_two_body_dm_ab_core_act(r1,r2) -end - -double precision function compute_extra_diag_two_body_dm_ab_act(r1,r2) - implicit none - BEGIN_DOC -! compute the extra diagonal contribution to the two body density at r1, r2 -! involving ONLY THE ACTIVE PART, which means that the four index of the excitations -! involved in the two body density matrix are ACTIVE - END_DOC - PROVIDE n_act_orb - double precision, intent(in) :: r1(3),r2(3) - integer :: i,j,k,l - double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb) - double precision :: contrib - double precision :: contrib_tmp -!print*,'n_act_orb = ',n_act_orb - compute_extra_diag_two_body_dm_ab_act = 0.d0 - call give_all_act_mos_at_r(r1,mos_array_r1) - call give_all_act_mos_at_r(r2,mos_array_r2) - do l = 1, n_act_orb ! p2 - do k = 1, n_act_orb ! h2 - do j = 1, n_act_orb ! p1 - do i = 1,n_act_orb ! h1 - contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l) - compute_extra_diag_two_body_dm_ab_act += two_body_dm_ab_big_array_act(i,j,k,l) * contrib_tmp - enddo - enddo - enddo - enddo - -end - -double precision function compute_extra_diag_two_body_dm_ab_core_act(r1,r2) - implicit none - BEGIN_DOC -! compute the extra diagonal contribution to the two body density at r1, r2 -! involving ONLY THE ACTIVE PART, which means that the four index of the excitations -! involved in the two body density matrix are ACTIVE - END_DOC - double precision, intent(in) :: r1(3),r2(3) - integer :: i,j,k,l - double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb) - double precision :: mos_array_core_r1(n_core_orb),mos_array_core_r2(n_core_orb) - double precision :: contrib_core_1,contrib_core_2 - double precision :: contrib_act_1,contrib_act_2 - double precision :: contrib_tmp - compute_extra_diag_two_body_dm_ab_core_act = 0.d0 - call give_all_act_mos_at_r(r1,mos_array_act_r1) - call give_all_act_mos_at_r(r2,mos_array_act_r2) - call give_all_core_mos_at_r(r1,mos_array_core_r1) - call give_all_core_mos_at_r(r2,mos_array_core_r2) - do i = 1, n_act_orb ! h1 - do j = 1, n_act_orb ! p1 - contrib_act_1 = mos_array_act_r1(i) * mos_array_act_r1(j) - contrib_act_2 = mos_array_act_r2(i) * mos_array_act_r2(j) - do k = 1,n_core_orb ! h2 - contrib_core_1 = mos_array_core_r1(k) * mos_array_core_r1(k) - contrib_core_2 = mos_array_core_r2(k) * mos_array_core_r2(k) - contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_core_2 + contrib_act_2 * contrib_core_1) - compute_extra_diag_two_body_dm_ab_core_act += two_body_dm_ab_big_array_core_act(k,i,j) * contrib_tmp - enddo - enddo - enddo - -end - -double precision function compute_diag_two_body_dm_ab_core(r1,r2) - implicit none - double precision :: r1(3),r2(3) - integer :: i,j,k,l - double precision :: mos_array_r1(n_core_orb_allocate),mos_array_r2(n_core_orb_allocate) - double precision :: contrib,contrib_tmp - compute_diag_two_body_dm_ab_core = 0.d0 - call give_all_core_mos_at_r(r1,mos_array_r1) - call give_all_core_mos_at_r(r2,mos_array_r2) - do l = 1, n_core_orb ! - contrib = mos_array_r2(l)*mos_array_r2(l) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - do k = 1, n_core_orb ! - contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - compute_diag_two_body_dm_ab_core += two_body_dm_ab_diag_core(k,l) * contrib_tmp - enddo - enddo - -end - - -double precision function compute_diag_two_body_dm_ab_act(r1,r2) - implicit none - double precision :: r1(3),r2(3) - integer :: i,j,k,l - double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb) - double precision :: contrib,contrib_tmp - compute_diag_two_body_dm_ab_act = 0.d0 - call give_all_act_mos_at_r(r1,mos_array_r1) - call give_all_act_mos_at_r(r2,mos_array_r2) - do l = 1, n_act_orb ! - contrib = mos_array_r2(l)*mos_array_r2(l) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - do k = 1, n_act_orb ! - contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - compute_diag_two_body_dm_ab_act += two_body_dm_ab_diag_act(k,l) * contrib_tmp - enddo - enddo -end - -double precision function compute_diag_two_body_dm_ab_core_act(r1,r2) - implicit none - double precision :: r1(3),r2(3) - integer :: i,j,k,l - double precision :: mos_array_core_r1(n_core_orb_allocate),mos_array_core_r2(n_core_orb_allocate) - double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb) - double precision :: contrib_core_1,contrib_core_2 - double precision :: contrib_act_1,contrib_act_2 - double precision :: contrib_tmp - compute_diag_two_body_dm_ab_core_act = 0.d0 - call give_all_act_mos_at_r(r1,mos_array_act_r1) - call give_all_act_mos_at_r(r2,mos_array_act_r2) - call give_all_core_mos_at_r(r1,mos_array_core_r1) - call give_all_core_mos_at_r(r2,mos_array_core_r2) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - do k = 1, n_act_orb ! - contrib_act_1 = mos_array_act_r1(k) * mos_array_act_r1(k) - contrib_act_2 = mos_array_act_r2(k) * mos_array_act_r2(k) - contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_act_2 + contrib_act_2 * contrib_act_1) -! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - do l = 1, n_core_orb ! - contrib_core_1 = mos_array_core_r1(l) * mos_array_core_r1(l) - contrib_core_2 = mos_array_core_r2(l) * mos_array_core_r2(l) - compute_diag_two_body_dm_ab_core_act += two_body_dm_diag_core_act(l,k) * contrib_tmp - enddo - enddo -end - -double precision function compute_diag_two_body_dm_ab(r1,r2) - implicit none - double precision,intent(in) :: r1(3),r2(3) - double precision :: compute_diag_two_body_dm_ab_act,compute_diag_two_body_dm_ab_core - double precision :: compute_diag_two_body_dm_ab_core_act - compute_diag_two_body_dm_ab = compute_diag_two_body_dm_ab_act(r1,r2)+compute_diag_two_body_dm_ab_core(r1,r2) & - + compute_diag_two_body_dm_ab_core_act(r1,r2) -end