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

New perturbation theory is working

Browse Source
This commit is contained in:
Emmanuel Giner 2016-08-26 18:00:49 +02:00
parent 78254c68c9
commit 0075d01bd9
13 changed files with 374 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
# 0 : Deactivate # 0 : Deactivate
# #
[OPTION] [OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py CACHE : 1 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags OPENMP : 1 ; Append OpenMP flags

7
plugins/Full_CI/H_apply.irp.f
View File

@ -7,7 +7,12 @@ s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_openmp() s.unset_openmp()
print s print s
s = H_apply_zmq("FCI_PT2") s = H_apply("FCI_PT2_new")
s.set_perturbation("decontracted")
s.unset_openmp()
print s
s = H_apply("FCI_PT2")
s.set_perturbation("epstein_nesbet_2x2") s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp() s.unset_openmp()
print s print s

4
plugins/MRPT_Utils/MRPT_Utils.main.irp.f
View File

@ -10,7 +10,9 @@ end
subroutine routine_3 subroutine routine_3
implicit none implicit none
provide one_creation !provide fock_virt_total_spin_trace
provide energy_cas_dyall
print*, 'nuclear_reuplsion = ',nuclear_repulsion
end end

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

@ -3,20 +3,21 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
integer :: i integer :: i
double precision :: energies(N_states_diag) double precision :: energies(N_states_diag)
do i = 1, N_states 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) call u0_H_dyall_u0(energies,psi_det,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
energy_cas_dyall(i) = energies(i) energy_cas_dyall(i) = energies(i)
print*, 'energy_cas_dyall(i)', energy_cas_dyall(i)
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)] BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin integer :: ispin
integer :: orb, hole_particle,spin_exc integer :: orb, hole_particle,spin_exc
double precision :: norm_out(N_states_diag) double precision :: norm_out(N_states_diag)
integer(bit_kind) :: psi_in_out(N_int,2,n_det) integer(bit_kind) :: psi_in_out(N_int,2,N_det)
double precision :: psi_in_out_coef(n_det,N_states_diag) double precision :: psi_in_out_coef(N_det,N_states_diag)
use bitmasks use bitmasks
integer :: iorb integer :: iorb
@ -30,8 +31,8 @@ BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)]
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
integer :: state_target integer :: state_target
@ -40,13 +41,13 @@ BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)]
call apply_exc_to_psi(orb,hole_particle,spin_exc, & 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) 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) 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_creation(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target) one_creat(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)] BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin integer :: ispin
@ -67,8 +68,8 @@ BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)]
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
integer :: state_target integer :: state_target
@ -76,14 +77,25 @@ BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)]
double precision :: energies(n_states_diag) double precision :: energies(n_states_diag)
call apply_exc_to_psi(orb,hole_particle,spin_exc, & 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) norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
! do j = 1, n_det
! print*, 'psi_in_out_coef'
! print*, psi_in_out_coef(j,1)
! call debug_det(psi_in_out(1,1,j),N_int)
! enddo
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) 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_anhilation(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target) ! print*,'energy_cas_dyall(state_target)'
! print*, energy_cas_dyall(state_target)
! print*,'energies(state_target)'
! print*, energies(state_target)
one_anhil(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target)
! print*,'one_anhil(iorb,ispin)'
! print*, one_anhil(iorb,ispin)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)] BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
@ -113,8 +125,8 @@ BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)]
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
@ -122,7 +134,7 @@ BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)]
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 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) 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) 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_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) two_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -130,7 +142,7 @@ BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)] BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
@ -160,8 +172,8 @@ BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)]
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
@ -169,7 +181,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)]
call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & 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) 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) 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_anhilation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) two_anhil(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -177,7 +189,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act_orb,2,2)] BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
@ -207,16 +219,16 @@ BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
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, & 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) 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) 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_anhilation_one_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) one_anhil_one_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -224,7 +236,7 @@ BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -261,18 +273,19 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
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, & 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) 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, & 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) 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) 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_anhilation_one_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -282,7 +295,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -319,18 +332,18 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
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, & 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) 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, & 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) 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) 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_creation_one_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -340,7 +353,7 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -377,8 +390,8 @@ BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_or
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
@ -388,7 +401,7 @@ BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_or
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 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) 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) 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_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) three_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo
@ -398,7 +411,7 @@ BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_or
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, three_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -435,8 +448,8 @@ BEGIN_PROVIDER [ double precision, three_anhilation, (n_act_orb,n_act_orb,n_act_
psi_in_out_coef(i,j) = psi_coef(i,j) psi_in_out_coef(i,j) = psi_coef(i,j)
enddo enddo
do j = 1, N_int do j = 1, N_int
psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,1,i) = psi_det(j,1,i)
psi_in_out(j,1,i) = psi_active(j,2,i) psi_in_out(j,2,i) = psi_det(j,2,i)
enddo enddo
enddo enddo
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
@ -446,7 +459,7 @@ BEGIN_PROVIDER [ double precision, three_anhilation, (n_act_orb,n_act_orb,n_act_
call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & 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) 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) 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_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) three_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target)
enddo enddo
enddo enddo
enddo enddo

