From a946fc615b5304a2e59b45e6180636b095b4bb3b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Tue, 29 Nov 2016 16:48:24 +0100 Subject: [PATCH] Beginning to merge MRCC and MRPT --- plugins/MRCC_Utils/mrcc_utils.irp.f | 11 +- plugins/MRPT_Utils/H_apply.irp.f | 2 + plugins/MRPT_Utils/energies_cas.irp.f | 88 ++++-- plugins/MRPT_Utils/mrpt_dress.irp.f | 38 ++- plugins/MRPT_Utils/mrpt_utils.irp.f | 225 +++++++------ plugins/MRPT_Utils/psi_active_prov.irp.f | 385 ++++++++++++++++++++--- 6 files changed, 549 insertions(+), 200 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 92425eb0..9cf0330b 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -634,8 +634,11 @@ END_PROVIDER integer :: i, j, k double precision :: coef_mrpt(N_States),coef_array(N_states),hij,delta_e(N_states) double precision :: hij_array(N_det_Ref),delta_e_array(N_det_ref,N_states) + integer :: number_of_holes, number_of_particles,nh,np do i = 1, N_det_non_ref print*,'i',i + nh = number_of_holes(psi_non_ref(1,1,i)) + np = number_of_particles(psi_non_ref(1,1,i)) do j = 1, N_det_ref do k = 1, N_States coef_array(k) = psi_ref_coef(j,k) @@ -653,7 +656,9 @@ END_PROVIDER coef_mrpt(k) += psi_ref_coef(j,k) * hij_array(j) / delta_e_array(j,k) enddo enddo + write(*,'(A7,X,100(F16.10,x))')'coef ',psi_non_ref_coef(i,1) , coef_mrpt(1),psi_non_ref_coef(i,2) , coef_mrpt(2) + print*, nh,np do k = 1, N_States if(dabs(coef_mrpt(k)) .le.1.d-10)then rho_mrpt(i,k) = 0.d0 @@ -666,6 +671,7 @@ END_PROVIDER endif enddo print*,'rho',rho_mrpt(i,:) + write(33,*)i,rho_mrpt(i,:) enddo END_PROVIDER @@ -1011,7 +1017,7 @@ END_PROVIDER double precision function get_dij_index(II, i, s, Nint) integer, intent(in) :: II, i, s, Nint double precision, external :: get_dij - double precision :: HIi, phase + double precision :: HIi, phase,delta_e_final(N_states) if(lambda_type == 0) then call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) @@ -1026,7 +1032,8 @@ double precision function get_dij_index(II, i, s, Nint) get_dij_index = get_dij_index else if(lambda_type == 3) then call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) - get_dij_index = HIi * rho_mrpt(i, s) + call get_delta_e_dyall_fast(psi_ref(1,1,II),psi_non_ref(1,1,i),delta_e_final) + get_dij_index = HIi * rho_mrpt(i, s) / delta_e_final(s) end if end function diff --git a/plugins/MRPT_Utils/H_apply.irp.f b/plugins/MRPT_Utils/H_apply.irp.f index 6f17ab05..56f8a0c7 100644 --- a/plugins/MRPT_Utils/H_apply.irp.f +++ b/plugins/MRPT_Utils/H_apply.irp.f @@ -23,6 +23,7 @@ print s s = H_apply("mrpt_1h") s.filter_only_1h() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -63,6 +64,7 @@ print s s = H_apply("mrpt_1h1p") s.filter_only_1h1p() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index 33baeb6a..02ff8302 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -188,12 +188,14 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) psi_in_out(j,2,i) = psi_active(j,2,i) enddo enddo - call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) - call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) - 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) - two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + do state_target = 1 , N_states + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + 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) + two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -319,7 +321,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a END_PROVIDER -BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] +BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -336,45 +338,69 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer :: korb integer :: state_target double precision :: energies(n_states) + double precision :: norm_spins(2,2,N_states), energies_spins(2,2,N_states) + double precision :: thresh_norm + thresh_norm = 1.d-10 do iorb = 1,n_act_orb - do ispin = 1,2 orb_i = list_act(iorb) hole_particle_i = 1 - spin_exc_i = ispin do jorb = 1, n_act_orb - do jspin = 1,2 orb_j = list_act(jorb) hole_particle_j = 1 - spin_exc_j = jspin do korb = 1, n_act_orb - do kspin = 1,2 orb_k = list_act(korb) hole_particle_k = -1 - spin_exc_k = kspin - do i = 1, n_det_ref - do j = 1, n_states - psi_in_out_coef(i,j) = psi_ref_coef(i,j) - enddo - do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,2,i) = psi_active(j,2,i) + + ! loop on the spins + ! By definition, orb_i is the particle of spin ispin + ! a^+_{ispin , orb_i} + do ispin = 1, 2 + do jspin = 1, 2 + ! By definition, orb_j and orb_k are the couple of particle/hole of spin jspin + ! a^+_{jspin , orb_j} a_{jspin , orb_k} + ! norm_spins(ispin,jspin) :: norm of the wave function a^+_{ispin , orb_i} a^+_{jspin , orb_j} a_{jspin , orb_k} | Psi > + ! energies_spins(ispin,jspin) :: Dyall energu of the wave function a^+_{ispin , orb_i} a^+_{jspin , orb_j} a_{jspin , orb_k} | Psi > + do i = 1, n_det_ref + do j = 1, n_states + psi_in_out_coef(i,j) = psi_ref_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,2,i) = psi_active(j,2,i) + enddo + enddo + do state_target = 1, N_states + ! hole :: hole_particle_k, jspin + call apply_exc_to_psi(orb_k,hole_particle_k,jspin, & + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call apply_exc_to_psi(orb_j,hole_particle_j,jspin, & + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call apply_exc_to_psi(orb_i,hole_particle_i,ispin, & + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + if(dabs(norm_out(state_target)).lt.thresh_norm)then + norm_spins(ispin,jspin,state_target) = 0.d0 + else + norm_spins(ispin,jspin,state_target) = 1.d0 + endif + 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_spins(ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + enddo enddo enddo + integer :: icount + ! averaging over all possible spin permutations with Heaviside norm do state_target = 1, N_states - call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) - call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) - call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) - 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) - two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) + icount = 0 + do jspin = 1, 2 + do ispin = 1, 2 + icount += 1 + two_creat_one_anhil(iorb,jorb,korb,state_target) = energies_spins(ispin,jspin,state_target) * norm_spins(ispin,jspin,state_target) + enddo + enddo + two_creat_one_anhil(iorb,jorb,korb,state_target) = two_creat_one_anhil(iorb,jorb,korb,state_target) / dble(icount) enddo - enddo enddo - enddo enddo - enddo enddo deallocate(psi_in_out,psi_in_out_coef) diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index eccdae0a..c50e4221 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -48,18 +48,18 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip integer :: exc(0:2,2,2),degree - leng = max(N_det_ref, N_det_ref) + leng = max(N_det_generators, N_det_generators) allocate(miniList(Nint, 2, leng), idx_miniList(leng)) !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) - call create_minilist_find_previous(key_mask, psi_ref, miniList, i_generator-1, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) if(fullMatch) then return end if - call find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + call find_connections_previous(n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) if(N_tq > 0) then call create_minilist(key_mask, psi_ref, miniList, idx_miniList, N_det_ref, N_minilist, Nint) @@ -88,17 +88,28 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,index_i),degree_scalar,N_int) if(dabs(hialpha).le.1.d-20)then - delta_e = 1.d+20 - else - call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) - endif - if(degree .ne. 2)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) + do i_state = 1, N_states + if(isnan(delta_e(i_state)))then + print*, 'i_state',i_state + call debug_det(psi_ref(1,1,index_i),N_int) + call debug_det(tq(1,1,i_alpha),N_int) + print*, delta_e(:) + stop + endif + enddo endif +! if(degree_scalar .ne. 1)then +! do i_state = 1, N_states +! delta_e(i_state) = 1.d+20 +! enddo +! endif hij_array(index_i) = hialpha - call get_excitation(psi_ref(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) +! call get_excitation(psi_ref(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) enddo @@ -134,12 +145,12 @@ end END_PROVIDER -subroutine find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) +subroutine find_connections_previous(n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) use bitmasks implicit none - integer, intent(in) :: i_generator,n_selected, Nint + integer, intent(in) :: n_selected, Nint integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m @@ -180,8 +191,3 @@ subroutine find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N end - - - - - diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index cb007199..39cf46db 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -39,127 +39,126 @@ enddo print*, '1h = ',accu -!! 1p -!delta_ij_tmp = 0.d0 -!call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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_1p(i_state) = accu(i_state) -!enddo -!print*, '1p = ',accu + ! 1p + delta_ij_tmp = 0.d0 + call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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_1p(i_state) = accu(i_state) + enddo + print*, '1p = ',accu -!! 1h1p -!delta_ij_tmp = 0.d0 -!call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref) -!double precision :: e_corr_from_1h1p_singles(N_states) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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 = ',accu + ! 1h1p + delta_ij_tmp = 0.d0 + call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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 = ',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 -! 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 + 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 -!call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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_2h(i_state) = accu(i_state) -!enddo -!print*, '2h = ',accu + ! 2h + delta_ij_tmp = 0.d0 + call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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_2h(i_state) = accu(i_state) + enddo + print*, '2h = ',accu -!! 2p -!delta_ij_tmp = 0.d0 -!call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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_2p(i_state) = accu(i_state) -!enddo -!print*, '2p = ',accu + ! 2p + delta_ij_tmp = 0.d0 + call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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_2p(i_state) = accu(i_state) + enddo + print*, '2p = ',accu -!! 1h2p -!delta_ij_tmp = 0.d0 -!call give_1h2p_contrib(delta_ij_tmp) -!!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det_ref) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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_1h2p(i_state) = accu(i_state) -!enddo -!print*, '1h2p = ',accu + ! 1h2p + delta_ij_tmp = 0.d0 + call give_1h2p_contrib(delta_ij_tmp) + !!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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_1h2p(i_state) = accu(i_state) + enddo + print*, '1h2p = ',accu -!! 2h1p -!delta_ij_tmp = 0.d0 -!call give_2h1p_contrib(delta_ij_tmp) -!!!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref) -!accu = 0.d0 -!do i_state = 1, N_states -!do i = 1, N_det_ref -! 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_2h1p(i_state) = accu(i_state) -!enddo -!print*, '2h1p = ',accu + ! 2h1p + delta_ij_tmp = 0.d0 + call give_2h1p_contrib(delta_ij_tmp) + !!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det_ref + 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_2h1p(i_state) = accu(i_state) + enddo + print*, '2h1p = ',accu -!! 2h2p + ! 2h2p -!double precision :: contrib_2h2p(N_states) -!call give_2h2p(contrib_2h2p) -!do i_state = 1, N_states -! do i = 1, N_det_ref -! delta_ij(i,i,i_state) += contrib_2h2p(i_state) -! enddo -!second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) -!enddo -!print*, '2h2p = ',contrib_2h2p(:) + double precision :: contrib_2h2p(N_states) + call give_2h2p(contrib_2h2p) + do i_state = 1, N_states + do i = 1, N_det_ref + delta_ij(i,i,i_state) += contrib_2h2p(i_state) + enddo + second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) + enddo + print*, '2h2p = ',contrib_2h2p(:) ! total diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index e3b0986a..55e8aefb 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -9,10 +9,10 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] integer :: i,j,k,l provide cas_bitmask !print*, 'psi_active ' - do i = 1, N_det + do i = 1, N_det_ref do j = 1, N_int - psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1)) - psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1)) + 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 enddo END_PROVIDER @@ -184,7 +184,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) call get_excitation_degree(det_1,det_2,degree,N_int) if(degree>2)then - delta_e_final = -1.d+10 + do i_state = 1, N_States + delta_e_final(i_state) = -1.d+10 + enddo return endif @@ -198,7 +200,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) list_holes_inact(n_holes_total,1) = i_hole_inact list_holes_inact(n_holes_total,2) = 1 do i_state = 1, N_states - delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) + delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) enddo enddo @@ -223,14 +225,14 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) do i = 1, n_particles_spin(1) i_part_virt = particles_list(i,1) do i_state = 1, N_states - delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state) + delta_e_virt(i_state) += fock_virt_total_spin_trace(i_part_virt,i_state) enddo enddo do i = 1, n_particles_spin(2) i_part_virt = particles_list(i,2) do i_state = 1, N_states - delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state) + delta_e_virt(i_state) += fock_virt_total_spin_trace(i_part_virt,i_state) enddo enddo @@ -382,40 +384,27 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) ! integer :: spin_hole_inact, spin_hole_part_act spin_hole_inact = list_holes_inact(1,2) -! spin_hole_part_act = - if(jspin == spin_hole_inact )then - kspin = spin_hole_part_act - ispin = spin_hole_part_act - else - jspin = spin_hole_part_act - ispin = spin_hole_part_act - endif - ! by convention, you first make a movement in the cas - ! first hole + +! ! by convention, you first make a movement in the cas +! ! first hole i_hole_act = hole_list_practical(2,1) - jspin = spin_hole_inact - ! first particle - i_particle_act = particle_list_practical(2,1) - ! second particle - j_particle_act = particle_list_practical(2,2) - call get_excitation_degree(det_1,det_2,degree,N_int) - if(degree == 2)then - print*, '' - call debug_det(det_1,N_int) - call debug_det(det_2,N_int) - call get_excitation(det_1,det_2,exc,degree,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - print*, s1,h1,p1 - print*, s2,h2,p2 - print*, '---' - print*, ispin, i_hole_act - print*, jspin, i_particle_act - print*, kspin, j_particle_act - pause + if(particle_list_practical(1,1) == spin_hole_inact)then +! ! first particle + i_particle_act = particle_list_practical(2,2) +! ! second particle + j_particle_act = particle_list_practical(1,2) + else if (particle_list_practical(1,2) == spin_hole_inact)then +! ! first particle + i_particle_act = particle_list_practical(1,2) +! ! second particle + j_particle_act = particle_list_practical(2,2) + else + print*, 'pb in n_holes_act == 1 .and. n_particles_act == 2 !!' + stop endif do i_state = 1, N_states - delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state) + delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,i_state) enddo else if (n_holes_act == 3 .and. n_particles_act == 0) then @@ -464,7 +453,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) enddo endif else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then - delta_e_act = -10000000.d0 + do i = 1, N_states + delta_e_act(i_state) = -10000000.d0 + enddo endif !print*, 'one_anhil_spin_trace' @@ -478,3 +469,321 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) end + +subroutine get_delta_e_dyall_fast(det_1,det_2,delta_e_final) + BEGIN_DOC + ! routine that returns the delta_e with the Moller Plesset and Dyall operators + ! + ! with det_1 being a determinant from the cas, and det_2 being a perturber + ! + ! Delta_e(det_1,det_2) = sum (hole) epsilon(hole) + sum(part) espilon(part) + delta_e(act) + ! + ! where hole is necessary in the inactive, part necessary in the virtuals + ! + ! and delta_e(act) is obtained from the contracted application of the excitation + ! + ! operator in the active space that lead from det_1 to det_2 + END_DOC + implicit none + use bitmasks + double precision, intent(out) :: delta_e_final(N_states) + integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + integer :: i,j,k,l + integer :: i_state + + integer :: n_holes_spin(2) + integer :: n_holes + integer :: holes_list(N_int*bit_kind_size,2) + + + double precision :: delta_e_inactive(N_states) + integer :: i_hole_inact, list_holes_inact(n_inact_orb,2) + + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree>2)then + do i_state = 1, N_States + delta_e_final(i_state) = -1.d+10 + enddo + return + endif + + call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) + delta_e_inactive = 0.d0 + integer :: n_holes_total + n_holes_total = 0 + do i = 1, n_holes_spin(1) + i_hole_inact = holes_list(i,1) + n_holes_total +=1 + list_holes_inact(n_holes_total,1) = i_hole_inact + list_holes_inact(n_holes_total,2) = 1 + do i_state = 1, N_states + delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) + enddo + enddo + + do i = 1, n_holes_spin(2) + i_hole_inact = holes_list(i,2) + n_holes_total +=1 + list_holes_inact(n_holes_total,1) = i_hole_inact + list_holes_inact(n_holes_total,2) = 2 + do i_state = 1, N_states + delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) + enddo + enddo + + double precision :: delta_e_virt(N_states) + integer :: i_part_virt + integer :: n_particles_spin(2) + integer :: n_particles + integer :: particles_list(N_int*bit_kind_size,2) + + call give_particles_in_virt_space(det_2,n_particles_spin,n_particles,particles_list) + delta_e_virt = 0.d0 + do i = 1, n_particles_spin(1) + i_part_virt = particles_list(i,1) + do i_state = 1, N_states + delta_e_virt(i_state) += fock_virt_total_spin_trace(i_part_virt,i_state) + enddo + enddo + + do i = 1, n_particles_spin(2) + i_part_virt = particles_list(i,2) + do i_state = 1, N_states + delta_e_virt(i_state) += fock_virt_total_spin_trace(i_part_virt,i_state) + enddo + enddo + + + integer :: n_holes_spin_act(2),n_particles_spin_act(2) + integer :: n_holes_act,n_particles_act + integer :: holes_active_list(2*n_act_orb,2) + integer :: holes_active_list_spin_traced(4*n_act_orb) + integer :: particles_active_list(2*n_act_orb,2) + integer :: particles_active_list_spin_traced(4*n_act_orb) + double precision :: delta_e_act(N_states) + delta_e_act = 0.d0 + call give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin_act,n_particles_spin_act, & + n_holes_act,n_particles_act,holes_active_list,particles_active_list) + integer :: icount,icountbis + integer :: hole_list_practical(2,elec_num_tab(1)+elec_num_tab(2)), particle_list_practical(2,elec_num_tab(1)+elec_num_tab(2)) + icount = 0 + icountbis = 0 + do i = 1, n_holes_spin_act(1) + icount += 1 + icountbis += 1 + hole_list_practical(1,icountbis) = 1 + hole_list_practical(2,icountbis) = holes_active_list(i,1) + holes_active_list_spin_traced(icount) = holes_active_list(i,1) + enddo + do i = 1, n_holes_spin_act(2) + icount += 1 + icountbis += 1 + hole_list_practical(1,icountbis) = 2 + hole_list_practical(2,icountbis) = holes_active_list(i,2) + holes_active_list_spin_traced(icount) = holes_active_list(i,2) + enddo + if(icount .ne. n_holes_act) then + print*,'' + print*, icount, n_holes_act + print * , 'pb in holes_active_list_spin_traced !!' + stop + endif + + icount = 0 + icountbis = 0 + do i = 1, n_particles_spin_act(1) + icount += 1 + icountbis += 1 + particle_list_practical(1,icountbis) = 1 + particle_list_practical(2,icountbis) = particles_active_list(i,1) + particles_active_list_spin_traced(icount) = particles_active_list(i,1) + enddo + do i = 1, n_particles_spin_act(2) + icount += 1 + icountbis += 1 + particle_list_practical(1,icountbis) = 2 + particle_list_practical(2,icountbis) = particles_active_list(i,2) + particles_active_list_spin_traced(icount) = particles_active_list(i,2) + enddo + if(icount .ne. n_particles_act) then + print*, icount, n_particles_act + print * , 'pb in particles_active_list_spin_traced !!' + stop + endif + + + integer :: i_hole_act, j_hole_act, k_hole_act + integer :: i_particle_act, j_particle_act, k_particle_act + + + integer :: ispin,jspin,kspin + if (n_holes_act == 0 .and. n_particles_act == 1) then + ispin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 1)then + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + i_hole = list_inact_reverse(h1) + i_part = list_act_reverse(p1) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) + enddo + else if (degree == 2)then + do i_state = 1, N_states + delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) + enddo + endif + + else if (n_holes_act == 1 .and. n_particles_act == 0) then + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 1)then + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + i_hole = list_act_reverse(h1) + i_part = list_virt_reverse(p1) + do i_state = 1, N_states + delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state) + enddo + else if (degree == 2)then + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) + enddo + endif + + else if (n_holes_act == 1 .and. n_particles_act == 1) then + ! first hole + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + ! first particle + jspin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin,i_state) + enddo + + else if (n_holes_act == 2 .and. n_particles_act == 0) then + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + jspin = hole_list_practical(1,2) + j_hole_act = hole_list_practical(2,2) + do i_state = 1, N_states + delta_e_act(i_state) += two_anhil(i_hole_act,j_hole_act,ispin,jspin,i_state) + enddo + + else if (n_holes_act == 0 .and. n_particles_act == 2) then + ispin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + jspin = particle_list_practical(1,2) + j_particle_act = particle_list_practical(2,2) + do i_state = 1, N_states + delta_e_act(i_state) += two_creat(i_particle_act,j_particle_act,ispin,jspin,i_state) + enddo + + else if (n_holes_act == 2 .and. n_particles_act == 1) then + ! first hole + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + ! second hole + jspin = hole_list_practical(1,2) + j_hole_act = hole_list_practical(2,2) + ! first particle + kspin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + do i_state = 1, N_states + delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state) + enddo + + else if (n_holes_act == 1 .and. n_particles_act == 2) then + ! First find the particle that has been added from the inactive + ! + integer :: spin_hole_inact, spin_hole_part_act + spin_hole_inact = list_holes_inact(1,2) + +! ! by convention, you first make a movement in the cas +! ! first hole + i_hole_act = hole_list_practical(2,1) + if(particle_list_practical(1,1) == spin_hole_inact)then +! ! first particle + i_particle_act = particle_list_practical(2,2) +! ! second particle + j_particle_act = particle_list_practical(1,2) + else if (particle_list_practical(1,2) == spin_hole_inact)then +! ! first particle + i_particle_act = particle_list_practical(1,2) +! ! second particle + j_particle_act = particle_list_practical(2,2) + else + print*, 'pb in n_holes_act == 1 .and. n_particles_act == 2 !!' + stop + endif + + do i_state = 1, N_states + delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,i_state) + enddo + + else if (n_holes_act == 3 .and. n_particles_act == 0) then + ! first hole + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + ! second hole + jspin = hole_list_practical(1,2) + j_hole_act = hole_list_practical(2,2) + ! third hole + kspin = hole_list_practical(1,3) + k_hole_act = hole_list_practical(2,3) + do i_state = 1, N_states + delta_e_act(i_state) += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin,i_state) + enddo + + else if (n_holes_act == 0 .and. n_particles_act == 3) then + ! first particle + ispin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + ! second particle + jspin = particle_list_practical(1,2) + j_particle_act = particle_list_practical(2,2) + ! second particle + kspin = particle_list_practical(1,3) + k_particle_act = particle_list_practical(2,3) + do i_state = 1, N_states + delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state) + enddo + + else if (n_holes_act .eq. 0 .and. n_particles_act .eq.0)then + integer :: degree + integer(bit_kind) :: det_1_active(N_int,2) + integer :: h1,h2,p1,p2,s1,s2 + integer :: exc(0:2,2,2) + integer :: i_hole, i_part + double precision :: phase + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 1)then + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + i_hole = list_inact_reverse(h1) + i_part = list_virt_reverse(p1) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) + enddo + endif + else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then + do i = 1, N_states + delta_e_act(i_state) = -10000000.d0 + enddo + endif + +!print*, 'one_anhil_spin_trace' +!print*, one_anhil_spin_trace(1), one_anhil_spin_trace(2) + + + do i_state = 1, n_states + delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state) + enddo +!write(*,'(100(f16.10,X))'), delta_e_final(1) , delta_e_act(1) , delta_e_inactive(1) , delta_e_virt(1) + +end + +