From dbf894a99a71c4bddf4f21ee5effd01b7bd71773 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 1 Sep 2016 17:43:33 +0200 Subject: [PATCH] mrpt new with multi state version --- config/ifort.cfg | 2 +- plugins/MRPT_Utils/MRPT_Utils.main.irp.f | 32 +- plugins/MRPT_Utils/energies_cas.irp.f | 167 +++++---- .../energies_cas_spin_averaged.irp.f | 190 ---------- plugins/MRPT_Utils/fock_like_operators.irp.f | 145 ++++---- plugins/MRPT_Utils/mrpt_dress.irp.f | 17 +- plugins/MRPT_Utils/mrpt_utils.irp.f | 79 ++-- plugins/MRPT_Utils/psi_active_prov.irp.f | 337 ++++-------------- src/Bitmask/bitmask_cas_routines.irp.f | 6 - src/Determinants/density_matrix.irp.f | 61 ++-- 10 files changed, 319 insertions(+), 717 deletions(-) delete mode 100644 plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f 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/MRPT_Utils/MRPT_Utils.main.irp.f b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f index 7e4aa80b..9ee42820 100644 --- a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f +++ b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f @@ -27,43 +27,17 @@ subroutine routine_2 implicit none integer :: i do i = 1, n_core_inact_orb - print*,fock_core_inactive_total(i,1),fock_core_inactive(i) + print*,fock_core_inactive_total(i,1,1),fock_core_inactive(i) enddo double precision :: accu accu = 0.d0 do i = 1, n_act_orb integer :: j_act_orb j_act_orb = list_act(i) - accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb) - print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb),one_body_dm_mo_beta(j_act_orb,j_act_orb) + accu += one_body_dm_mo_alpha(j_act_orb,j_act_orb,1) + print*,one_body_dm_mo_alpha(j_act_orb,j_act_orb,1),one_body_dm_mo_beta(j_act_orb,j_act_orb,1) enddo print*,'accu = ',accu end -subroutine routine - implicit none - integer :: i,j - integer :: orb, spin_exc - integer :: hole_particle - double precision, allocatable :: norm_out(:) - allocate(norm_out(N_states_diag)) - - orb = list_virt(10) - hole_particle = -1 - spin_exc = 1 - - call apply_exc_to_psi(orb,hole_particle,spin_exc, & - norm_out,psi_det,psi_coef, n_det,psi_det_size,psi_det_size,N_states_diag) - do i = 1, N_det - if(psi_coef(i,1).ne.0.d0)then - print*, '' - call debug_det(psi_det(1,1,i),N_int) - print*, 'coef = ',psi_coef(i,1) - endif - enddo - print*,'norm_out = ',norm_out - - deallocate(norm_out) - -end diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index 0c104be6..8644bfa8 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -10,7 +10,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] END_PROVIDER -BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)] +BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] implicit none integer :: i,j integer :: ispin @@ -21,6 +21,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)] use bitmasks integer :: iorb + integer :: state_target + double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) @@ -35,19 +37,18 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2)] psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo - integer :: state_target - state_target = 1 - 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) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - one_creat(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1,N_states + call apply_exc_to_psi(orb,hole_particle,spin_exc, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)] +BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] implicit none integer :: i,j integer :: ispin @@ -58,6 +59,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)] use bitmasks integer :: iorb + integer :: state_target + double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) @@ -72,30 +75,18 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2)] psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo - integer :: state_target - state_target = 1 - 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) -! 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) + do state_target = 1, N_states + call apply_exc_to_psi(orb,hole_particle,spin_exc, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)] +BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin @@ -108,7 +99,6 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)] integer :: iorb,jorb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -129,12 +119,14 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)] 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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - two_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1 , N_states + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -142,7 +134,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2)] END_PROVIDER -BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)] +BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin @@ -181,7 +173,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (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_anhil(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo @@ -189,7 +181,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2)] END_PROVIDER -BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2)] +BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2,N_States)] implicit none integer :: i,j integer :: ispin,jspin @@ -202,7 +194,6 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 integer :: iorb,jorb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -223,12 +214,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 psi_in_out(j,2,i) = psi_det(j,2,i) enddo enddo - 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_anhil_one_creat(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1, N_states + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -236,7 +229,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 END_PROVIDER -BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (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,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -251,7 +244,6 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a integer :: iorb,jorb integer :: korb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -278,14 +270,16 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a enddo enddo - call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) - call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) - call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1, N_states + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -295,7 +289,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a END_PROVIDER -BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (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,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -310,7 +304,6 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer :: iorb,jorb integer :: korb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -336,14 +329,16 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a 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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1, N_states + call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -353,7 +348,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a END_PROVIDER -BEGIN_PROVIDER [ double precision, three_creat, (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,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -368,7 +363,6 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 integer :: iorb,jorb integer :: korb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -394,14 +388,16 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - three_creat(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1, N_states + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo @@ -411,7 +407,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 END_PROVIDER -BEGIN_PROVIDER [ double precision, three_anhil, (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,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin @@ -426,7 +422,6 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 integer :: iorb,jorb integer :: korb integer :: state_target - state_target = 1 double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 @@ -452,14 +447,16 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) - three_anhil(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + do state_target = 1, N_states + call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + enddo enddo enddo enddo diff --git a/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f b/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f deleted file mode 100644 index f6b542bd..00000000 --- a/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f +++ /dev/null @@ -1,190 +0,0 @@ - -BEGIN_PROVIDER [ double precision, one_creat_spin_trace, (n_act_orb)] - implicit none - integer :: i - do i = 1, n_act_orb - 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_anhil_spin_trace, (n_act_orb)] - implicit none - integer :: i - do i = 1, n_act_orb - 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_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_creat_spin_trace(j,i) = 0.d0 - counting = 0.d0 - do ispin = 1, 2 - do jspin = 1,2 - two_creat_spin_trace(j,i) += two_creat(j,i,ispin,jspin) - counting += 1.d0 - enddo - enddo - two_creat_spin_trace(j,i) = two_creat_spin_trace(j,i) / counting - enddo - enddo -END_PROVIDER - -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_anhil_spin_trace(j,i) = 0.d0 - counting = 0.d0 - do ispin = 1, 2 - do jspin = 1,2 - two_anhil_spin_trace(j,i) += two_anhil(j,i,ispin,jspin) - counting += 1.d0 - enddo - enddo - two_anhil_spin_trace(j,i) = two_anhil_spin_trace(j,i) / counting - enddo - enddo -END_PROVIDER - - -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_anhil_one_creat_spin_trace(j,i) = 0.d0 - counting = 0.d0 - do ispin = 1, 2 - do jspin = 1,2 - one_anhil_one_creat_spin_trace(j,i) += one_anhil_one_creat(j,i,jspin,ispin) - counting += 1.d0 - enddo - enddo - 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_anhil_one_creat_spin_trace, (n_act_orb,n_act_orb,n_act_orb)] - implicit none - integer :: i,j,k - integer :: ispin,jspin,kspin - double precision :: counting - - do i = 1, n_act_orb - do j = 1, n_act_orb - do k = 1, n_act_orb - 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_anhil_one_creat_spin_trace(k,j,i) += two_anhil_one_creat(k,j,i,kspin,jspin,ispin) - counting += 1.d0 - enddo - enddo - two_anhil_one_creat_spin_trace(k,j,i) = two_anhil_one_creat_spin_trace(k,j,i) / counting - enddo - enddo - enddo - enddo - -END_PROVIDER - -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 - double precision :: counting - - do i = 1, n_act_orb - do j = 1, n_act_orb - do k = 1, n_act_orb - 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_creat_one_anhil_spin_trace(k,j,i) += two_creat_one_anhil(k,j,i,kspin,jspin,ispin) - counting += 1.d0 - enddo - enddo - two_creat_one_anhil_spin_trace(k,j,i) = two_creat_one_anhil_spin_trace(k,j,i) / counting - enddo - enddo - enddo - enddo - -END_PROVIDER - - -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 - double precision :: counting - - do i = 1, n_act_orb - do j = 1, n_act_orb - do k = 1, n_act_orb - 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_creat_spin_trace(k,j,i) += three_creat(k,j,i,kspin,jspin,ispin) - counting += 1.d0 - enddo - enddo - three_creat_spin_trace(k,j,i) = three_creat_spin_trace(k,j,i) / counting - enddo - enddo - enddo - enddo - -END_PROVIDER - - -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 - double precision :: counting - - do i = 1, n_act_orb - do j = 1, n_act_orb - do k = 1, n_act_orb - 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_anhil_spin_trace(k,j,i) += three_anhil(k,j,i,kspin,jspin,ispin) - counting += 1.d0 - enddo - enddo - three_anhil_spin_trace(k,j,i) = three_anhil_spin_trace(k,j,i) / counting - enddo - enddo - enddo - enddo - -END_PROVIDER - diff --git a/plugins/MRPT_Utils/fock_like_operators.irp.f b/plugins/MRPT_Utils/fock_like_operators.irp.f index 44c16da6..5900516e 100644 --- a/plugins/MRPT_Utils/fock_like_operators.irp.f +++ b/plugins/MRPT_Utils/fock_like_operators.irp.f @@ -44,7 +44,7 @@ enddo END_PROVIDER - BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2)] + BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2,N_states)] BEGIN_DOC ! inactive part of the fock operator with contributions only from the active END_DOC @@ -55,39 +55,42 @@ double precision :: coulomb, exchange double precision :: get_mo_bielec_integral_schwartz integer :: j_act_orb,k_act_orb,i_inact_core_orb + integer :: i_state - do i = 1, n_core_inact_orb - accu_coulomb = 0.d0 - accu_exchange = 0.d0 - i_inact_core_orb = list_core_inact(i) - do j = 1, n_act_orb - j_act_orb = list_act(j) - na = one_body_dm_mo_alpha(j_act_orb,j_act_orb) - nb = one_body_dm_mo_beta(j_act_orb,j_act_orb) - ntot = na + nb - coulomb = mo_bielec_integral_jj(i_inact_core_orb,j_act_orb) - exchange = mo_bielec_integral_jj_exchange(i_inact_core_orb,j_act_orb) - 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) + do i_state = 1,N_states + do i = 1, n_core_inact_orb + accu_coulomb = 0.d0 + accu_exchange = 0.d0 + i_inact_core_orb = list_core_inact(i) + do j = 1, n_act_orb + j_act_orb = list_act(j) + na = one_body_dm_mo_alpha(j_act_orb,j_act_orb,i_state) + nb = one_body_dm_mo_beta(j_act_orb,j_act_orb,i_state) 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 + coulomb = mo_bielec_integral_jj(i_inact_core_orb,j_act_orb) + exchange = mo_bielec_integral_jj_exchange(i_inact_core_orb,j_act_orb) + 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,i_state) + nb = one_body_dm_mo_beta(j_act_orb,k_act_orb,i_state) + 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,i_state) = accu_coulomb + accu_exchange(1) + fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb + accu_exchange(2) 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) enddo END_PROVIDER - BEGIN_PROVIDER [double precision, fock_virt_from_act, (mo_tot_num,2)] + BEGIN_PROVIDER [double precision, fock_virt_from_act, (mo_tot_num,2,N_states)] BEGIN_DOC ! virtual part of the fock operator with contributions only from the active END_DOC @@ -98,67 +101,77 @@ double precision :: coulomb, exchange double precision :: get_mo_bielec_integral_schwartz integer :: j_act_orb,i_virt_orb,k_act_orb + integer :: i_state + ! TODO : inverse loop of i_state - do i = 1, n_virt_orb - accu_coulomb = 0.d0 - accu_exchange = 0.d0 - i_virt_orb = list_virt(i) - do j = 1, n_act_orb - j_act_orb = list_act(j) - na = one_body_dm_mo_alpha(j_act_orb,j_act_orb) - nb = one_body_dm_mo_beta(j_act_orb,j_act_orb) - ntot = na + nb - coulomb = mo_bielec_integral_jj(i_virt_orb,j_act_orb) - exchange = mo_bielec_integral_jj_exchange(i_virt_orb,j_act_orb) - 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) + do i_state = 1, N_states + do i = 1, n_virt_orb + accu_coulomb = 0.d0 + accu_exchange = 0.d0 + i_virt_orb = list_virt(i) + do j = 1, n_act_orb + j_act_orb = list_act(j) + na = one_body_dm_mo_alpha(j_act_orb,j_act_orb,i_state) + nb = one_body_dm_mo_beta(j_act_orb,j_act_orb,i_state) 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 + coulomb = mo_bielec_integral_jj(i_virt_orb,j_act_orb) + exchange = mo_bielec_integral_jj_exchange(i_virt_orb,j_act_orb) + 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,i_state) + nb = one_body_dm_mo_beta(j_act_orb,k_act_orb,i_state) + 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,i_state) = accu_coulomb + accu_exchange(1) + fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb + accu_exchange(2) 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) enddo END_PROVIDER - BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2)] -&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_trace, (mo_tot_num)] + BEGIN_PROVIDER [double precision, fock_core_inactive_total, (mo_tot_num,2,N_states)] +&BEGIN_PROVIDER [double precision, fock_core_inactive_total_spin_trace, (mo_tot_num,N_states)] BEGIN_DOC ! inactive part of the fock operator END_DOC implicit none integer :: i integer :: i_inact_core_orb - do i = 1, n_core_inact_orb - 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_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)) + integer :: i_state + do i_state = 1, N_states + do i = 1, n_core_inact_orb + i_inact_core_orb = list_core_inact(i) + fock_core_inactive_total(i_inact_core_orb,1,i_state) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,1,i_state) + fock_core_inactive_total(i_inact_core_orb,2,i_state) = fock_core_inactive(i_inact_core_orb) + fock_core_inactive_from_act(i_inact_core_orb,2,i_state) + fock_core_inactive_total_spin_trace(i_inact_core_orb,i_state) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1,i_state) + fock_core_inactive_total(i_inact_core_orb,2,i_state)) + enddo enddo END_PROVIDER - BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)] -&BEGIN_PROVIDER [double precision, fock_virt_total_spin_trace, (mo_tot_num)] + BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2,N_states)] +&BEGIN_PROVIDER [double precision, fock_virt_total_spin_trace, (mo_tot_num,N_states)] BEGIN_DOC ! inactive part of the fock operator END_DOC implicit none integer :: i 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)+ 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) ) + integer :: i_state + do i_state = 1, N_states + do i = 1, n_virt_orb + i_virt_orb= list_virt(i) + fock_virt_total(i_virt_orb,1,i_state) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,1,i_state)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb) + fock_virt_total(i_virt_orb,2,i_state) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2,i_state)+ mo_mono_elec_integral(i_virt_orb,i_virt_orb) + fock_virt_total_spin_trace(i_virt_orb,i_state) = 0.5d0 * ( fock_virt_total(i_virt_orb,1,i_state) + fock_virt_total(i_virt_orb,2,i_state) ) + enddo enddo END_PROVIDER diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index c04b14fa..512158cf 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -30,19 +30,19 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip integer :: degree_alpha(psi_det_size) logical :: fullMatch - double precision :: delta_e_array(psi_det_size) + double precision :: delta_e_inv_array(psi_det_size,N_states) double precision :: hij_array(psi_det_size) integer(bit_kind) :: tq(Nint,2,n_selected) integer :: N_tq - double precision :: hialpha + double precision :: hialpha,hij integer :: i_state, i_alpha integer(bit_kind),allocatable :: miniList(:,:,:) integer,allocatable :: idx_miniList(:) integer :: N_miniList, leng - double precision :: delta_e_final,hij_tmp + double precision :: delta_e(N_states),hij_tmp integer :: index_i,index_j @@ -76,12 +76,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip ! coef_pert = 0.d0 do i = 1,idx_alpha(0) index_i = idx_alpha(i) - call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e_final) + call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) - delta_e_array(index_i) = 1.d0/delta_e_final hij_array(index_i) = hialpha - ! ihpsi0 += hialpha * psi_coef(index_i,1) - ! coef_pert += hialpha * psi_coef(index_i,1) * delta_e_array(index_i) + do i_state = 1,N_states + delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) + enddo enddo do i=1,idx_alpha(0) @@ -91,7 +91,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip do j = 1, idx_alpha(0) index_j = idx_alpha(j) do i_state=1,N_states - delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_array(index_j) +! standard dressing first order + delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state) enddo enddo call omp_unset_lock( psi_ref_bis_lock(index_i)) diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 19a44640..3b832b58 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -14,7 +14,8 @@ ! Dressing matrix in N_det basis END_DOC integer :: i,j,m - double precision :: accu + integer :: i_state + double precision :: accu(N_states) double precision, allocatable :: delta_ij_tmp(:,:,:) @@ -27,118 +28,136 @@ delta_ij_tmp = 0.d0 call H_apply_mrpt_1h(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo + enddo + second_order_pt_new_1h(i_state) = accu(i_state) enddo print*, '1h = ',accu - second_order_pt_new_1h(1) = accu ! 1p delta_ij_tmp = 0.d0 call H_apply_mrpt_1p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_1p(i_state) = accu(i_state) + enddo print*, '1p = ',accu - second_order_pt_new_1p(1) = accu ! 1h1p delta_ij_tmp = 0.d0 call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_1h1p(i_state) = accu(i_state) + enddo print*, '1h1p = ',accu - second_order_pt_new_1h1p(1) = accu ! 2h delta_ij_tmp = 0.d0 call H_apply_mrpt_2h(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_2h(i_state) = accu(i_state) + enddo print*, '2h = ',accu - second_order_pt_new_2h(1) = accu ! 2p delta_ij_tmp = 0.d0 call H_apply_mrpt_2p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_2p(i_state) = accu(i_state) + enddo print*, '2p = ',accu - second_order_pt_new_2p(1) = accu ! 1h2p delta_ij_tmp = 0.d0 call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_1h2p(i_state) = accu(i_state) + enddo print*, '1h2p = ',accu - second_order_pt_new_1h2p(1) = accu ! 2h1p delta_ij_tmp = 0.d0 call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_2h1p(i_state) = accu(i_state) + enddo print*, '2h1p = ',accu - second_order_pt_new_2h1p(1) = accu ! 2h2p delta_ij_tmp = 0.d0 call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det do j = 1, N_det - accu += delta_ij_tmp(j,i,1) * psi_coef(i,1) * psi_coef(j,1) - delta_ij(j,i,1) += delta_ij_tmp(j,i,1) + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo + second_order_pt_new_2h2p(i_state) = accu(i_state) + enddo print*, '2h2p = ',accu - second_order_pt_new_2h2p(1) = accu ! total accu = 0.d0 + do i_state = 1, N_states do i = 1, N_det - do j = 1, N_det - accu += delta_ij(j,i,1) * psi_coef(i,1) * psi_coef(j,1) + write(*,'(1000(F16.10,x))')delta_ij(i,:,:) + do j = i_state, N_det + accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) enddo enddo - print*, 'total= ',accu - second_order_pt_new(1) = accu + second_order_pt_new(i_state) = accu(i_state) + print*, 'total= ',accu(i_state) + enddo -! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) END_PROVIDER diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 02d11244..8d705deb 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -152,33 +152,51 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti end subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) + BEGIN_DOC + ! routine that returns the delta_e with the Moller Plesset and Dyall operators + ! + ! with det_1 being a determinant from the cas, and det_2 being a perturber + ! + ! Delta_e(det_1,det_2) = sum (hole) epsilon(hole) + sum(part) espilon(part) + delta_e(act) + ! + ! where hole is necessary in the inactive, part necessary in the virtuals + ! + ! and delta_e(act) is obtained from the contracted application of the excitation + ! + ! operator in the active space that lead from det_1 to det_2 + END_DOC implicit none use bitmasks - double precision, intent(out) :: delta_e_final + double precision, intent(out) :: delta_e_final(N_states) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) integer :: i,j,k,l + integer :: i_state integer :: n_holes_spin(2) integer :: n_holes integer :: holes_list(N_int*bit_kind_size,2) - double precision :: delta_e_inactive + double precision :: delta_e_inactive(N_states) integer :: i_hole_inact call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) 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_trace(i_hole_inact) + do i_state = 1, N_states + delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) + enddo enddo do i = 1, n_holes_spin(2) i_hole_inact = holes_list(i,2) - delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact) + do i_state = 1, N_states + delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) + enddo enddo - double precision :: delta_e_virt + double precision :: delta_e_virt(N_states) integer :: i_part_virt integer :: n_particles_spin(2) integer :: n_particles @@ -188,12 +206,16 @@ 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_trace(i_part_virt) + do i_state = 1, N_states + delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state) + enddo 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) + do i_state = 1, N_states + delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state) + enddo enddo @@ -203,7 +225,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) integer :: holes_active_list_spin_traced(4*n_act_orb) integer :: particles_active_list(2*n_act_orb,2) integer :: particles_active_list_spin_traced(4*n_act_orb) - double precision :: delta_e_act + double precision :: delta_e_act(N_states) 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) @@ -265,14 +287,18 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) ! 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) + do i_state = 1, N_states + delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) + enddo 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) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) + enddo else if (n_holes_act == 1 .and. n_particles_act == 1) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -284,7 +310,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) ! 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) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil_one_creat(i_particle_act,i_hole_act,jspin,ispin,i_state) + enddo else if (n_holes_act == 2 .and. n_particles_act == 0) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -294,7 +322,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) 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) + do i_state = 1, N_states + delta_e_act(i_state) += two_anhil(i_hole_act,j_hole_act,ispin,jspin,i_state) + enddo else if (n_holes_act == 0 .and. n_particles_act == 2) then ! i_particle_act = particles_active_list_spin_traced(1) @@ -304,7 +334,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) 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) + do i_state = 1, N_states + delta_e_act(i_state) += two_creat(i_particle_act,j_particle_act,ispin,jspin,i_state) + enddo else if (n_holes_act == 2 .and. n_particles_act == 1) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -324,7 +356,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) ! 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) + do i_state = 1, N_states + delta_e_act(i_state) += two_anhil_one_creat(i_particle_act,i_hole_act,j_hole_act,kspin,ispin,jspin,i_state) + enddo else if (n_holes_act == 1 .and. n_particles_act == 2) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -342,7 +376,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) 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) + do i_state = 1, N_states + delta_e_act(i_state) += two_creat_one_anhil(i_particle_act,j_particle_act,i_hole_act,jspin,kspin,ispin,i_state) + enddo else if (n_holes_act == 3 .and. n_particles_act == 0) then ! i_hole_act = holes_active_list_spin_traced(1) @@ -359,7 +395,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) ! 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) + do i_state = 1, N_states + delta_e_act(i_state) += three_anhil(i_hole_act,j_hole_act,k_hole_act,ispin,jspin,kspin,i_state) + enddo else if (n_holes_act == 0 .and. n_particles_act == 3) then ! i_particle_act = particles_active_list_spin_traced(1) @@ -376,7 +414,9 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) 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) + do i_state = 1, N_states + delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state) + enddo else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then @@ -388,265 +428,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) !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 .or. delta_e_final == 0.d0)then -!if(delta_e_final == 0.d0)then - if(.False.)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 - delta_e_final = -1000.d0 -!stop - - endif - -end - -subroutine get_delta_e_dyall_verbose(det_1,det_2,delta_e_final) - implicit none - use bitmasks - double precision, intent(out) :: delta_e_final - integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - integer :: i,j,k,l - - integer :: n_holes_spin(2) - integer :: n_holes - integer :: holes_list(N_int*bit_kind_size,2) - - - double precision :: delta_e_inactive - integer :: i_hole_inact - - call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) - 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_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_trace(i_hole_inact) - enddo - - double precision :: delta_e_virt - integer :: i_part_virt - integer :: n_particles_spin(2) - integer :: n_particles - integer :: particles_list(N_int*bit_kind_size,2) - - call give_particles_in_virt_space(det_2,n_particles_spin,n_particles,particles_list) - 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_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 - - - integer :: n_holes_spin_act(2),n_particles_spin_act(2) - integer :: n_holes_act,n_particles_act - integer :: holes_active_list(2*n_act_orb,2) - integer :: holes_active_list_spin_traced(4*n_act_orb) - integer :: particles_active_list(2*n_act_orb,2) - integer :: particles_active_list_spin_traced(4*n_act_orb) - double precision :: delta_e_act - 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,icountbis - integer :: hole_list_practical(2,elec_num_tab(1)+elec_num_tab(2)), particle_list_practical(2,elec_num_tab(1)+elec_num_tab(2)) - 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_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_act) then - print*, icount, n_particles_act - print * , 'pb in particles_active_list_spin_traced !!' - stop - endif - - - integer :: i_hole_act, j_hole_act, k_hole_act - integer :: i_particle_act, j_particle_act, k_particle_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 == 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_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) -! 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_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_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_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 - -!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 .or. delta_e_final == 0.d0)then -!if(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 - delta_e_final = -1000.d0 + do i_state = 1, n_states + delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state) + enddo end diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 6619b125..5cd09aa2 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -7,12 +7,6 @@ use bitmasks integer :: i number_of_holes = 0 - do i = 1, N_int - number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,1),iand(key_in(i,1),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) ) - enddo - return if(N_int == 1)then number_of_holes = number_of_holes & + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 2253c33c..118bbdf7 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -1,5 +1,22 @@ - BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] -&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_average, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_average, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: i + + one_body_dm_mo_alpha_average = 0.d0 + one_body_dm_mo_beta_average = 0.d0 + do i = 1,N_states + one_body_dm_mo_alpha_average(:,:) += one_body_dm_mo_alpha(:,:,i) * state_average_weight(i) + one_body_dm_mo_beta_average(:,:) += one_body_dm_mo_beta(:,:,i) * state_average_weight(i) + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num,N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num,N_states) ] implicit none BEGIN_DOC ! Alpha and beta one-body density matrix for each state @@ -11,36 +28,31 @@ double precision :: phase integer :: h1,h2,p1,p2,s1,s2, degree integer :: exc(0:2,2,2),n_occ(2) - double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:) - if(only_single_double_dm)then - print*,'ONLY DOUBLE DM' - one_body_dm_mo_alpha = one_body_single_double_dm_mo_alpha - one_body_dm_mo_beta = one_body_single_double_dm_mo_beta - else one_body_dm_mo_alpha = 0.d0 one_body_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & !$OMP tmp_a, tmp_b, n_occ)& - !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,& !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& !$OMP mo_tot_num) - allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) ) tmp_a = 0.d0 tmp_b = 0.d0 !$OMP DO SCHEDULE(dynamic) do k=1,N_det call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int) do m=1,N_states - ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + ck = psi_coef(k,m)*psi_coef(k,m) do l=1,elec_alpha_num j = occ(l,1) - tmp_a(j,j) += ck + tmp_a(j,j,m) += ck enddo do l=1,elec_beta_num j = occ(l,2) - tmp_b(j,j) += ck + tmp_b(j,j,m) += ck enddo enddo do l=1,k-1 @@ -51,28 +63,27 @@ call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) do m=1,N_states - ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + ckl = psi_coef(k,m) * psi_coef(l,m) * phase if (s1==1) then - tmp_a(h1,p1) += ckl - tmp_a(p1,h1) += ckl + tmp_a(h1,p1,m) += ckl + tmp_a(p1,h1,m) += ckl else - tmp_b(h1,p1) += ckl - tmp_b(p1,h1) += ckl + tmp_b(h1,p1,m) += ckl + tmp_b(p1,h1,m) += ckl endif enddo enddo enddo !$OMP END DO NOWAIT !$OMP CRITICAL - one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:) !$OMP END CRITICAL !$OMP CRITICAL - one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + one_body_dm_mo_beta(:,:,:) = one_body_dm_mo_beta(:,:,:) + tmp_b(:,:,:) !$OMP END CRITICAL deallocate(tmp_a,tmp_b) !$OMP END PARALLEL - endif END_PROVIDER BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] @@ -163,7 +174,7 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) BEGIN_DOC ! One-body density matrix END_DOC - one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta + one_body_dm_mo = one_body_dm_mo_alpha_average + one_body_dm_mo_beta_average END_PROVIDER BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,mo_tot_num) ] @@ -171,7 +182,7 @@ BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,m BEGIN_DOC ! rho(alpha) - rho(beta) END_DOC - one_body_spin_density_mo = one_body_dm_mo_alpha - one_body_dm_mo_beta + one_body_spin_density_mo = one_body_dm_mo_alpha_average - one_body_dm_mo_beta_average END_PROVIDER subroutine set_natural_mos @@ -246,8 +257,8 @@ END_PROVIDER do l = 1, ao_num do i = 1, mo_tot_num do j = 1, mo_tot_num - mo_alpha = one_body_dm_mo_alpha(j,i) - mo_beta = one_body_dm_mo_beta(j,i) + mo_alpha = one_body_dm_mo_alpha_average(j,i) + mo_beta = one_body_dm_mo_beta_average(j,i) ! if(dabs(dm_mo).le.1.d-10)cycle one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_beta