From eda249e631c50e058d429e6ac555fb406a949de0 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 5 Dec 2016 15:10:53 +0100 Subject: [PATCH] final version of MRPT, at least I hope --- plugins/MRPT/MRPT_Utils.main.irp.f | 6 +- plugins/MRPT/print_1h2p.irp.f | 10 +- plugins/MRPT_Utils/H_apply.irp.f | 6 + plugins/MRPT_Utils/energies_cas.irp.f | 188 +++++++++----- plugins/MRPT_Utils/mrpt_dress.irp.f | 64 +++-- plugins/MRPT_Utils/mrpt_utils.irp.f | 303 +++++++++++++---------- plugins/MRPT_Utils/new_way.irp.f | 140 +++++------ plugins/MRPT_Utils/psi_active_prov.irp.f | 88 +++---- 8 files changed, 485 insertions(+), 320 deletions(-) diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/MRPT_Utils.main.irp.f index 72750f8c..e5d925a3 100644 --- a/plugins/MRPT/MRPT_Utils.main.irp.f +++ b/plugins/MRPT/MRPT_Utils.main.irp.f @@ -10,7 +10,7 @@ end subroutine routine_3 implicit none - integer :: i + integer :: i,j !provide fock_virt_total_spin_trace provide delta_ij @@ -23,6 +23,10 @@ subroutine routine_3 write(*,'(A12,X,I3,A3,XX,F16.09)') ' E+PT2 ', i,' = ', CI_energy(i)+second_order_pt_new(i) write(*,'(A12,X,I3,A3,XX,F16.09)') ' E dressed ', i,' = ', CI_dressed_pt2_new_energy(i) write(*,'(A12,X,I3,A3,XX,F16.09)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i) + print*,'coef before and after' + do j = 1, N_det_ref + print*,psi_ref_coef(j,i),CI_dressed_pt2_new_eigenvectors(j,i) + enddo enddo end diff --git a/plugins/MRPT/print_1h2p.irp.f b/plugins/MRPT/print_1h2p.irp.f index a3500e49..2739340b 100644 --- a/plugins/MRPT/print_1h2p.irp.f +++ b/plugins/MRPT/print_1h2p.irp.f @@ -2,16 +2,18 @@ program print_1h2p implicit none read_wf = .True. touch read_wf - call routine + call routine_2 end subroutine routine_2 implicit none - integer :: i,j + integer :: i,j,degree + double precision :: hij +!provide one_creat_virt do i =1, n_act_orb -!do i =1, 2 - write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(i,:,:,:,1) + write(*,'(I3,x,100(F16.10,X))')i,one_creat(i,:,1) ! write(*,'(I3,x,100(F16.10,X))')i,one_anhil_one_creat(1,4,1,2,1) +! enddo diff --git a/plugins/MRPT_Utils/H_apply.irp.f b/plugins/MRPT_Utils/H_apply.irp.f index 56f8a0c7..a7adc480 100644 --- a/plugins/MRPT_Utils/H_apply.irp.f +++ b/plugins/MRPT_Utils/H_apply.irp.f @@ -44,6 +44,7 @@ print s s = H_apply("mrpt_1p") s.filter_only_1p() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -85,6 +86,7 @@ print s s = H_apply("mrpt_2p") s.filter_only_2p() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -105,6 +107,7 @@ print s s = H_apply("mrpt_2h") s.filter_only_2h() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -126,6 +129,7 @@ print s s = H_apply("mrpt_1h2p") s.filter_only_1h2p() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -146,6 +150,7 @@ print s s = H_apply("mrpt_2h1p") s.filter_only_2h1p() +s.unset_skip() s.data["parameters"] = ", delta_ij_, Ndet" s.data["declarations"] += """ integer, intent(in) :: Ndet @@ -166,6 +171,7 @@ print s s = H_apply("mrpt_2h2p") s.filter_only_2h2p() +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 02ff8302..8f6a7eb6 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -321,8 +321,8 @@ 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,N_states)] - implicit none +BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] +implicit none integer :: i,j integer :: ispin,jspin,kspin integer :: orb_i, hole_particle_i,spin_exc_i @@ -332,80 +332,142 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb 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) + enddo + enddo - ! 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 - 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) + do state_target = 1, 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 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 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) enddo + enddo enddo + enddo enddo + enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER +!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 +!integer :: orb_i, hole_particle_i,spin_exc_i +!integer :: orb_j, hole_particle_j,spin_exc_j +!integer :: orb_k, hole_particle_k,spin_exc_k +!double precision :: norm_out(N_states) +!integer(bit_kind), allocatable :: psi_in_out(:,:,:) +!double precision, allocatable :: psi_in_out_coef(:,:) +!use bitmasks +!allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) + +!integer :: iorb,jorb +!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 +! orb_i = list_act(iorb) +! hole_particle_i = 1 +! do jorb = 1, n_act_orb +! orb_j = list_act(jorb) +! hole_particle_j = 1 +! do korb = 1, n_act_orb +! orb_k = list_act(korb) +! hole_particle_k = -1 + +! ! 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 +! 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 +!deallocate(psi_in_out,psi_in_out_coef) + +!END_PROVIDER + BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j @@ -767,6 +829,9 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State norm_bis = 0.d0 do ispin = 1,2 do i = 1, n_det_ref +! if(orb_a == 6 .and. orb_v == 12)then +! print*, 'i ref = ',i +! endif do j = 1, N_int psi_in_out(j,1,i) = psi_ref(j,1,i) psi_in_out(j,2,i) = psi_ref(j,2,i) @@ -778,15 +843,25 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State enddo else call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij) + if(orb_a == 6 .and. orb_v == 12)then + call debug_det(psi_ref(1,1,i),N_int) + call debug_det(psi_in_out(1,1,i),N_int) + print*, hij + endif do j = 1, n_states - double precision :: coef,contrib - coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j) - psi_in_out_coef(i,j) = sign(coef,psi_ref_coef(i,j)) * hij + double precision :: contrib + psi_in_out_coef(i,j) = psi_ref_coef(i,j) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) + !if(orb_a == 6 .and. orb_v == 12)then + ! print*, j,psi_ref_coef(i,j),psi_in_out_coef(i,j) + !endif enddo endif enddo do j = 1, N_states +! if(orb_a == 6 .and. orb_v == 12)then +! print*, 'norm',norm(j,ispin) +! endif if (dabs(norm(j,ispin)) .le. thresh_norm)then norm(j,ispin) = 0.d0 norm_no_inv(j,ispin) = norm(j,ispin) @@ -822,6 +897,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) else +! one_creat_virt(aorb,vorb,state_target) = 0.5d0 * (one_anhil(aorb, 1,state_target) + one_anhil(aorb, 2,state_target) ) one_creat_virt(aorb,vorb,state_target) = 0.d0 endif ! print*, '********' diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index c50e4221..f241d35a 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -66,21 +66,47 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip end if + double precision :: coef_array(N_states) do i_alpha=1,N_tq +! do i = 1, N_det_ref +! do i_state = 1, N_states +! coef_array(i_state) = psi_ref_coef(i,i_state) +! enddo +! call i_H_j(psi_ref(1,1,i),tq(1,1,i_alpha),n_int,hialpha) +! 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,i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) +! endif +! hij_array(i) = hialpha +! do i_state = 1,N_states +! delta_e_inv_array(i,i_state) = 1.d0/delta_e(i_state) +! enddo +! enddo +! do i = 1, N_det_ref +! do j = 1, N_det_ref +! do i_state = 1, N_states +! delta_ij_(i,j,i_state) += hij_array(i) * hij_array(j)* delta_e_inv_array(j,i_state) +! enddo +! enddo +! enddo +! cycle + + + + ! call get_excitation_degree_vector(psi_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_ref,idx_alpha) call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) do j=1,idx_alpha(0) idx_alpha(j) = idx_miniList(idx_alpha(j)) enddo -! double precision :: ihpsi0,coef_pert -! ihpsi0 = 0.d0 -! coef_pert = 0.d0 phase_array =0.d0 do i = 1,idx_alpha(0) index_i = idx_alpha(i) call i_h_j(tq(1,1,i_alpha),psi_ref(1,1,index_i),Nint,hialpha) - double precision :: coef_array(N_states) do i_state = 1, N_states coef_array(i_state) = psi_ref_coef(index_i,i_state) enddo @@ -91,25 +117,21 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip do i_state = 1, N_states delta_e(i_state) = 1.d+20 enddo - else + !else if(degree_scalar== 1)then + 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 + !if(dabs(delta_e(2)) .le. dabs(0.01d0))then + !print*, delta_e(2) + !call debug_det(psi_ref(1,1,index_i),N_int) + !call debug_det(tq(1,1,i_alpha),N_int) + !endif + + !else + !do i_state = 1, N_states + ! delta_e(i_state) = 1.d+20 + !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) do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) enddo @@ -122,7 +144,7 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip do j = 1, idx_alpha(0) index_j = idx_alpha(j) do i_state=1,N_states -! standard dressing first order + ! 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) enddo enddo diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 39cf46db..09efc536 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -30,6 +30,7 @@ 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) @@ -39,141 +40,174 @@ 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 - - ! 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) + ! 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 + 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 - 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 - - ! 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 - - ! 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 - - 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(:) - + 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) + 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 = ',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 + + ! 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 + 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_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 + 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_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 + 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_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 + 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_2h1p(i_state) = accu(i_state) + enddo + print*, '2h1p = ',accu + + ! 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 - accu = 0.d0 + accu = 0.d0 + print*, 'naked matrix' + double precision, allocatable :: hmatrix(:,:) + double precision:: hij,h00 + allocate(hmatrix(N_det_ref, N_det_ref)) + call i_h_j(psi_ref(1,1,1),psi_ref(1,1,1),N_int,h00) + do i = 1, N_det_ref + do j = 1, N_det_Ref + call i_h_j(psi_ref(1,1,i),psi_ref(1,1,j),N_int,hij) + hmatrix(i,j) = hij + enddo + print*, hmatrix(i,i), h00 + hmatrix(i,i) += - h00 + enddo + do i = 1, N_det_ref + write(*,'(1000(F16.10,x))')hmatrix(i,:) + enddo + print*, '' + print*, '' + print*, '' do i_state = 1, N_states print*,'state ',i_state do i = 1, N_det_ref write(*,'(1000(F16.10,x))')delta_ij(i,:,i_state) - do j = i , N_det_ref + do j = 1 , N_det_ref accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) + hmatrix(i,j) += delta_ij(j,i,i_state) enddo enddo second_order_pt_new(i_state) = accu(i_state) print*, 'total= ',accu(i_state) + + do i = 1, N_det_ref + write(*,'(1000(F16.10,x))')hmatrix(i,:) + enddo + enddo + deallocate(hmatrix) @@ -206,7 +240,7 @@ END_PROVIDER call i_h_j(psi_ref(1,1,j),psi_ref(1,1,i),N_int,hij) Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = hij & + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) ) - Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) +! Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo enddo enddo @@ -260,8 +294,8 @@ END_PROVIDER allocate (hmatrix_tmp(N_det_ref,N_det_ref)) allocate (iorder(N_det_ref)) allocate (psi_tmp(N_det_ref)) - print*,'' - print*,'***************************' + print*,'' + print*,'***************************' do i_state = 1, N_states !! Big loop over states print*,'' print*,'Diagonalizing with the dressing for state',i_state @@ -305,7 +339,26 @@ END_PROVIDER call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det_ref,psi_det,N_int,1,N_det_ref) print*,'S^2 = ', CI_dressed_pt2_new_eigenvectors_s2(i_state) enddo + !else if(state_average)then + ! print*,'' + ! print*,'***************************' + ! print*,'' + ! print*,'Doing state average dressings' + ! allocate (hmatrix_tmp(N_det_ref,N_det_ref)) + ! hmatrix_tmp = 0.d0 + ! do i_state = 1, N_states !! Big loop over states + ! do i = 1, N_det_ref + ! do j = 1, N_det_ref + ! hmatrix_tmp(j,i) += Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) + ! enddo + ! enddo + ! enddo + + + ! deallocate(hmatrix_tmp) + else + call lapack_diag(eigenvalues,eigenvectors, & Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det_ref,N_det_ref) CI_electronic_dressed_pt2_new_energy(:) = 0.d0 diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index a4bbe93a..dc921551 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -1,7 +1,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) use bitmasks implicit none - double precision , intent(inout) :: matrix_2h1p(N_det,N_det,*) + double precision , intent(inout) :: matrix_2h1p(N_det_ref,N_det_ref,*) integer :: i,j,r,a,b integer :: iorb, jorb, rorb, aorb, borb integer :: ispin,jspin @@ -22,8 +22,8 @@ subroutine give_2h1p_contrib(matrix_2h1p) elec_num_tab_local = 0 do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + 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, n_inact_orb ! First inactive iorb = list_inact(i) @@ -38,14 +38,14 @@ subroutine give_2h1p_contrib(matrix_2h1p) active_int(a,2) = get_mo_bielec_integral(iorb,jorb,aorb,rorb,mo_integrals_map) ! exchange enddo - integer :: degree(N_det_Ref) - integer :: idx(0:N_det_Ref) + integer :: degree(N_det_ref) + integer :: idx(0:N_det_ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate - integer :: index_orb_act_mono(N_det,3) + integer :: index_orb_act_mono(N_det_ref,3) - do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do idet = 1, N_det_ref + call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det_ref,idx) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a) @@ -53,8 +53,8 @@ subroutine give_2h1p_contrib(matrix_2h1p) do a = 1, n_act_orb ! First active aorb = list_act(a) do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) + det_tmp(inint,1) = psi_ref(inint,1,idet) + det_tmp(inint,2) = psi_ref(inint,2,idet) enddo ! Do the excitation inactive -- > virtual call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin @@ -64,7 +64,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin - ! Check if the excitation is possible or not on psi_det(idet) + ! Check if the excitation is possible or not on psi_ref(idet) accu_elec= 0 do inint = 1, N_int accu_elec+= popcnt(det_tmp(inint,jspin)) @@ -81,7 +81,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1) perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2) enddo - call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int) perturb_dets_phase(a,jspin,ispin) = phase do istate = 1, N_states delta_e(a,jspin,istate) = one_creat(a,jspin,istate) & @@ -109,7 +109,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) !!!!!!!!!!!!!!!!!!!!!!!!!!!! do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int) if (exc(0,1,1) == 1) then ! Mono alpha index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_a @@ -150,7 +150,7 @@ subroutine give_2h1p_contrib(matrix_2h1p) ! you determine the interaction between the excited determinant and the other parent | Jdet > ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet > ! hja = < det_tmp | H | Jdet > - call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int) if(kspin == ispin)then hja = phase * (active_int(borb,2) - active_int(borb,1) ) else @@ -195,7 +195,7 @@ end subroutine give_1h2p_contrib(matrix_1h2p) use bitmasks implicit none - double precision , intent(inout) :: matrix_1h2p(N_det,N_det,*) + double precision , intent(inout) :: matrix_1h2p(N_det_ref,N_det_ref,*) integer :: i,v,r,a,b integer :: iorb, vorb, rorb, aorb, borb integer :: ispin,jspin @@ -216,8 +216,8 @@ subroutine give_1h2p_contrib(matrix_1h2p) elec_num_tab_local = 0 do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + 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, n_inact_orb ! First inactive iorb = list_inact(i) @@ -232,14 +232,14 @@ subroutine give_1h2p_contrib(matrix_1h2p) active_int(a,2) = get_mo_bielec_integral(iorb,aorb,vorb,rorb,mo_integrals_map) ! exchange enddo - integer :: degree(N_det_Ref) - integer :: idx(0:N_det_Ref) + integer :: degree(N_det_ref) + integer :: idx(0:N_det_ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate - integer :: index_orb_act_mono(N_det,3) + integer :: index_orb_act_mono(N_det_ref,3) - do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do idet = 1, N_det_ref + call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det_ref,idx) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb) do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb) @@ -247,8 +247,8 @@ subroutine give_1h2p_contrib(matrix_1h2p) aorb = list_act(a) if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) + det_tmp(inint,1) = psi_ref(inint,1,idet) + det_tmp(inint,2) = psi_ref(inint,2,idet) enddo ! Do the excitation inactive -- > virtual call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin @@ -258,7 +258,7 @@ subroutine give_1h2p_contrib(matrix_1h2p) call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin - ! Check if the excitation is possible or not on psi_det(idet) + ! Check if the excitation is possible or not on psi_ref(idet) accu_elec= 0 do inint = 1, N_int accu_elec+= popcnt(det_tmp(inint,jspin)) @@ -280,7 +280,7 @@ subroutine give_1h2p_contrib(matrix_1h2p) det_tmp(inint,2) = perturb_dets(inint,2,a,jspin,ispin) enddo - call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idet),det_tmp,exc,phase,N_int) perturb_dets_phase(a,jspin,ispin) = phase do istate = 1, N_states delta_e(a,jspin,istate) = one_anhil(a,jspin,istate) & @@ -308,7 +308,7 @@ subroutine give_1h2p_contrib(matrix_1h2p) !!!!!!!!!!!!!!!!!!!!!!!!!!!! do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int) if (exc(0,1,1) == 1) then ! Mono alpha index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a @@ -350,7 +350,7 @@ subroutine give_1h2p_contrib(matrix_1h2p) ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{borb,kspin} a_{iorb,ispin} | Jdet > ! hja = < det_tmp | H | Jdet > - call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int) if(kspin == ispin)then hja = phase * (active_int(borb,1) - active_int(borb,2) ) else @@ -396,7 +396,7 @@ end subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) use bitmasks implicit none - double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*) + double precision , intent(inout) :: matrix_1h1p(N_det_ref,N_det_ref,*) integer :: i,j,r,a,b integer :: iorb, jorb, rorb, aorb, borb,s,sorb integer :: ispin,jspin @@ -413,8 +413,8 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) double precision :: get_mo_bielec_integral double precision :: active_int(n_act_orb,2) double precision :: hij,phase - integer :: degree(N_det_Ref) - integer :: idx(0:N_det_Ref) + integer :: degree(N_det_ref) + integer :: idx(0:N_det_ref) integer :: istate double precision :: hja,delta_e_inact_virt(N_states) integer :: kspin,degree_scalar @@ -422,13 +422,13 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) elec_num_tab_local = 0 do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1)) + elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1)) enddo double precision :: himono,delta_e(N_states),coef_mono(N_states) integer :: state_target - do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do idet = 1, N_det_ref + call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det_ref,idx) do i = 1, n_inact_orb ! First inactive iorb = list_inact(i) do r = 1, n_virt_orb ! First virtual @@ -443,13 +443,13 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) - fock_virt_total_spin_trace(rorb,j) enddo do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) + det_tmp(inint,1) = psi_ref(inint,1,idet) + det_tmp(inint,2) = psi_ref(inint,2,idet) enddo ! Do the excitation inactive -- > virtual call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin - call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono) + call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,himono) do inint = 1, N_int det_pert(inint,1,i,r,ispin) = det_tmp(inint,1) det_pert(inint,2,i,r,ispin) = det_tmp(inint,2) @@ -510,7 +510,7 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) do jdet = 1, idx(0) ! if(idx(jdet).ne.idet)then - call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + call get_mono_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int) if (exc(0,1,1) == 1) then ! Mono alpha aorb = (exc(1,2,1)) !!! a^{\dagger}_a @@ -522,24 +522,24 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) jspin = 2 endif - call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int) + call get_excitation_degree(psi_ref(1,1,idx(jdet)),det_tmp,degree_scalar,N_int) if(degree_scalar .ne. 2)then print*, 'pb !!!' print*, degree_scalar - call debug_det(psi_det(1,1,idx(jdet)),N_int) + call debug_det(psi_ref(1,1,idx(jdet)),N_int) call debug_det(det_tmp,N_int) stop endif - call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idx(jdet)),det_tmp,exc,phase,N_int) double precision :: hij_test hij_test = 0.d0 - call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test) + call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hij_test) do state_target = 1, N_states matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2) enddo else hij_test = 0.d0 - call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij_test) + call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hij_test) do state_target = 1, N_states matrix_1h1p(idet,idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2) enddo @@ -556,7 +556,7 @@ end subroutine give_1p_sec_order_singles_contrib(matrix_1p) use bitmasks implicit none - double precision , intent(inout) :: matrix_1p(N_det,N_det,*) + double precision , intent(inout) :: matrix_1p(N_det_ref,N_det_ref,*) integer :: i,j,r,a,b integer :: iorb, jorb, rorb, aorb, borb,s,sorb integer :: ispin,jspin @@ -572,8 +572,8 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) integer :: accu_elec double precision :: get_mo_bielec_integral double precision :: hij,phase - integer :: degree(N_det_Ref) - integer :: idx(0:N_det_Ref) + integer :: degree(N_det_ref) + integer :: idx(0:N_det_ref) integer :: istate double precision :: hja,delta_e_act_virt(N_states) integer :: kspin,degree_scalar @@ -581,13 +581,13 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) elec_num_tab_local = 0 do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + elec_num_tab_local(1) += popcnt(psi_ref(inint,1,1)) + elec_num_tab_local(2) += popcnt(psi_ref(inint,2,1)) enddo double precision :: himono,delta_e(N_states),coef_mono(N_states) integer :: state_target - do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do idet = 1, N_det_ref + call get_excitation_degree_vector_mono(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det_ref,idx) do i = 1, n_act_orb ! First active iorb = list_act(i) do r = 1, n_virt_orb ! First virtual @@ -601,8 +601,8 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) delta_e_act_virt(j) = - fock_virt_total_spin_trace(rorb,j) enddo do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) + det_tmp(inint,1) = psi_ref(inint,1,idet) + det_tmp(inint,2) = psi_ref(inint,2,idet) enddo ! Do the excitation active -- > virtual call do_mono_excitation(det_tmp,iorb,rorb,ispin,i_ok) @@ -619,7 +619,7 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) enddo cycle endif - call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono) + call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,himono) do inint = 1, N_int det_pert(inint,1,i,r,ispin) = det_tmp(inint,1) det_pert(inint,2,i,r,ispin) = det_tmp(inint,2) @@ -681,10 +681,10 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) det_tmp(inint,1) = det_pert(inint,1,i,r,ispin) det_tmp(inint,2) = det_pert(inint,2,i,r,ispin) enddo - do jdet = 1,N_det + do jdet = 1,N_det_ref double precision :: coef_array(N_states),hij_test - call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono) - call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,coef_array,hij_test,delta_e) + call i_H_j(det_tmp,psi_ref(1,1,jdet),N_int,himono) + call get_delta_e_dyall(psi_ref(1,1,jdet),det_tmp,coef_array,hij_test,delta_e) do state_target = 1, N_states ! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1) matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target) @@ -702,7 +702,7 @@ end subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) use bitmasks implicit none - double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*) + double precision , intent(inout) :: matrix_1h1p(N_det_ref,N_det_ref,*) integer :: i,j,r,a,b integer :: iorb, jorb, rorb, aorb, borb integer :: ispin,jspin @@ -715,8 +715,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) double precision :: get_mo_bielec_integral double precision :: active_int(n_act_orb,2) double precision :: hij,phase - integer :: degree(N_det_Ref) - integer :: idx(0:N_det_Ref) + integer :: degree(N_det_ref) + integer :: idx(0:N_det_ref) integer :: istate double precision :: hja,delta_e_inact_virt(N_states) integer(bit_kind) :: pert_det(N_int,2,n_act_orb,n_act_orb,2) @@ -730,8 +730,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) elec_num_tab_local = 0 do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + 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, n_inact_orb ! First inactive iorb = list_inact(i) @@ -741,8 +741,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(iorb,j) & - fock_virt_total_spin_trace(rorb,j) enddo - do idet = 1, N_det - call get_excitation_degree_vector_double_alpha_beta(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do idet = 1, N_det_ref + call get_excitation_degree_vector_double_alpha_beta(psi_ref,psi_ref(1,1,idet),degree,N_int,N_det_ref,idx) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations do ispin = 1, 2 @@ -752,8 +752,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) do b = 1, n_act_orb borb = list_act(b) do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) + det_tmp(inint,1) = psi_ref(inint,1,idet) + det_tmp(inint,2) = psi_ref(inint,2,idet) enddo ! Do the excitation (i-->a)(ispin) + (b-->r)(other_spin(ispin)) integer :: i_ok,corb,dorb @@ -784,7 +784,7 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) pert_det(inint,2,a,b,ispin) = det_tmp(inint,2) enddo - call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hidouble) + call i_H_j(psi_ref(1,1,idet),det_tmp,N_int,hidouble) do state_target = 1, N_states delta_e(state_target) = one_anhil_one_creat(a,b,ispin,jspin,state_target) + delta_e_inact_virt(state_target) pert_det_coef(a,b,ispin,state_target) = hidouble / delta_e(state_target) @@ -795,7 +795,7 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) enddo do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + call get_double_excitation(psi_ref(1,1,idet),psi_ref(1,1,idx(jdet)),exc,phase,N_int) integer :: c,d,state_target integer(bit_kind) :: det_tmp_bis(N_int,2) ! excitation from I --> J @@ -815,8 +815,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) det_tmp_bis(inint,2) = pert_det(inint,2,c,d,2) enddo double precision :: hjdouble_1,hjdouble_2 - call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1) - call i_H_j(psi_det(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2) + call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1) + call i_H_j(psi_ref(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2) do state_target = 1, N_states matrix_1h1p(idx(jdet),idet,state_target) += (pert_det_coef(c,d,1,state_target) * hjdouble_1 + pert_det_coef(c,d,2,state_target) * hjdouble_2 ) enddo diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 55e8aefb..bd31dc1d 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -380,32 +380,46 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) 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) - + ! first hole + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + ! first particle + kspin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + ! first particle + 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_one_anhil(i_particle_act,j_particle_act,i_hole_act,kspin,jspin,ispin,i_state) + enddo + + +! ! 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 +! 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) +! i_particle_act = particle_list_practical(1,2) ! ! second particle - j_particle_act = particle_list_practical(1,2) - else if (particle_list_practical(1,2) == spin_hole_inact)then +! j_particle_act = particle_list_practical(2,2) +! else if (particle_list_practical(1,2) == spin_hole_inact)then ! ! first particle - i_particle_act = particle_list_practical(1,2) +! i_particle_act = particle_list_practical(2,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 +! j_particle_act = particle_list_practical(1,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 +! 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 @@ -466,6 +480,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) 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) +!write(*,'(100(f16.10,X))'), delta_e_final(2) , delta_e_act(2) , delta_e_inactive(2) , delta_e_virt(2) end @@ -697,31 +712,18 @@ subroutine get_delta_e_dyall_fast(det_1,det_2,delta_e_final) 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 + ! first hole + ispin = hole_list_practical(1,1) + i_hole_act = hole_list_practical(2,1) + ! first particle + kspin = particle_list_practical(1,1) + i_particle_act = particle_list_practical(2,1) + ! first particle + 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_one_anhil(i_particle_act,j_particle_act,i_hole_act,i_state) + delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,kspin,jspin,ispin,i_state) enddo else if (n_holes_act == 3 .and. n_particles_act == 0) then @@ -782,7 +784,7 @@ subroutine get_delta_e_dyall_fast(det_1,det_2,delta_e_final) 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) +!write(*,'(100(f16.10,X))'), delta_e_final(2) , delta_e_act(2) , delta_e_inactive(2) , delta_e_virt(2) end