10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-16 18:25:27 +02:00

Trying do really fo sin free multiple excitations

This commit is contained in:
Emmanuel Giner 2016-11-28 14:51:12 +01:00
parent 8c6bb03a23
commit 38ccfc0cf1
10 changed files with 235 additions and 271 deletions

View File

@ -12,11 +12,6 @@ s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()
print s
s = H_apply("FCI_PT2_new")
s.set_perturbation("decontracted")
s.unset_openmp()
print s
s = H_apply("FCI_no_skip")
s.set_selection_pt2("epstein_nesbet_2x2")

View File

@ -635,34 +635,37 @@ END_PROVIDER
double precision :: coef_mrpt(N_States),coef_array(N_states),hij,delta_e(N_states)
double precision :: hij_array(N_det_Ref),delta_e_array(N_det_ref,N_states)
do i = 1, N_det_non_ref
print*,'i',i
do j = 1, N_det_ref
do k = 1, N_States
coef_array(j) = psi_ref_coef(j,k)
coef_array(k) = psi_ref_coef(j,k)
enddo
call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,i), N_int, Hij_array(j))
call get_delta_e_dyall(psi_ref(1,1,j),psi_non_ref(1,1,i),coef_array,hij_array(j),delta_e)
print*,delta_e(:)
! write(*,'(A7,x,100(F16.10,x))')'delta_e',delta_e(:)
do k = 1, N_states
delta_e_Array(j,k) = delta_e(k)
enddo
enddo
coef_mrpt = 0.d0
do k = 1, N_states
do j = 1, N_det_Ref
coef_mrpt(k) += psi_ref_coef(j,k) * hij_array(j) / delta_e_array(j,k)
enddo
enddo
write(*,'(A7,X,100(F16.10,x))')'coef ',psi_non_ref_coef(i,1) , coef_mrpt(1),psi_non_ref_coef(i,2) , coef_mrpt(2)
do k = 1, N_States
if(dabs(coef_mrpt(k)) .le.1.d-10)then
rho_mrpt(i,k) = 0.d0
exit
endif
print*,k,psi_non_ref_coef(i,k) , coef_mrpt(k)
if(psi_non_ref_coef(i,k) / coef_mrpt(k) .lt.0d0)then
rho_mrpt(i,k) = 1.d0
else
rho_mrpt(i,k) = psi_non_ref_coef(i,k) / coef_mrpt(k)
endif
enddo
print*,'rho',rho_mrpt(i,:)
enddo
END_PROVIDER

View File

@ -2,7 +2,7 @@ program print_1h2p
implicit none
read_wf = .True.
touch read_wf
call routine_2
call routine
end
subroutine routine_2
@ -20,44 +20,20 @@ end
subroutine routine
implicit none
double precision,allocatable :: matrix_1h2p(:,:,:)
allocate (matrix_1h2p(N_det,N_det,N_states))
double precision :: accu(2)
allocate (matrix_1h2p(N_det_ref,N_det_ref,N_states))
integer :: i,j,istate
do i = 1, N_det
do j = 1, N_det
do istate = 1, N_states
matrix_1h2p(i,j,istate) = 0.d0
accu = 0.d0
matrix_1h2p = 0.d0
call H_apply_mrpt_2p(matrix_1h2p,N_det_ref)
do istate = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(istate) += matrix_1h2p(i,j,istate) * psi_coef(i,istate) * psi_coef(j,istate)
enddo
enddo
print*,accu(istate)
enddo
if(.False.)then
call give_1h2p_contrib(matrix_1h2p)
double precision :: accu
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
print*, 'second order ', accu
endif
if(.True.)then
do i = 1, N_det
do j = 1, N_det
do istate = 1, N_states
matrix_1h2p(i,j,istate) = 0.d0
enddo
enddo
enddo
call give_1h2p_new(matrix_1h2p)
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
endif
print*, 'third order ', accu
deallocate (matrix_1h2p)
end

View File