68
plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f
View File

@ -1,88 +1,88 @@
BEGIN_PROVIDER [ double precision, one_creation_spin_averaged, (n_act_orb)] BEGIN_PROVIDER [ double precision, one_creat_spin_trace, (n_act_orb)]
implicit none implicit none
integer :: i integer :: i
do i = 1, n_act_orb do i = 1, n_act_orb
one_creation_spin_averaged(i) = one_creation(i,1) + one_creation(i,2) one_creat_spin_trace(i) = one_creat(i,1) + one_creat(i,2)
one_creation_spin_averaged(i) = 0.5d0 * one_creation_spin_averaged(i) one_creat_spin_trace(i) = 0.5d0 * one_creat_spin_trace(i)
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_anhilation_spin_averaged, (n_act_orb)] BEGIN_PROVIDER [ double precision, one_anhil_spin_trace, (n_act_orb)]
implicit none implicit none
integer :: i integer :: i
do i = 1, n_act_orb do i = 1, n_act_orb
one_anhilation_spin_averaged(i) = one_anhilation(i,1) + one_anhilation(i,2) one_anhil_spin_trace(i) = one_anhil(i,1) + one_anhil(i,2)
one_anhilation_spin_averaged(i) = 0.5d0 * one_anhilation_spin_averaged(i) one_anhil_spin_trace(i) = 0.5d0 * one_anhil_spin_trace(i)
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creation_spin_averaged, (n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, two_creat_spin_trace, (n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
double precision :: counting double precision :: counting
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
two_creation_spin_averaged(j,i) = 0.d0 two_creat_spin_trace(j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
two_creation_spin_averaged(j,i) += two_creation(j,i,ispin,jspin) two_creat_spin_trace(j,i) += two_creat(j,i,ispin,jspin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
two_creation_spin_averaged(j,i) = two_creation_spin_averaged(j,i) / counting two_creat_spin_trace(j,i) = two_creat_spin_trace(j,i) / counting
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_anhilation_spin_averaged, (n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, two_anhil_spin_trace, (n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
double precision :: counting double precision :: counting
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
two_anhilation_spin_averaged(j,i) = 0.d0 two_anhil_spin_trace(j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
two_anhilation_spin_averaged(j,i) += two_anhilation(j,i,ispin,jspin) two_anhil_spin_trace(j,i) += two_anhil(j,i,ispin,jspin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
two_anhilation_spin_averaged(j,i) = two_anhilation_spin_averaged(j,i) / counting two_anhil_spin_trace(j,i) = two_anhil_spin_trace(j,i) / counting
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, one_anhil_one_creat_spin_trace, (n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j integer :: i,j
integer :: ispin,jspin integer :: ispin,jspin
double precision :: counting double precision :: counting
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
one_anhilation_one_creation_spin_averaged(j,i) = 0.d0 one_anhil_one_creat_spin_trace(j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
one_anhilation_one_creation_spin_averaged(j,i) += one_anhilation_one_creation(j,i,jspin,ispin) one_anhil_one_creat_spin_trace(j,i) += one_anhil_one_creat(j,i,jspin,ispin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
one_anhilation_one_creation_spin_averaged(j,i) = one_anhilation_one_creation_spin_averaged(j,i) / counting one_anhil_one_creat_spin_trace(j,i) = one_anhil_one_creat_spin_trace(j,i) / counting
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, two_anhil_one_creat_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j,k integer :: i,j,k
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -91,16 +91,16 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
do k = 1, n_act_orb do k = 1, n_act_orb
two_anhilation_one_creation_spin_averaged(k,j,i) = 0.d0 two_anhil_one_creat_spin_trace(k,j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
do kspin = 1,2 do kspin = 1,2
two_anhilation_one_creation_spin_averaged(k,j,i) += two_anhilation_one_creation(k,j,i,kspin,jspin,ispin) two_anhil_one_creat_spin_trace(k,j,i) += two_anhil_one_creat(k,j,i,kspin,jspin,ispin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
two_anhilation_one_creation_spin_averaged(k,j,i) = two_anhilation_one_creation_spin_averaged(k,j,i) / counting two_anhil_one_creat_spin_trace(k,j,i) = two_anhil_one_creat_spin_trace(k,j,i) / counting
enddo enddo
enddo enddo
enddo enddo
@ -108,7 +108,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, two_creat_one_anhil_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j,k integer :: i,j,k
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -117,16 +117,16 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (n
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
do k = 1, n_act_orb do k = 1, n_act_orb
two_creation_one_anhilation_spin_averaged(k,j,i) = 0.d0 two_creat_one_anhil_spin_trace(k,j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
do kspin = 1,2 do kspin = 1,2
two_creation_one_anhilation_spin_averaged(k,j,i) += two_creation_one_anhilation(k,j,i,kspin,jspin,ispin) two_creat_one_anhil_spin_trace(k,j,i) += two_creat_one_anhil(k,j,i,kspin,jspin,ispin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
two_creation_one_anhilation_spin_averaged(k,j,i) = two_creation_one_anhilation_spin_averaged(k,j,i) / counting two_creat_one_anhil_spin_trace(k,j,i) = two_creat_one_anhil_spin_trace(k,j,i) / counting
enddo enddo
enddo enddo
enddo enddo
@ -135,7 +135,7 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (n
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, three_creat_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j,k integer :: i,j,k
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -144,16 +144,16 @@ BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (n_act_orb,n_ac
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
do k = 1, n_act_orb do k = 1, n_act_orb
three_creation_spin_averaged(k,j,i) = 0.d0 three_creat_spin_trace(k,j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
do kspin = 1,2 do kspin = 1,2
three_creation_spin_averaged(k,j,i) += three_creation(k,j,i,kspin,jspin,ispin) three_creat_spin_trace(k,j,i) += three_creat(k,j,i,kspin,jspin,ispin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
three_creation_spin_averaged(k,j,i) = three_creation_spin_averaged(k,j,i) / counting three_creat_spin_trace(k,j,i) = three_creat_spin_trace(k,j,i) / counting
enddo enddo
enddo enddo
enddo enddo
@ -162,7 +162,7 @@ BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (n_act_orb,n_ac
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, three_anhilation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)] BEGIN_PROVIDER [ double precision, three_anhil_spin_trace, (n_act_orb,n_act_orb,n_act_orb)]
implicit none implicit none
integer :: i,j,k integer :: i,j,k
integer :: ispin,jspin,kspin integer :: ispin,jspin,kspin
@ -171,16 +171,16 @@ BEGIN_PROVIDER [ double precision, three_anhilation_spin_averaged, (n_act_orb,n_
do i = 1, n_act_orb do i = 1, n_act_orb
do j = 1, n_act_orb do j = 1, n_act_orb
do k = 1, n_act_orb do k = 1, n_act_orb
three_anhilation_spin_averaged(k,j,i) = 0.d0 three_anhil_spin_trace(k,j,i) = 0.d0
counting = 0.d0 counting = 0.d0
do ispin = 1, 2 do ispin = 1, 2
do jspin = 1,2 do jspin = 1,2
do kspin = 1,2 do kspin = 1,2
three_anhilation_spin_averaged(k,j,i) += three_anhilation(k,j,i,kspin,jspin,ispin) three_anhil_spin_trace(k,j,i) += three_anhil(k,j,i,kspin,jspin,ispin)
counting += 1.d0 counting += 1.d0
enddo enddo
enddo enddo
three_anhilation_spin_averaged(k,j,i) = three_anhilation_spin_averaged(k,j,i) / counting three_anhil_spin_trace(k,j,i) = three_anhil_spin_trace(k,j,i) / counting
enddo enddo
enddo enddo
enddo enddo

63
plugins/MRPT_Utils/fock_like_operators.irp.f
View File

@ -40,7 +40,7 @@
accu += 2.d0 * mo_bielec_integral_jj(i_virt_orb,j_inact_core_orb) & accu += 2.d0 * mo_bielec_integral_jj(i_virt_orb,j_inact_core_orb) &
- mo_bielec_integral_jj_exchange(i_virt_orb,j_inact_core_orb) - mo_bielec_integral_jj_exchange(i_virt_orb,j_inact_core_orb)
enddo enddo
fock_virt_from_core_inact(i_virt_orb) = accu + mo_mono_elec_integral(i_virt_orb,i_virt_orb) fock_virt_from_core_inact(i_virt_orb) = accu
enddo enddo
END_PROVIDER END_PROVIDER
@ -49,11 +49,12 @@
! inactive part of the fock operator with contributions only from the active ! inactive part of the fock operator with contributions only from the active
END_DOC END_DOC
implicit none implicit none
integer :: i,j integer :: i,j,k
double precision :: accu_coulomb,accu_exchange(2) double precision :: accu_coulomb,accu_exchange(2)
double precision :: na,nb,ntot double precision :: na,nb,ntot
double precision :: coulomb, exchange double precision :: coulomb, exchange
integer :: j_act_orb,i_inact_core_orb double precision :: get_mo_bielec_integral_schwartz
integer :: j_act_orb,k_act_orb,i_inact_core_orb
do i = 1, n_core_inact_orb do i = 1, n_core_inact_orb
accu_coulomb = 0.d0 accu_coulomb = 0.d0
@ -69,6 +70,17 @@
accu_coulomb += ntot * coulomb accu_coulomb += ntot * coulomb
accu_exchange(1) += na * exchange accu_exchange(1) += na * exchange
accu_exchange(2) += nb * exchange accu_exchange(2) += nb * exchange
do k = j+1, n_act_orb
k_act_orb = list_act(k)
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb)
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb)
ntot = na + nb
coulomb = get_mo_bielec_integral_schwartz(j_act_orb,i_inact_core_orb,k_act_orb,i_inact_core_orb,mo_integrals_map)
exchange = get_mo_bielec_integral_schwartz(j_act_orb,k_act_orb,i_inact_core_orb,i_inact_core_orb,mo_integrals_map)
accu_coulomb += 2.d0 * ntot * coulomb
accu_exchange(1) += 2.d0 * na * exchange
accu_exchange(2) += 2.d0 * nb * exchange
enddo
enddo enddo
fock_core_inactive_from_act(i_inact_core_orb,1) = accu_coulomb + accu_exchange(1) fock_core_inactive_from_act(i_inact_core_orb,1) = accu_coulomb + accu_exchange(1)
fock_core_inactive_from_act(i_inact_core_orb,2) = accu_coulomb + accu_exchange(2) fock_core_inactive_from_act(i_inact_core_orb,2) = accu_coulomb + accu_exchange(2)
@ -80,11 +92,12 @@
! virtual part of the fock operator with contributions only from the active ! virtual part of the fock operator with contributions only from the active
END_DOC END_DOC
implicit none implicit none
integer :: i,j integer :: i,j,k
double precision :: accu_coulomb,accu_exchange(2) double precision :: accu_coulomb,accu_exchange(2)
double precision :: na,nb,ntot double precision :: na,nb,ntot
double precision :: coulomb, exchange double precision :: coulomb, exchange
integer :: j_act_orb,i_virt_orb double precision :: get_mo_bielec_integral_schwartz
integer :: j_act_orb,i_virt_orb,k_act_orb
do i = 1, n_virt_orb do i = 1, n_virt_orb
accu_coulomb = 0.d0 accu_coulomb = 0.d0
@ -100,14 +113,26 @@
accu_coulomb += ntot * coulomb accu_coulomb += ntot * coulomb
accu_exchange(1) += na * exchange accu_exchange(1) += na * exchange
accu_exchange(2) += nb * exchange accu_exchange(2) += nb * exchange
do k = j+1, n_act_orb
k_act_orb = list_act(k)
na = one_body_dm_mo_alpha(j_act_orb,k_act_orb)
nb = one_body_dm_mo_beta(j_act_orb,k_act_orb)
ntot = na + nb
coulomb = get_mo_bielec_integral_schwartz(j_act_orb,i_virt_orb,k_act_orb,i_virt_orb,mo_integrals_map)
exchange = get_mo_bielec_integral_schwartz(j_act_orb,k_act_orb,i_virt_orb,i_virt_orb,mo_integrals_map)
accu_coulomb += 2.d0 * ntot * coulomb
accu_exchange(1) += 2.d0 * na * exchange
accu_exchange(2) += 2.d0 * nb * exchange
enddo
enddo enddo
fock_virt_from_act(i_virt_orb,1) = accu_coulomb + accu_exchange(1) fock_virt_from_act(i_virt_orb,1) = accu_coulomb + accu_exchange(1)
fock_virt_from_act(i_virt_orb,2) = accu_coulomb + accu_exchange(2) fock_virt_from_act(i_virt_orb,2) = accu_coulomb + accu_exchange(2)
print*, fock_virt_from_act(i_virt_orb,1) , fock_virt_from_act(i_virt_orb,2)
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2)] BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2)]
&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_averaged, (mo_tot_num)] &BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_trace, (mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! inactive part of the fock operator ! inactive part of the fock operator
END_DOC END_DOC
@ -118,12 +143,12 @@
i_inact_core_orb = list_core_inact(i) i_inact_core_orb = list_core_inact(i)
fock_core_inactive_total(i_inact_core_orb,1) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1) fock_core_inactive_total(i_inact_core_orb,1) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1)
fock_core_inactive_total(i_inact_core_orb,2) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2) fock_core_inactive_total(i_inact_core_orb,2) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2)
fock_core_inactive_total_spin_averaged(i_inact_core_orb) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1) + fock_core_inactive_total(i_inact_core_orb,2)) fock_core_inactive_total_spin_trace(i_inact_core_orb) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1) + fock_core_inactive_total(i_inact_core_orb,2))
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)] BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)]
&BEGIN_PROVIDER [double precision, fock_virt_total_spin_averaged, (mo_tot_num)] &BEGIN_PROVIDER [double precision, fock_virt_total_spin_trace, (mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! inactive part of the fock operator ! inactive part of the fock operator
END_DOC END_DOC
@ -132,9 +157,10 @@
integer :: i_virt_orb integer :: i_virt_orb
do i = 1, n_virt_orb do i = 1, n_virt_orb
i_virt_orb= list_virt(i) i_virt_orb= list_virt(i)
fock_virt_total(i_virt_orb,1) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1) fock_virt_total(i_virt_orb,1) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
fock_virt_total(i_virt_orb,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2) fock_virt_total(i_virt_orb,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb)
fock_virt_total_spin_averaged(i_virt_orb) = 0.5d0 * ( fock_virt_total(i_virt_orb,1) + fock_virt_total(i_virt_orb,2) ) fock_virt_total_spin_trace(i_virt_orb) = 0.5d0 * ( fock_virt_total(i_virt_orb,1) + fock_virt_total(i_virt_orb,2) )
print*, fock_virt_total_spin_trace(i_virt_orb)
enddo enddo
END_PROVIDER END_PROVIDER
@ -142,24 +168,29 @@
BEGIN_PROVIDER [double precision, fock_operator_active_from_core_inact, (n_act_orb,n_act_orb)] BEGIN_PROVIDER [double precision, fock_operator_active_from_core_inact, (mo_tot_num,mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! active part of the fock operator with contributions only from the inactive ! active part of the fock operator with contributions only from the inactive
END_DOC END_DOC
implicit none implicit none
integer :: i,j,k integer :: i,j,k,k_inact_core_orb
integer :: iorb,jorb
double precision :: accu double precision :: accu
double precision :: get_mo_bielec_integral,coulomb, exchange double precision :: get_mo_bielec_integral,coulomb, exchange
PROVIDE mo_bielec_integrals_in_map PROVIDE mo_bielec_integrals_in_map
fock_operator_active_from_core_inact = 0.d0
do i = 1, n_act_orb do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_act_orb do j = 1, n_act_orb
jorb = list_act(j)
accu = 0.d0 accu = 0.d0
do k = 1, n_core_inact_orb do k = 1, n_core_inact_orb
coulomb = get_mo_bielec_integral(k,i,k,j,mo_integrals_map) k_inact_core_orb = list_core_inact(k)
exchange = get_mo_bielec_integral(k,i,i,k,mo_integrals_map) coulomb = get_mo_bielec_integral(k_inact_core_orb,iorb,k_inact_core_orb,jorb,mo_integrals_map)
exchange = get_mo_bielec_integral(k_inact_core_orb,jorb,iorb,k_inact_core_orb,mo_integrals_map)
accu += 2.d0 * coulomb - exchange accu += 2.d0 * coulomb - exchange
enddo enddo
fock_operator_active_from_core_inact(i,j) = accu fock_operator_active_from_core_inact(iorb,jorb) = accu
enddo enddo
enddo enddo

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

@ -7,6 +7,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
implicit none implicit none
use bitmasks use bitmasks
integer :: i,j,k,l integer :: i,j,k,l
provide cas_bitmask
do i = 1, N_det do i = 1, N_det
do j = 1, N_int do j = 1, N_int
psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1)) psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1))
@ -138,8 +139,8 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti
call give_virt_part_determinant(det_1,det_tmp_1) call give_virt_part_determinant(det_1,det_tmp_1)
do i = 1, N_int do i = 1, N_int
particles(i,1) = iand(virt_bitmask(i,1),xor(det_tmp_1(i,1),virt_bitmask(i,1))) particles(i,1) = iand(virt_bitmask(i,1),det_tmp_1(i,1))
particles(i,2) = iand(virt_bitmask(i,2),xor(det_tmp_1(i,2),virt_bitmask(i,2))) particles(i,2) = iand(virt_bitmask(i,2),det_tmp_1(i,2))
enddo enddo
particles_list = 0 particles_list = 0
@ -169,12 +170,12 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
delta_e_inactive = 0.d0 delta_e_inactive = 0.d0
do i = 1, n_holes_spin(1) do i = 1, n_holes_spin(1)
i_hole_inact = holes_list(i,1) i_hole_inact = holes_list(i,1)
delta_e_inactive += fock_core_inactive_total_spin_averaged(i_hole_inact) delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
enddo enddo
do i = 1, n_holes_spin(2) do i = 1, n_holes_spin(2)
i_hole_inact = holes_list(i,2) i_hole_inact = holes_list(i,2)
delta_e_inactive += fock_core_inactive_total_spin_averaged(i_hole_inact) delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact)
enddo enddo
double precision :: delta_e_virt double precision :: delta_e_virt
@ -187,7 +188,12 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
delta_e_virt = 0.d0 delta_e_virt = 0.d0
do i = 1, n_particles_spin(1) do i = 1, n_particles_spin(1)
i_part_virt = particles_list(i,1) i_part_virt = particles_list(i,1)
delta_e_virt += fock_virt_total_spin_averaged(i_part_virt) delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
enddo
do i = 1, n_particles_spin(2)
i_part_virt = particles_list(i,2)
delta_e_virt += fock_virt_total_spin_trace(i_part_virt)
enddo enddo
@ -201,31 +207,49 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
delta_e_act = 0.d0 delta_e_act = 0.d0
call give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin_act,n_particles_spin_act, & call give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin_act,n_particles_spin_act, &
n_holes_act,n_particles_act,holes_active_list,particles_active_list) n_holes_act,n_particles_act,holes_active_list,particles_active_list)
integer :: icount integer :: icount,icountbis
integer :: hole_list_practical(2,elec_num_tab(1)), particle_list_practical(2,elec_num_tab(1))
icount = 0 icount = 0
icountbis = 0
do i = 1, n_holes_spin_act(1) do i = 1, n_holes_spin_act(1)
icount += 1 icount += 1
icountbis += 1
hole_list_practical(1,icountbis) = 1
hole_list_practical(2,icountbis) = holes_active_list(i,1)
holes_active_list_spin_traced(icount) = holes_active_list(i,1) holes_active_list_spin_traced(icount) = holes_active_list(i,1)
enddo enddo
do i = 1, n_holes_spin_act(2) do i = 1, n_holes_spin_act(2)
icount += 1 icount += 1
icountbis += 1
hole_list_practical(1,icountbis) = 2
hole_list_practical(2,icountbis) = holes_active_list(i,2)
holes_active_list_spin_traced(icount) = holes_active_list(i,2) holes_active_list_spin_traced(icount) = holes_active_list(i,2)
enddo enddo
if(icount .ne. n_holes) then if(icount .ne. n_holes_act) then
print*,''
print*, icount, n_holes_act
print * , 'pb in holes_active_list_spin_traced !!' print * , 'pb in holes_active_list_spin_traced !!'
stop stop
endif endif
icount = 0 icount = 0
icountbis = 0
do i = 1, n_particles_spin_act(1) do i = 1, n_particles_spin_act(1)
icount += 1 icount += 1
icountbis += 1
particle_list_practical(1,icountbis) = 1
particle_list_practical(2,icountbis) = particles_active_list(i,1)
particles_active_list_spin_traced(icount) = particles_active_list(i,1) particles_active_list_spin_traced(icount) = particles_active_list(i,1)
enddo enddo
do i = 1, n_particles_spin_act(2) do i = 1, n_particles_spin_act(2)
icount += 1 icount += 1
icountbis += 1
particle_list_practical(1,icountbis) = 2
particle_list_practical(2,icountbis) = particles_active_list(i,2)
particles_active_list_spin_traced(icount) = particles_active_list(i,2) particles_active_list_spin_traced(icount) = particles_active_list(i,2)
enddo enddo
if(icount .ne. n_particles) then if(icount .ne. n_particles_act) then
print*, icount, n_particles_act
print * , 'pb in particles_active_list_spin_traced !!' print * , 'pb in particles_active_list_spin_traced !!'
stop stop
endif endif
@ -235,45 +259,148 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
integer :: i_particle_act, j_particle_act, k_particle_act integer :: i_particle_act, j_particle_act, k_particle_act
if (n_holes_act == 1 .and. n_particles_act == 0) then integer :: ispin,jspin,kspin
i_hole_act = holes_active_list_spin_traced(1) if (n_holes_act == 0 .and. n_particles_act == 1) then
delta_e_act += one_creation_spin_averaged(i_hole_act) ! i_particle_act = particles_active_list_spin_traced(1)
! delta_e_act += one_creat_spin_trace(i_particle_act )
ispin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
delta_e_act += one_creat(i_particle_act,ispin)
else if (n_holes_act == 0 .and. n_particles_act == 1) then else if (n_holes_act == 1 .and. n_particles_act == 0) then
i_particle_act = particles_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
delta_e_act += one_anhilation_spin_averaged(i_particle_act) ! delta_e_act += one_anhil_spin_trace(i_hole_act )
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
delta_e_act += one_anhil(i_hole_act , ispin)
else if (n_holes_act == 1 .and. n_particles_act == 1) then else if (n_holes_act == 1 .and. n_particles_act == 1) then
i_hole_act = holes_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
i_particle_act = particles_active_list_spin_traced(1) ! i_particle_act = particles_active_list_spin_traced(1)
delta_e_act += one_anhilation_one_creation_spin_averaged(i_hole_act,i_particle_act) ! delta_e_act += one_anhil_one_creat_spin_trace(i_hole_act,i_particle_act)
! first hole
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
! first particle
jspin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
delta_e_act += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin)
else if (n_holes_act == 2 .and. n_particles_act == 0) then
! i_hole_act = holes_active_list_spin_traced(1)
! j_hole_act = holes_active_list_spin_traced(1)
! delta_e_act += two_anhil_spin_trace(i_hole_act,j_hole_act)
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
jspin = hole_list_practical(1,2)
j_hole_act = hole_list_practical(2,2)
delta_e_act += two_anhil(i_hole_act,j_hole_act,ispin,jspin)
else if (n_holes_act == 0 .and. n_particles_act == 2) then
! i_particle_act = particles_active_list_spin_traced(1)
! j_particle_act = particles_active_list_spin_traced(2)
! delta_e_act += two_creat_spin_trace(i_particle_act,j_particle_act)
ispin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
jspin = particle_list_practical(1,2)
j_particle_act = particle_list_practical(2,2)
delta_e_act += two_creat(i_particle_act,j_particle_act,ispin,jspin)
else if (n_holes_act == 2 .and. n_particles_act == 1) then else if (n_holes_act == 2 .and. n_particles_act == 1) then
i_hole_act = holes_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
j_hole_act = holes_active_list_spin_traced(2) ! j_hole_act = holes_active_list_spin_traced(2)
i_particle_act = particles_active_list_spin_traced(1) ! i_particle_act = particles_active_list_spin_traced(1)
delta_e_act += two_anhilation_one_creation_spin_averaged(i_hole_act,j_hole_act,i_particle_act) ! print*, 'i_hole_act,j_hole_act,i_particle_act'
! print*, i_hole_act,j_hole_act,i_particle_act
! print*, two_anhil_one_creat_spin_trace(i_hole_act,j_hole_act,i_particle_act)
! delta_e_act += two_anhil_one_creat_spin_trace(i_hole_act,j_hole_act,i_particle_act)
! first hole
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
! second hole
jspin = hole_list_practical(1,2)
j_hole_act = hole_list_practical(2,2)
! first particle
kspin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
delta_e_act += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin)
else if (n_holes_act == 1 .and. n_particles_act == 2) then else if (n_holes_act == 1 .and. n_particles_act == 2) then
i_hole_act = holes_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
i_particle_act = particles_active_list_spin_traced(1) ! i_particle_act = particles_active_list_spin_traced(1)
j_particle_act = particles_active_list_spin_traced(2) ! j_particle_act = particles_active_list_spin_traced(2)
delta_e_act += two_creation_one_anhilation_spin_averaged(i_hole_act,i_particle_act,j_particle_act) ! delta_e_act += two_creat_one_anhil_spin_trace(i_hole_act,i_particle_act,j_particle_act)
! first hole
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
! 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)
delta_e_act += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin)
else if (n_holes_act == 3 .and. n_particles_act == 0) then else if (n_holes_act == 3 .and. n_particles_act == 0) then
i_hole_act = holes_active_list_spin_traced(1) ! i_hole_act = holes_active_list_spin_traced(1)
j_hole_act = holes_active_list_spin_traced(2) ! j_hole_act = holes_active_list_spin_traced(2)
k_hole_act = holes_active_list_spin_traced(3) ! k_hole_act = holes_active_list_spin_traced(3)
delta_e_act += three_anhilation_spin_averaged(i_hole_act,j_hole_act,k_hole_act) ! delta_e_act += three_anhil_spin_trace(i_hole_act,j_hole_act,k_hole_act)
! first hole
ispin = hole_list_practical(1,1)
i_hole_act = hole_list_practical(2,1)
! second hole
jspin = hole_list_practical(1,2)
j_hole_act = hole_list_practical(2,2)
! third hole
kspin = hole_list_practical(1,3)
k_hole_act = hole_list_practical(2,3)
delta_e_act += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin)
else if (n_holes_act == 0 .and. n_particles_act == 3) then else if (n_holes_act == 0 .and. n_particles_act == 3) then
i_particle_act = particles_active_list_spin_traced(1) ! i_particle_act = particles_active_list_spin_traced(1)
j_particle_act = particles_active_list_spin_traced(2) ! j_particle_act = particles_active_list_spin_traced(2)
k_particle_act = particles_active_list_spin_traced(3) ! k_particle_act = particles_active_list_spin_traced(3)
delta_e_act += three_creation_spin_averaged(i_particle_act,j_particle_act,k_particle_act) ! delta_e_act += three_creat_spin_trace(i_particle_act,j_particle_act,k_particle_act)
! first particle
ispin = particle_list_practical(1,1)
i_particle_act = particle_list_practical(2,1)
! second particle
jspin = particle_list_practical(1,2)
j_particle_act = particle_list_practical(2,2)
! second particle
kspin = particle_list_practical(1,3)
k_particle_act = particle_list_practical(2,3)
delta_e_act += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin)
endif endif
delta_e_final = delta_e_act + delta_e_virt + delta_e_inactive !print*, 'one_anhil_spin_trace'
!print*, one_anhil_spin_trace(1), one_anhil_spin_trace(2)
delta_e_final = delta_e_act + delta_e_inactive - delta_e_virt
if(delta_e_final .le. -100d0.or.delta_e_final > 0.d0)then
call debug_det(det_1,N_int)
call debug_det(det_2,N_int)
print*, 'n_holes_act,n_particles_act'
print*, n_holes_act,n_particles_act
print*, 'delta_e_act,delta_e_inactive,delta_e_vir'
print*, delta_e_act,delta_e_inactive,delta_e_virt
stop
endif
!if(delta_e_final > 0.d0)then
!print*, delta_e_final
!stop
!endif
end end

