LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity
b2a5d497a0
quantum_package/plugins/MRPT_Utils/energies_cas.irp.f
Anthony Scemama ca973a1e92 Revert "Bugs to fix (#50)" (#51) 
This reverts commit 94f01c0892.
2017-04-20 08:45:56 +02:00

1115 lines
40 KiB
Fortran
Raw Blame History

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
Reference in New Issue View Git Blame Copy Permalink