mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-18 12:03:57 +01:00
1115 lines
40 KiB
Fortran
1115 lines
40 KiB
Fortran
BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
|
|
implicit none
|
|
integer :: i
|
|
double precision :: energies(N_states_diag)
|
|
do i = 1, N_states
|
|
call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
|
|
energy_cas_dyall(i) = energies(i)
|
|
print*, 'energy_cas_dyall(i)', energy_cas_dyall(i)
|
|
enddo
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
|
|
implicit none
|
|
integer :: i
|
|
double precision :: energies(N_states_diag)
|
|
do i = 1, N_states
|
|
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
|
|
energy_cas_dyall_no_exchange(i) = energies(i)
|
|
print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i)
|
|
enddo
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
|
|
implicit none
|
|
integer :: i,j
|
|
integer :: ispin
|
|
integer :: orb, hole_particle,spin_exc
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
use bitmasks
|
|
|
|
integer :: iorb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
do iorb = 1,n_act_orb
|
|
do ispin = 1,2
|
|
orb = list_act(iorb)
|
|
hole_particle = 1
|
|
spin_exc = ispin
|
|
do i = 1, n_det
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
|
|
implicit none
|
|
integer :: i,j
|
|
integer :: ispin
|
|
integer :: orb, hole_particle,spin_exc
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
do iorb = 1,n_act_orb
|
|
do ispin = 1,2
|
|
orb = list_act(iorb)
|
|
hole_particle = -1
|
|
spin_exc = ispin
|
|
do i = 1, n_det
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb,hole_particle,spin_exc, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)]
|
|
implicit none
|
|
integer :: i,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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 i = 1, n_det
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)]
|
|
implicit none
|
|
integer :: i,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: state_target
|
|
state_target = 1
|
|
double precision :: energies(n_states_diag)
|
|
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 i = 1, n_det
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2,N_States)]
|
|
implicit none
|
|
integer :: i,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
integer :: iorb,jorb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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 i = 1, n_det
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
if(orb_i == orb_j .and. ispin .ne. jspin)then
|
|
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
|
|
else
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (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
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: korb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(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,2,2,2,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_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: korb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(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, three_creat, (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
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: korb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(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, three_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
|
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb
|
|
integer :: korb
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
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
|
|
do j = 1, n_states_diag
|
|
psi_in_out_coef(i,j) = psi_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
|
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, &
|
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(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, one_anhil_one_creat_inact_virt, (n_inact_orb,n_virt_orb,N_States)]
|
|
&BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_norm, (n_inact_orb,n_virt_orb,N_States,2)]
|
|
implicit none
|
|
integer :: i,vorb,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i
|
|
integer :: orb_v
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb,i_ok
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
double precision :: hij
|
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2),norm_bis(N_states,2)
|
|
double precision :: energies_alpha_beta(N_states,2)
|
|
|
|
|
|
double precision :: thresh_norm
|
|
|
|
thresh_norm = 1.d-10
|
|
|
|
|
|
|
|
do vorb = 1,n_virt_orb
|
|
orb_v = list_virt(vorb)
|
|
do iorb = 1, n_inact_orb
|
|
orb_i = list_inact(iorb)
|
|
norm = 0.d0
|
|
norm_bis = 0.d0
|
|
do ispin = 1,2
|
|
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 j = 1, N_int
|
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
|
enddo
|
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
|
if(i_ok.ne.1)then
|
|
print*, orb_i,orb_v
|
|
call debug_det(psi_in_out,N_int)
|
|
print*, 'pb, i_ok ne 0 !!!'
|
|
endif
|
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
|
do j = 1, n_states
|
|
double precision :: coef,contrib
|
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
|
enddo
|
|
enddo
|
|
do j = 1, N_states
|
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
|
norm(j,ispin) = 0.d0
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0
|
|
else
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
|
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
|
endif
|
|
enddo
|
|
do i = 1, N_det
|
|
do j = 1, N_states
|
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
|
norm_bis(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_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
|
|
energies_alpha_beta(state_target, ispin) = - mo_bielec_integral_jj_exchange(orb_i,orb_v)
|
|
! energies_alpha_beta(state_target, ispin) = 0.d0
|
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
|
endif
|
|
enddo
|
|
enddo ! ispin
|
|
do state_target = 1, N_states
|
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
|
! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.5d0 * &
|
|
! ( energy_cas_dyall(state_target) - energies_alpha_beta(state_target,1) + &
|
|
! energy_cas_dyall(state_target) - energies_alpha_beta(state_target,2) )
|
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
|
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
|
else
|
|
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_States)]
|
|
implicit none
|
|
integer :: i,iorb,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: jorb,i_ok,aorb,orb_a
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
double precision :: hij
|
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
|
|
double precision :: energies_alpha_beta(N_states,2)
|
|
double precision :: norm_alpha_beta(N_states,2)
|
|
|
|
double precision :: thresh_norm
|
|
|
|
thresh_norm = 1.d-10
|
|
|
|
do aorb = 1,n_act_orb
|
|
orb_a = list_act(aorb)
|
|
do iorb = 1, n_inact_orb
|
|
orb_i = list_inact(iorb)
|
|
do state_target = 1, N_states
|
|
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
|
enddo
|
|
norm_alpha_beta = 0.d0
|
|
norm = 0.d0
|
|
norm_bis = 0.d0
|
|
do ispin = 1,2
|
|
do i = 1, n_det
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
|
enddo
|
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_a,ispin,i_ok)
|
|
if(i_ok.ne.1)then
|
|
do j = 1, n_states
|
|
psi_in_out_coef(i,j) = 0.d0
|
|
enddo
|
|
else
|
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
|
do j = 1, n_states
|
|
double precision :: coef,contrib
|
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
|
enddo
|
|
endif
|
|
enddo
|
|
do j = 1, N_states
|
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
|
norm(j,ispin) = 0.d0
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
else
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
|
endif
|
|
enddo
|
|
double precision :: norm_bis(N_states,2)
|
|
do i = 1, N_det
|
|
do j = 1, N_states
|
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
|
norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j)
|
|
enddo
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
|
|
psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
|
|
enddo
|
|
enddo
|
|
do state_target = 1, N_states
|
|
energies_alpha_beta(state_target, ispin) = 0.d0
|
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
|
endif
|
|
enddo
|
|
enddo ! ispin
|
|
do state_target = 1, N_states
|
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
|
one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall(state_target) - &
|
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
|
else
|
|
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
|
endif
|
|
! print*, '********'
|
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
|
! print*, one_anhil_inact(iorb,aorb,state_target)
|
|
! print*, one_creat(aorb,1,state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_States)]
|
|
implicit none
|
|
integer :: i,vorb,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i
|
|
integer :: orb_v
|
|
double precision :: norm_out(N_states_diag)
|
|
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,N_states_diag))
|
|
|
|
integer :: iorb,jorb,i_ok,aorb,orb_a
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
double precision :: hij
|
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
|
|
double precision :: energies_alpha_beta(N_states,2)
|
|
double precision :: norm_alpha_beta(N_states,2)
|
|
|
|
double precision :: thresh_norm
|
|
|
|
thresh_norm = 1.d-10
|
|
|
|
do aorb = 1,n_act_orb
|
|
orb_a = list_act(aorb)
|
|
do vorb = 1, n_virt_orb
|
|
orb_v = list_virt(vorb)
|
|
do state_target = 1, N_states
|
|
one_creat_virt(aorb,vorb,state_target) = 0.d0
|
|
enddo
|
|
norm_alpha_beta = 0.d0
|
|
norm = 0.d0
|
|
norm_bis = 0.d0
|
|
do ispin = 1,2
|
|
do i = 1, n_det
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
|
enddo
|
|
call do_mono_excitation(psi_in_out(1,1,i),orb_a,orb_v,ispin,i_ok)
|
|
if(i_ok.ne.1)then
|
|
do j = 1, n_states
|
|
psi_in_out_coef(i,j) = 0.d0
|
|
enddo
|
|
else
|
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
|
do j = 1, n_states
|
|
double precision :: coef,contrib
|
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
|
enddo
|
|
endif
|
|
enddo
|
|
do j = 1, N_states
|
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
|
norm(j,ispin) = 0.d0
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
else
|
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
|
endif
|
|
enddo
|
|
double precision :: norm_bis(N_states,2)
|
|
do i = 1, N_det
|
|
do j = 1, N_states
|
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
|
norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j)
|
|
enddo
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
|
|
psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
|
|
enddo
|
|
enddo
|
|
do state_target = 1, N_states
|
|
energies_alpha_beta(state_target, ispin) = 0.d0
|
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
|
! print*, energies(state_target)
|
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
|
endif
|
|
enddo
|
|
enddo ! ispin
|
|
do state_target = 1, N_states
|
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
|
one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
|
else
|
|
one_creat_virt(aorb,vorb,state_target) = 0.d0
|
|
endif
|
|
! print*, '********'
|
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
|
! print*, one_creat_virt(aorb,vorb,state_target)
|
|
! print*, one_anhil(aorb,1,state_target)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef)
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_bis, (n_inact_orb,n_virt_orb,N_det,N_States)]
|
|
&BEGIN_PROVIDER [ double precision, corr_e_from_1h1p, (N_States)]
|
|
implicit none
|
|
integer :: i,vorb,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i
|
|
integer :: orb_v
|
|
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
|
|
double precision :: delta_e_inact_virt(N_states)
|
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
|
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:)
|
|
use bitmasks
|
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1))
|
|
allocate (eigenvectors(size(H_matrix,1),N_det+1))
|
|
allocate (eigenvalues(N_det+1))
|
|
|
|
integer :: iorb,jorb,i_ok
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
double precision :: hij
|
|
double precision :: energies_alpha_beta(N_states,2)
|
|
|
|
|
|
double precision :: accu(N_states),norm
|
|
double precision :: amplitudes_alpha_beta(N_det,2)
|
|
double precision :: delta_e_alpha_beta(N_det,2)
|
|
|
|
corr_e_from_1h1p = 0.d0
|
|
do vorb = 1,n_virt_orb
|
|
orb_v = list_virt(vorb)
|
|
do iorb = 1, n_inact_orb
|
|
orb_i = list_inact(iorb)
|
|
! print*, '---------------------------------------------------------------------------'
|
|
do j = 1, N_states
|
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) &
|
|
- fock_virt_total_spin_trace(orb_v,j)
|
|
enddo
|
|
do ispin = 1,2
|
|
do i = 1, n_det
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
|
enddo
|
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
|
if(i_ok.ne.1)then
|
|
print*, orb_i,orb_v
|
|
call debug_det(psi_in_out,N_int)
|
|
print*, 'pb, i_ok ne 0 !!!'
|
|
endif
|
|
interact_psi0(i) = 0.d0
|
|
do j = 1 , N_det
|
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
|
|
interact_psi0(i) += hij * psi_coef(j,1)
|
|
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
|
|
call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij)
|
|
diag_elem(i) = hij
|
|
enddo
|
|
do state_target = 1, N_states
|
|
! Building the Hamiltonian matrix
|
|
H_matrix(1,1) = energy_cas_dyall(state_target)
|
|
do i = 1, N_det
|
|
! interaction with psi0
|
|
H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target)
|
|
H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target)
|
|
! diagonal elements
|
|
H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target)
|
|
! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1)
|
|
do j = i+1, N_det
|
|
call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij)
|
|
H_matrix(i+1,j+1) = hij !0.d0 !
|
|
H_matrix(j+1,i+1) = hij !0.d0 !
|
|
enddo
|
|
enddo
|
|
print*, '***'
|
|
do i = 1, N_det+1
|
|
write(*,'(100(F16.10,1X))')H_matrix(i,:)
|
|
enddo
|
|
call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1)
|
|
corr_e_from_1h1p(state_target) += eigenvalues(1) - energy_cas_dyall(state_target)
|
|
norm = 0.d0
|
|
do i = 1, N_det
|
|
psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1)
|
|
!! if(dabs(psi_coef(i,state_target)*) .gt. 1.d-8)then
|
|
if(dabs(psi_in_out_coef(i,state_target)) .gt. 1.d-8)then
|
|
! if(dabs(interact_psi0(i)) .gt. 1.d-8)then
|
|
delta_e_alpha_beta(i,ispin) = H_matrix(1,i+1) / psi_in_out_coef(i,state_target)
|
|
! delta_e_alpha_beta(i,ispin) = interact_psi0(i) / psi_in_out_coef(i,state_target)
|
|
amplitudes_alpha_beta(i,ispin) = psi_in_out_coef(i,state_target) / psi_coef(i,state_target)
|
|
else
|
|
amplitudes_alpha_beta(i,ispin) = 0.d0
|
|
delta_e_alpha_beta(i,ispin) = delta_e_inact_virt(state_target)
|
|
endif
|
|
!! one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,ispin,state_target) = amplitudes_alpha_beta(i,ispin)
|
|
norm += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target)
|
|
enddo
|
|
print*, 'Coef '
|
|
write(*,'(100(1X,F16.10))')psi_coef(1:N_det,state_target)
|
|
write(*,'(100(1X,F16.10))')psi_in_out_coef(:,state_target)
|
|
double precision :: coef_tmp(N_det)
|
|
do i = 1, N_det
|
|
coef_tmp(i) = psi_coef(i,1) * interact_psi0(i) / delta_e_alpha_beta(i,ispin)
|
|
enddo
|
|
write(*,'(100(1X,F16.10))')coef_tmp(:)
|
|
print*, 'naked interactions'
|
|
write(*,'(100(1X,F16.10))')interact_psi0(:)
|
|
print*, ''
|
|
|
|
print*, 'norm ',norm
|
|
norm = 1.d0/(norm)
|
|
accu(state_target) = 0.d0
|
|
do i = 1, N_det
|
|
accu(state_target) += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target) * H_matrix(i+1,i+1)
|
|
do j = i+1, N_det
|
|
accu(state_target) += 2.d0 * psi_in_out_coef(i,state_target) * psi_in_out_coef(j,state_target) * H_matrix(i+1,j+1)
|
|
enddo
|
|
enddo
|
|
accu(state_target) = accu(state_target) * norm
|
|
print*, delta_e_inact_virt(state_target)
|
|
print*, eigenvalues(1),accu(state_target),eigenvectors(1,1)
|
|
print*, energy_cas_dyall(state_target) - accu(state_target), one_anhil_one_creat_inact_virt(iorb,vorb,state_target) + delta_e_inact_virt(state_target)
|
|
|
|
enddo
|
|
enddo ! ispin
|
|
do state_target = 1, N_states
|
|
do i = 1, N_det
|
|
one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,state_target) = 0.5d0 * &
|
|
( delta_e_alpha_beta(i,1) + delta_e_alpha_beta(i,1))
|
|
enddo
|
|
enddo
|
|
print*, '***'
|
|
write(*,'(100(1X,F16.10))')
|
|
write(*,'(100(1X,F16.10))')delta_e_alpha_beta(:,2)
|
|
! write(*,'(100(1X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,1,:)
|
|
! write(*,'(100(1X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,2,:)
|
|
print*, '---------------------------------------------------------------------------'
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues)
|
|
print*, 'corr_e_from_1h1p,',corr_e_from_1h1p(:)
|
|
|
|
END_PROVIDER
|
|
|
|
subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles)
|
|
implicit none
|
|
double precision , intent(inout) :: matrix_1h1p(N_det,N_det,N_states)
|
|
double precision , intent(out) :: e_corr_from_1h1p_singles(N_states)
|
|
integer :: i,vorb,j
|
|
integer :: ispin,jspin
|
|
integer :: orb_i, hole_particle_i
|
|
integer :: orb_v
|
|
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
|
|
double precision :: delta_e_inact_virt(N_states)
|
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
|
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,N_states_diag),H_matrix(N_det+1,N_det+1))
|
|
allocate (eigenvectors(size(H_matrix,1),N_det+1))
|
|
allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det))
|
|
allocate (delta_e_det(N_det,N_det))
|
|
|
|
integer :: iorb,jorb,i_ok
|
|
integer :: state_target
|
|
double precision :: energies(n_states_diag)
|
|
double precision :: hij
|
|
double precision :: energies_alpha_beta(N_states,2)
|
|
double precision :: lamda_pt2(N_det)
|
|
|
|
|
|
double precision :: accu(N_states),norm
|
|
double precision :: amplitudes_alpha_beta(N_det,2)
|
|
double precision :: delta_e_alpha_beta(N_det,2)
|
|
double precision :: coef_array(N_states)
|
|
double precision :: coef_perturb(N_det)
|
|
double precision :: coef_perturb_bis(N_det)
|
|
|
|
do vorb = 1,n_virt_orb
|
|
orb_v = list_virt(vorb)
|
|
do iorb = 1, n_inact_orb
|
|
orb_i = list_inact(iorb)
|
|
do j = 1, N_states
|
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) &
|
|
- fock_virt_total_spin_trace(orb_v,j)
|
|
enddo
|
|
do ispin = 1,2
|
|
do i = 1, n_det
|
|
do j = 1, N_int
|
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
|
enddo
|
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
|
if(i_ok.ne.1)then
|
|
print*, orb_i,orb_v
|
|
call debug_det(psi_in_out,N_int)
|
|
print*, 'pb, i_ok ne 0 !!!'
|
|
endif
|
|
interact_psi0(i) = 0.d0
|
|
do j = 1 , N_det
|
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
|
|
call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j))
|
|
interact_cas(i,j) = hij
|
|
interact_psi0(i) += hij * psi_coef(j,1)
|
|
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
|
|
call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij)
|
|
diag_elem(i) = hij
|
|
enddo
|
|
do state_target = 1, N_states
|
|
! Building the Hamiltonian matrix
|
|
H_matrix(1,1) = energy_cas_dyall(state_target)
|
|
do i = 1, N_det
|
|
! interaction with psi0
|
|
H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target)
|
|
H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target)
|
|
! diagonal elements
|
|
H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target)
|
|
! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1)
|
|
do j = i+1, N_det
|
|
call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij)
|
|
H_matrix(i+1,j+1) = hij !0.d0 !
|
|
H_matrix(j+1,i+1) = hij !0.d0 !
|
|
enddo
|
|
enddo
|
|
call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1)
|
|
e_corr_from_1h1p_singles(state_target) += eigenvalues(1) - energy_cas_dyall(state_target)
|
|
|
|
do i = 1, N_det
|
|
psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1)
|
|
coef_perturb(i) = 0.d0
|
|
do j = 1, N_det
|
|
coef_perturb(i) += psi_coef(j,state_target) * interact_cas(i,j) *1.d0/delta_e_det(i,j)
|
|
enddo
|
|
coef_perturb_bis(i) = interact_psi0(i) / (eigenvalues(1) - H_matrix(i+1,i+1))
|
|
if(dabs(interact_psi0(i)) .gt. 1.d-12)then
|
|
lamda_pt2(i) = psi_in_out_coef(i,state_target) / interact_psi0(i)
|
|
else
|
|
lamda_pt2(i) =energy_cas_dyall(state_target) - H_matrix(i+1,i+1)
|
|
endif
|
|
enddo
|
|
if(dabs(eigenvalues(1) - energy_cas_dyall(state_target)).gt.1.d-10)then
|
|
print*, ''
|
|
do i = 1, N_det+1
|
|
write(*,'(100(F16.10))') H_matrix(i,:)
|
|
enddo
|
|
accu = 0.d0
|
|
do i = 1, N_det
|
|
accu(state_target) += psi_in_out_coef(i,state_target) * interact_psi0(i)
|
|
enddo
|
|
print*, ''
|
|
print*, 'e corr diagonal ',accu(state_target)
|
|
accu = 0.d0
|
|
do i = 1, N_det
|
|
accu(state_target) += coef_perturb(i) * interact_psi0(i)
|
|
enddo
|
|
print*, 'e corr perturb ',accu(state_target)
|
|
accu = 0.d0
|
|
do i = 1, N_det
|
|
accu(state_target) += coef_perturb_bis(i) * interact_psi0(i)
|
|
enddo
|
|
print*, 'e corr perturb EN',accu(state_target)
|
|
print*, ''
|
|
print*, 'coef diagonalized'
|
|
write(*,'(100(F16.10,1X))')psi_in_out_coef(:,state_target)
|
|
print*, 'coef_perturb'
|
|
write(*,'(100(F16.10,1X))')coef_perturb(:)
|
|
print*, 'coef_perturb EN'
|
|
write(*,'(100(F16.10,1X))')coef_perturb_bis(:)
|
|
endif
|
|
integer :: k
|
|
do k = 1, N_det
|
|
do i = 1, N_det
|
|
matrix_1h1p(i,i,state_target) += interact_cas(k,i) * interact_cas(k,i) * lamda_pt2(k)
|
|
do j = i+1, N_det
|
|
matrix_1h1p(i,j,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k)
|
|
matrix_1h1p(j,i,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo ! ispin
|
|
enddo
|
|
enddo
|
|
deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues,interact_cas)
|
|
deallocate(delta_e_det)
|
|
end
|