2
plugins/Perturbation/NEEDED_CHILDREN_MODULES
View File

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

51
plugins/Perturbation/pt2_equations.irp.f
View File

@ -45,6 +45,36 @@ subroutine pt2_epstein_nesbet ($arguments)
end 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)
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) subroutine pt2_epstein_nesbet_2x2 ($arguments)
use bitmasks use bitmasks
implicit none implicit none
@ -68,8 +98,8 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
ASSERT (Nint > 0) ASSERT (Nint > 0)
PROVIDE CI_electronic_energy PROVIDE CI_electronic_energy
!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(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) !call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st do i =1,N_st
@ -86,12 +116,29 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
c_pert(i) = 0.d0 c_pert(i) = 0.d0
endif endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i) H_pert_diag(i) = h*c_pert(i)*c_pert(i)
! print*, 'N_det,N_det_selectors = ',N_det,N_det_selectors
! print*, 'threshold_selectors',threshold_selectors
! print*, delta_e,i_H_psi_array(1)
! double precision :: hij,accu
! accu = 0.d0
! do j = 1, N_det
! call i_H_j(det_pert,psi_selectors(1,1,j),N_int,hij)
! print*, 'psi_selectors_coef(j,1 = ',psi_selectors_coef(j,1),psi_coef(j,1)
! call debug_det(psi_det(1,1,i),N_int)
! call debug_det(psi_selectors(1,1,i),N_int)
! accu += psi_selectors_coef(j,1) * hij
! enddo
! print*, 'accu,ihpsi0',accu,i_H_psi_array(1)
! stop
else else
e_2_pert(i) = 0.d0 e_2_pert(i) = 0.d0
c_pert(i) = 0.d0 c_pert(i) = 0.d0
H_pert_diag(i) = 0.d0 H_pert_diag(i) = 0.d0
endif endif
enddo enddo
! if( e_2_pert(1) .ne. 0.d0)then
! print*,' e_2_pert(1) ', e_2_pert(1)
! endif
end end