@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
double precision :: norm_out(N_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
use bitmasks
integer :: iorb
@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
orb = list_act(iorb)
hole_particle = 1
spin_exc = ispin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -54,8 +54,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
do state_target = 1,N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
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(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_creat(iorb,ispin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -72,7 +72,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb
integer :: state_target
@ -82,7 +82,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
orb = list_act(iorb)
hole_particle = -1
spin_exc = ispin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -94,8 +94,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
do state_target = 1, N_states
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
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(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
one_anhil(iorb,ispin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -113,7 +113,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
orb_j = list_act(jorb)
hole_particle_j = 1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -142,8 +142,8 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -163,7 +163,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
@ -179,7 +179,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
orb_j = list_act(jorb)
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -192,8 +192,8 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -213,7 +213,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
double precision, allocatable :: psi_in_out_coef(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: state_target
double precision :: energies(n_states)
@ -227,7 +227,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
orb_j = list_act(jorb)
hole_particle_j = -1
spin_exc_j = jspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -289,7 +289,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -306,8 +306,8 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
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(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -330,7 +330,7 @@ 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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -351,7 +351,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -367,8 +367,8 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states)
call u0_H_dyall_u0(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(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
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
@ -391,7 +391,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -412,7 +412,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
orb_k = list_act(korb)
hole_particle_k = 1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -428,8 +428,8 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -452,7 +452,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb
integer :: korb
@ -473,7 +473,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
orb_k = list_act(korb)
hole_particle_k = -1
spin_exc_k = kspin
do i = 1, n_det
do i = 1, n_det_ref
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)
enddo
@ -489,8 +489,8 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target)
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
enddo
enddo
enddo
@ -515,7 +515,7 @@ END_PROVIDER
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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb,i_ok
integer :: state_target
@ -541,7 +541,7 @@ END_PROVIDER
do state_target =1 , N_states
one_anhil_one_creat_inact_virt_norm(iorb,vorb,state_target,ispin) = 0.d0
enddo
do i = 1, n_det
do i = 1, n_det_ref
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)
@ -616,7 +616,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: jorb,i_ok,aorb,orb_a
integer :: state_target
@ -641,7 +641,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
norm = 0.d0
norm_bis = 0.d0
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
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)
@ -715,7 +715,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
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),psi_in_out_coef(n_det_ref,N_states))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states))
integer :: iorb,jorb,i_ok,aorb,orb_a
integer :: state_target
@ -740,7 +740,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
norm = 0.d0
norm_bis = 0.d0
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
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)
@ -825,7 +825,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:)
double precision, allocatable :: delta_e_det(:,:)
use bitmasks
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det_ref+1,N_det_ref+1))
allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det_ref+1,N_det_ref+1))
allocate (eigenvectors(size(H_matrix,1),N_det_ref+1))
allocate (eigenvalues(N_det_ref+1),interact_cas(N_det_ref,N_det_ref))
allocate (delta_e_det(N_det_ref,N_det_ref))
@ -854,7 +854,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from
- fock_virt_total_spin_trace(orb_v,j)
enddo
do ispin = 1,2
do i = 1, n_det
do i = 1, n_det_ref
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)

View File

@ -84,14 +84,21 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
do i_state = 1, N_states
coef_array(i_state) = psi_ref_coef(index_i,i_state)
enddo
if(dabs(hialpha).le.1.d-10)then
integer :: degree_scalar
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,index_i),degree_scalar,N_int)
if(dabs(hialpha).le.1.d-20)then
delta_e = 1.d+20
else
call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
endif
if(degree .ne. 2)then
do i_state = 1, N_states
delta_e(i_state) = 1.d+20
enddo
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)
! phase_array(index_i) = phase
do i_state = 1,N_states
delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state)
enddo
@ -103,16 +110,6 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
call omp_set_lock( psi_ref_bis_lock(index_i) )
do j = 1, idx_alpha(0)
index_j = idx_alpha(j)
! call get_excitation(psi_ref(1,1,index_i),psi_ref(1,1,index_i),exc,degree,phase,N_int)
! if(index_j.ne.index_i)then
! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then
! print*, phase_array(index_j) , phase_array(index_i) ,phase
! call debug_det(psi_ref(1,1,index_i),N_int)
! call debug_det(psi_ref(1,1,index_j),N_int)
! call debug_det(tq(1,1,i_alpha),N_int)
! stop
! endif
! endif
do i_state=1,N_states
! standard dressing first order
delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state)

View File

