diff --git a/config/ifort.cfg b/config/ifort.cfg index da414912..a738a83c 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [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 OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 1eb2d45a..921b9a1a 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -7,7 +7,12 @@ s.set_selection_pt2("epstein_nesbet_2x2") s.unset_openmp() 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.unset_openmp() print s diff --git a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f index c65e89d8..c8140f70 100644 --- a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f +++ b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f @@ -10,7 +10,9 @@ end subroutine routine_3 implicit none - provide one_creation +!provide fock_virt_total_spin_trace + provide energy_cas_dyall + print*, 'nuclear_reuplsion = ',nuclear_repulsion end diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index dda6e11e..0c104be6 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -3,20 +3,21 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] 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) + 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) + print*, 'energy_cas_dyall(i)', energy_cas_dyall(i) enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)] +BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)] implicit none integer :: i,j integer :: ispin integer :: orb, hole_particle,spin_exc double precision :: norm_out(N_states_diag) - integer(bit_kind) :: psi_in_out(N_int,2,n_det) - double precision :: psi_in_out_coef(n_det,N_states_diag) + integer(bit_kind) :: psi_in_out(N_int,2,N_det) + double precision :: psi_in_out_coef(N_det,N_states_diag) use bitmasks 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo 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, & 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_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 END_PROVIDER -BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)] +BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)] implicit none integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo integer :: state_target @@ -76,14 +77,25 @@ BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)] double precision :: energies(n_states_diag) 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) +! 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) - 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 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo 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, & 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_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 @@ -130,7 +142,7 @@ BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)] 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo 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, & 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_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 @@ -177,7 +189,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)] 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(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 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) - 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 @@ -224,7 +236,7 @@ BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(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 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_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 @@ -282,7 +295,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) 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, & 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) - 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 @@ -340,7 +353,7 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo 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, & 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_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 @@ -398,7 +411,7 @@ BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_or 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 integer :: i,j 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) enddo do j = 1, N_int - psi_in_out(j,1,i) = psi_active(j,1,i) - psi_in_out(j,1,i) = psi_active(j,2,i) + psi_in_out(j,1,i) = psi_det(j,1,i) + psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo 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, & 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_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 diff --git a/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f b/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f index 6ecffcce..f6b542bd 100644 --- a/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f +++ b/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f @@ -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 integer :: i do i = 1, n_act_orb - one_creation_spin_averaged(i) = one_creation(i,1) + one_creation(i,2) - one_creation_spin_averaged(i) = 0.5d0 * one_creation_spin_averaged(i) + one_creat_spin_trace(i) = one_creat(i,1) + one_creat(i,2) + one_creat_spin_trace(i) = 0.5d0 * one_creat_spin_trace(i) enddo 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 integer :: i do i = 1, n_act_orb - one_anhilation_spin_averaged(i) = one_anhilation(i,1) + one_anhilation(i,2) - one_anhilation_spin_averaged(i) = 0.5d0 * one_anhilation_spin_averaged(i) + one_anhil_spin_trace(i) = one_anhil(i,1) + one_anhil(i,2) + one_anhil_spin_trace(i) = 0.5d0 * one_anhil_spin_trace(i) enddo 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 integer :: i,j integer :: ispin,jspin double precision :: counting do i = 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 do ispin = 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 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 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 integer :: i,j integer :: ispin,jspin double precision :: counting do i = 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 do ispin = 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 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 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 integer :: i,j integer :: ispin,jspin double precision :: counting do i = 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 do ispin = 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 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 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 integer :: i,j,k 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 j = 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 do ispin = 1, 2 do jspin = 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 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 @@ -108,7 +108,7 @@ BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n 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 integer :: i,j,k 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 j = 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 do ispin = 1, 2 do jspin = 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 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 @@ -135,7 +135,7 @@ BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (n 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 integer :: i,j,k 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 j = 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 do ispin = 1, 2 do jspin = 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 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 @@ -162,7 +162,7 @@ BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (n_act_orb,n_ac 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 integer :: i,j,k 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 j = 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 do ispin = 1, 2 do jspin = 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 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 diff --git a/plugins/MRPT_Utils/fock_like_operators.irp.f b/plugins/MRPT_Utils/fock_like_operators.irp.f index c7e21d0c..bb561455 100644 --- a/plugins/MRPT_Utils/fock_like_operators.irp.f +++ b/plugins/MRPT_Utils/fock_like_operators.irp.f @@ -40,7 +40,7 @@ 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) 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 END_PROVIDER @@ -49,11 +49,12 @@ ! inactive part of the fock operator with contributions only from the active END_DOC implicit none - integer :: i,j + integer :: i,j,k double precision :: accu_coulomb,accu_exchange(2) double precision :: na,nb,ntot 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 accu_coulomb = 0.d0 @@ -69,6 +70,17 @@ accu_coulomb += ntot * coulomb accu_exchange(1) += na * 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 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) @@ -80,11 +92,12 @@ ! virtual part of the fock operator with contributions only from the active END_DOC implicit none - integer :: i,j + integer :: i,j,k double precision :: accu_coulomb,accu_exchange(2) double precision :: na,nb,ntot 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 accu_coulomb = 0.d0 @@ -100,14 +113,26 @@ accu_coulomb += ntot * coulomb accu_exchange(1) += na * 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 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) + print*, fock_virt_from_act(i_virt_orb,1) , fock_virt_from_act(i_virt_orb,2) enddo END_PROVIDER 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 ! inactive part of the fock operator END_DOC @@ -118,12 +143,12 @@ 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,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 END_PROVIDER 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 ! inactive part of the fock operator END_DOC @@ -132,9 +157,10 @@ integer :: i_virt_orb do i = 1, n_virt_orb 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,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2) - 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(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)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb) + 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 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 ! active part of the fock operator with contributions only from the inactive END_DOC implicit none - integer :: i,j,k + integer :: i,j,k,k_inact_core_orb + integer :: iorb,jorb double precision :: accu double precision :: get_mo_bielec_integral,coulomb, exchange PROVIDE mo_bielec_integrals_in_map + fock_operator_active_from_core_inact = 0.d0 do i = 1, n_act_orb + iorb = list_act(i) do j = 1, n_act_orb + jorb = list_act(j) accu = 0.d0 do k = 1, n_core_inact_orb - coulomb = get_mo_bielec_integral(k,i,k,j,mo_integrals_map) - exchange = get_mo_bielec_integral(k,i,i,k,mo_integrals_map) + k_inact_core_orb = list_core_inact(k) + 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 enddo - fock_operator_active_from_core_inact(i,j) = accu + fock_operator_active_from_core_inact(iorb,jorb) = accu enddo enddo diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index e5848814..e6865bd8 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -7,6 +7,7 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] implicit none use bitmasks integer :: i,j,k,l + provide cas_bitmask do i = 1, N_det do j = 1, N_int 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) 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,2) = iand(virt_bitmask(i,2),xor(det_tmp_1(i,2),virt_bitmask(i,2))) + particles(i,1) = iand(virt_bitmask(i,1),det_tmp_1(i,1)) + particles(i,2) = iand(virt_bitmask(i,2),det_tmp_1(i,2)) enddo particles_list = 0 @@ -169,12 +170,12 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) delta_e_inactive = 0.d0 do i = 1, n_holes_spin(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 do i = 1, n_holes_spin(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 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 do i = 1, n_particles_spin(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 @@ -201,31 +207,49 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) 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, & 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 + icountbis = 0 do i = 1, n_holes_spin_act(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) enddo do i = 1, n_holes_spin_act(2) 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) 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 !!' stop endif icount = 0 + icountbis = 0 do i = 1, n_particles_spin_act(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) enddo do i = 1, n_particles_spin_act(2) 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) 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 !!' stop 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 - if (n_holes_act == 1 .and. n_particles_act == 0) then - i_hole_act = holes_active_list_spin_traced(1) - delta_e_act += one_creation_spin_averaged(i_hole_act) + integer :: ispin,jspin,kspin + if (n_holes_act == 0 .and. n_particles_act == 1) then +! 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 - i_particle_act = particles_active_list_spin_traced(1) - delta_e_act += one_anhilation_spin_averaged(i_particle_act) + else if (n_holes_act == 1 .and. n_particles_act == 0) then +! i_hole_act = holes_active_list_spin_traced(1) +! 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 - i_hole_act = holes_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) +! i_hole_act = holes_active_list_spin_traced(1) +! i_particle_act = particles_active_list_spin_traced(1) +! 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 - i_hole_act = holes_active_list_spin_traced(1) - j_hole_act = holes_active_list_spin_traced(2) - 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) +! i_hole_act = holes_active_list_spin_traced(1) +! j_hole_act = holes_active_list_spin_traced(2) +! i_particle_act = particles_active_list_spin_traced(1) +! 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 - i_hole_act = holes_active_list_spin_traced(1) - i_particle_act = particles_active_list_spin_traced(1) - 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) +! i_hole_act = holes_active_list_spin_traced(1) +! i_particle_act = particles_active_list_spin_traced(1) +! j_particle_act = particles_active_list_spin_traced(2) +! 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 - i_hole_act = holes_active_list_spin_traced(1) - j_hole_act = holes_active_list_spin_traced(2) - 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) +! i_hole_act = holes_active_list_spin_traced(1) +! j_hole_act = holes_active_list_spin_traced(2) +! k_hole_act = holes_active_list_spin_traced(3) +! 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 - i_particle_act = particles_active_list_spin_traced(1) - j_particle_act = particles_active_list_spin_traced(2) - 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) +! i_particle_act = particles_active_list_spin_traced(1) +! j_particle_act = particles_active_list_spin_traced(2) +! k_particle_act = particles_active_list_spin_traced(3) +! 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 - 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 diff --git a/plugins/Perturbation/NEEDED_CHILDREN_MODULES b/plugins/Perturbation/NEEDED_CHILDREN_MODULES index e29a6721..6a9bca47 100644 --- a/plugins/Perturbation/NEEDED_CHILDREN_MODULES +++ b/plugins/Perturbation/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Properties Hartree_Fock +Determinants Properties Hartree_Fock MRPT_Utils diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 0086c67e..2621a9c6 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -45,6 +45,36 @@ subroutine pt2_epstein_nesbet ($arguments) end + +subroutine pt2_decontracted ($arguments) + use bitmasks + implicit none + $declarations + + BEGIN_DOC + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem_fock, h + double precision :: i_H_psi_array(N_st) + double precision :: coef_pert + PROVIDE selection_criterion + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_pert_new_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array,coef_pert) + + + c_pert(1) = coef_pert + e_2_pert(1) = coef_pert * i_H_psi_array(1) +! print*,coef_pert,i_H_psi_array(1) + +end + + + + subroutine pt2_epstein_nesbet_2x2 ($arguments) use bitmasks implicit none @@ -68,8 +98,8 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) ASSERT (Nint > 0) 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_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(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) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st @@ -86,12 +116,29 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) c_pert(i) = 0.d0 endif 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 e_2_pert(i) = 0.d0 c_pert(i) = 0.d0 H_pert_diag(i) = 0.d0 endif enddo +! if( e_2_pert(1) .ne. 0.d0)then +! print*,' e_2_pert(1) ', e_2_pert(1) +! endif end diff --git a/plugins/Selectors_full/e_corr_selectors.irp.f b/plugins/Selectors_full/e_corr_selectors.irp.f index 952e1c23..fec480f0 100644 --- a/plugins/Selectors_full/e_corr_selectors.irp.f +++ b/plugins/Selectors_full/e_corr_selectors.irp.f @@ -56,7 +56,7 @@ END_PROVIDER i_H_HF_per_selectors(i) = hij E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij 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 coef_hf_selector = psi_selectors_coef(i,1) E_corr_per_selectors(i) = -1000.d0 diff --git a/plugins/Selectors_full/selectors.irp.f b/plugins/Selectors_full/selectors.irp.f index 826dcc4b..71c3550e 100644 --- a/plugins/Selectors_full/selectors.irp.f +++ b/plugins/Selectors_full/selectors.irp.f @@ -18,8 +18,10 @@ BEGIN_PROVIDER [ integer, N_det_selectors] N_det_selectors = N_det do i=1,N_det norm = norm + psi_average_norm_contrib_sorted(i) + if (norm > threshold_selectors) then - N_det_selectors = i-1 +! N_det_selectors = i-1 + N_det_selectors = i exit endif enddo diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index ec786941..735ca8e7 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1076,6 +1076,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max, end + subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions) use bitmasks implicit none diff --git a/src/Determinants/test_3d.irp.f b/src/Determinants/test_3d.irp.f index 748890da..a5d09cd3 100644 --- a/src/Determinants/test_3d.irp.f +++ b/src/Determinants/test_3d.irp.f @@ -2,14 +2,14 @@ program test_3d implicit none integer :: i,npt double precision :: dx,domain,x_min,x,step_function_becke - domain = 5.d0 - npt = 100 - dx = domain/dble(npt) - x_min = -0.5d0 * domain - x = x_min - do i = 1, npt - write(33,*)x,step_function_becke(x) - x += dx - enddo +!domain = 5.d0 +!npt = 100 +!dx = domain/dble(npt) +!x_min = -0.5d0 * domain +!x = x_min +!do i = 1, npt +! write(33,*)x,step_function_becke(x) +! x += dx +!enddo end