LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-29 08:24:51 +02:00
Code Issues Releases Wiki Activity

final version of MRPT, at least I hope

Browse Source
This commit is contained in:
Emmanuel Giner 2016-12-05 15:10:53 +01:00
parent a946fc615b
commit eda249e631
8 changed files with 485 additions and 320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
plugins/MRPT/MRPT_Utils.main.irp.f
View File

@ -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

10
plugins/MRPT/print_1h2p.irp.f
View File

@ -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

6
plugins/MRPT_Utils/H_apply.irp.f
View File

@ -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

188
plugins/MRPT_Utils/energies_cas.irp.f
View File

@ -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*, '********'

64
plugins/MRPT_Utils/mrpt_dress.irp.f
View File

@ -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

303
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -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

140
plugins/MRPT_Utils/new_way.irp.f
View File

@ -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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a_{b} a^{\dagger}_a | Idet>
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)
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a^{\dagger}_b a_{a} | Idet>
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

88
plugins/MRPT_Utils/psi_active_prov.irp.f
View File

@ -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