@ -39,155 +39,141 @@
enddo
print*, '1h = ',accu
! 1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1p(i_state) = accu(i_state)
enddo
print*, '1p = ',accu
!! 1p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_1p(i_state) = accu(i_state)
!enddo
!print*, '1p = ',accu
! 1h1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref)
double precision :: e_corr_from_1h1p_singles(N_states)
!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles)
!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
enddo
print*, '1h1p = ',accu
!! 1h1p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref)
!double precision :: e_corr_from_1h1p_singles(N_states)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_1h1p(i_state) = accu(i_state)
!enddo
!print*, '1h1p = ',accu
! 1h1p third order
if(do_third_order_1h1p)then
delta_ij_tmp = 0.d0
call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
enddo
print*, '1h1p(3)',accu
endif
!! 1h1p third order
!if(do_third_order_1h1p)then
! delta_ij_tmp = 0.d0
! call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
! accu = 0.d0
! do i_state = 1, N_states
! do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
! enddo
! second_order_pt_new_1h1p(i_state) = accu(i_state)
! enddo
! print*, '1h1p(3)',accu
!endif
! 2h
delta_ij_tmp = 0.d0
call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2h(i_state) = accu(i_state)
enddo
print*, '2h = ',accu
!! 2h
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_2h(i_state) = accu(i_state)
!enddo
!print*, '2h = ',accu
! 2p
delta_ij_tmp = 0.d0
call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2p(i_state) = accu(i_state)
enddo
print*, '2p = ',accu
!! 2p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_2p(i_state) = accu(i_state)
!enddo
!print*, '2p = ',accu
! 1h2p
delta_ij_tmp = 0.d0
call give_1h2p_contrib(delta_ij_tmp)
!! 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
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_1h2p(i_state) = accu(i_state)
!enddo
!print*, '1h2p = ',accu
! 2h1p
delta_ij_tmp = 0.d0
call give_2h1p_contrib(delta_ij_tmp)
!!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2h1p(i_state) = accu(i_state)
enddo
print*, '2h1p = ',accu
!! 2h1p
!delta_ij_tmp = 0.d0
!call give_2h1p_contrib(delta_ij_tmp)
!!!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det_ref
! do j = 1, N_det_ref
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_2h1p(i_state) = accu(i_state)
!enddo
!print*, '2h1p = ',accu
! 2h2p
delta_ij_tmp = 0.d0
!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
enddo
enddo
second_order_pt_new_2h2p(i_state) = accu(i_state)
enddo
print*, '2h2p = ',accu
!! 2h2p
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(1)
!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
do i_state = 1, N_states
do i = 1, N_det_ref
! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
do j = i_state, N_det_ref
accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
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
accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state)
enddo
enddo
enddo
second_order_pt_new(i_state) = accu(i_state)
print*, 'total= ',accu(i_state)
second_order_pt_new(i_state) = accu(i_state)
print*, 'total= ',accu(i_state)
enddo
@ -217,7 +203,7 @@ END_PROVIDER
double precision :: hij
do i_state = 1, N_states
do i = 1,N_det_ref
do j = i,N_det_ref
do j = 1,N_det_ref
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) )
@ -284,9 +270,9 @@ END_PROVIDER
do j = 1, N_det_ref
hmatrix_tmp(j,i) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state)
enddo
! print*,i,hmatrix_tmp(i,i)+nuclear_repulsion
enddo
call lapack_diag(eigenvalues,eigenvectors, &
Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det_ref,N_det_ref)
call lapack_diag(eigenvalues,eigenvectors,hmatrix_tmp,N_det_ref,N_det_ref)
write(*,'(A86)')'Looking for the most overlapping state within all eigenvectors of the dressed matrix'
print*,''
print*,'Calculating the overlap for ...'
@ -303,7 +289,10 @@ END_PROVIDER
enddo
print*,''
print*,'Sorting the eigenvectors per overlap'
call dsort(overlap,iorder,n_states)
call dsort(overlap,iorder,n_det_ref)
do j = 1, N_det_ref
print*,overlap(j),iorder(j)
enddo
print*,''
print*,'The most overlapping state is the ',iorder(1)
print*,'with the overlap of ',dabs(overlap(1))

View File

@ -180,7 +180,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
double precision :: delta_e_inactive(N_states)
integer :: i_hole_inact
integer :: i_hole_inact, list_holes_inact(n_inact_orb,2)
call get_excitation_degree(det_1,det_2,degree,N_int)
if(degree>2)then
@ -190,8 +190,13 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list)
delta_e_inactive = 0.d0
integer :: n_holes_total
n_holes_total = 0
do i = 1, n_holes_spin(1)
i_hole_inact = holes_list(i,1)
n_holes_total +=1
list_holes_inact(n_holes_total,1) = i_hole_inact
list_holes_inact(n_holes_total,2) = 1
do i_state = 1, N_states
delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact,i_state)
enddo
@ -199,6 +204,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
do i = 1, n_holes_spin(2)
i_hole_inact = holes_list(i,2)
n_holes_total +=1
list_holes_inact(n_holes_total,1) = i_hole_inact
list_holes_inact(n_holes_total,2) = 2
do i_state = 1, N_states
delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state)
enddo
@ -370,15 +378,41 @@ 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)
! spin_hole_part_act =
if(jspin == spin_hole_inact )then
kspin = spin_hole_part_act
ispin = spin_hole_part_act
else
jspin = spin_hole_part_act
ispin = spin_hole_part_act
endif
! by convention, you first make a movement in the cas
! first hole
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
i_hole_act = hole_list_practical(2,1)
jspin = spin_hole_inact
! first particle
jspin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
! second particle
kspin = particle_list_practical(1,2)
j_particle_act = particle_list_practical(2,2)
call get_excitation_degree(det_1,det_2,degree,N_int)
if(degree == 2)then
print*, ''
call debug_det(det_1,N_int)
call debug_det(det_2,N_int)
call get_excitation(det_1,det_2,exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*, s1,h1,p1
print*, s2,h2,p2
print*, '---'
print*, ispin, i_hole_act
print*, jspin, i_particle_act
print*, kspin, j_particle_act
pause
endif
do i_state = 1, N_states
delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state)

View File

@ -1 +1 @@
Determinants Properties Hartree_Fock Davidson MRPT_Utils
Determinants Properties Hartree_Fock Davidson

View File

@ -46,36 +46,6 @@ subroutine pt2_epstein_nesbet ($arguments)
end
subroutine pt2_decontracted ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem_fock, h
double precision :: i_H_psi_array(N_st)
double precision :: coef_pert
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_pert_new_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array,coef_pert)
H_pert_diag = 0.d0
c_pert(1) = coef_pert
e_2_pert(1) = coef_pert * i_H_psi_array(1)
! print*,coef_pert,i_H_psi_array(1)
end
subroutine pt2_epstein_nesbet_2x2 ($arguments)
use bitmasks
implicit none