2
plugins/Selectors_full/e_corr_selectors.irp.f
View File

@ -56,7 +56,7 @@ END_PROVIDER
i_H_HF_per_selectors(i) = hij i_H_HF_per_selectors(i) = hij
E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij
E_corr_double_only += E_corr_per_selectors(i) E_corr_double_only += E_corr_per_selectors(i)
E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int)) ! E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int))
elseif(exc_degree_per_selectors(i) == 0)then elseif(exc_degree_per_selectors(i) == 0)then
coef_hf_selector = psi_selectors_coef(i,1) coef_hf_selector = psi_selectors_coef(i,1)
E_corr_per_selectors(i) = -1000.d0 E_corr_per_selectors(i) = -1000.d0

4
plugins/Selectors_full/selectors.irp.f
View File

@ -18,8 +18,10 @@ BEGIN_PROVIDER [ integer, N_det_selectors]
N_det_selectors = N_det N_det_selectors = N_det
do i=1,N_det do i=1,N_det
norm = norm + psi_average_norm_contrib_sorted(i) norm = norm + psi_average_norm_contrib_sorted(i)
if (norm > threshold_selectors) then if (norm > threshold_selectors) then
N_det_selectors = i-1 ! N_det_selectors = i-1
N_det_selectors = i
exit exit
endif endif
enddo enddo

1
src/Determinants/slater_rules.irp.f
View File

@ -1076,6 +1076,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,
end end
subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions) subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions)
use bitmasks use bitmasks
implicit none implicit none

18
src/Determinants/test_3d.irp.f
View File

@ -2,14 +2,14 @@ program test_3d
implicit none implicit none
integer :: i,npt integer :: i,npt
double precision :: dx,domain,x_min,x,step_function_becke double precision :: dx,domain,x_min,x,step_function_becke
domain = 5.d0 !domain = 5.d0
npt = 100 !npt = 100
dx = domain/dble(npt) !dx = domain/dble(npt)
x_min = -0.5d0 * domain !x_min = -0.5d0 * domain
x = x_min !x = x_min
do i = 1, npt !do i = 1, npt
write(33,*)x,step_function_becke(x) ! write(33,*)x,step_function_becke(x)
x += dx ! x += dx
enddo !enddo
end end