From 38ccfc0cf1aabfa504049da09d51cbe6623ad56f Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 28 Nov 2016 14:51:12 +0100 Subject: [PATCH] Trying do really fo sin free multiple excitations --- plugins/Full_CI/H_apply.irp.f | 5 - plugins/MRCC_Utils/mrcc_utils.irp.f | 9 +- plugins/MRPT/print_1h2p.irp.f | 46 +-- plugins/MRPT_Utils/energies_cas.irp.f | 82 +++--- plugins/MRPT_Utils/mrpt_dress.irp.f | 21 +- plugins/MRPT_Utils/mrpt_utils.irp.f | 267 +++++++++--------- plugins/MRPT_Utils/psi_active_prov.irp.f | 44 ++- .../pt2_new.irp.f | 0 plugins/Perturbation/NEEDED_CHILDREN_MODULES | 2 +- plugins/Perturbation/pt2_equations.irp.f | 30 -- 10 files changed, 235 insertions(+), 271 deletions(-) rename plugins/{Perturbation => MRPT_Utils}/pt2_new.irp.f (100%) diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 79599065..8977b7fd 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -12,11 +12,6 @@ s.set_perturbation("epstein_nesbet_2x2") s.unset_openmp() print s -s = H_apply("FCI_PT2_new") -s.set_perturbation("decontracted") -s.unset_openmp() -print s - s = H_apply("FCI_no_skip") s.set_selection_pt2("epstein_nesbet_2x2") diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index c1a277cf..92425eb0 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -635,34 +635,37 @@ END_PROVIDER double precision :: coef_mrpt(N_States),coef_array(N_states),hij,delta_e(N_states) double precision :: hij_array(N_det_Ref),delta_e_array(N_det_ref,N_states) do i = 1, N_det_non_ref + print*,'i',i do j = 1, N_det_ref do k = 1, N_States - coef_array(j) = psi_ref_coef(j,k) + coef_array(k) = psi_ref_coef(j,k) enddo call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,i), N_int, Hij_array(j)) call get_delta_e_dyall(psi_ref(1,1,j),psi_non_ref(1,1,i),coef_array,hij_array(j),delta_e) - print*,delta_e(:) +! write(*,'(A7,x,100(F16.10,x))')'delta_e',delta_e(:) do k = 1, N_states delta_e_Array(j,k) = delta_e(k) enddo enddo + coef_mrpt = 0.d0 do k = 1, N_states do j = 1, N_det_Ref coef_mrpt(k) += psi_ref_coef(j,k) * hij_array(j) / delta_e_array(j,k) enddo enddo + write(*,'(A7,X,100(F16.10,x))')'coef ',psi_non_ref_coef(i,1) , coef_mrpt(1),psi_non_ref_coef(i,2) , coef_mrpt(2) do k = 1, N_States if(dabs(coef_mrpt(k)) .le.1.d-10)then rho_mrpt(i,k) = 0.d0 exit endif - print*,k,psi_non_ref_coef(i,k) , coef_mrpt(k) if(psi_non_ref_coef(i,k) / coef_mrpt(k) .lt.0d0)then rho_mrpt(i,k) = 1.d0 else rho_mrpt(i,k) = psi_non_ref_coef(i,k) / coef_mrpt(k) endif enddo + print*,'rho',rho_mrpt(i,:) enddo END_PROVIDER diff --git a/plugins/MRPT/print_1h2p.irp.f b/plugins/MRPT/print_1h2p.irp.f index 747e2817..a3500e49 100644 --- a/plugins/MRPT/print_1h2p.irp.f +++ b/plugins/MRPT/print_1h2p.irp.f @@ -2,7 +2,7 @@ program print_1h2p implicit none read_wf = .True. touch read_wf - call routine_2 + call routine end subroutine routine_2 @@ -20,44 +20,20 @@ end subroutine routine implicit none double precision,allocatable :: matrix_1h2p(:,:,:) - allocate (matrix_1h2p(N_det,N_det,N_states)) + double precision :: accu(2) + allocate (matrix_1h2p(N_det_ref,N_det_ref,N_states)) integer :: i,j,istate - do i = 1, N_det - do j = 1, N_det - do istate = 1, N_states - matrix_1h2p(i,j,istate) = 0.d0 + accu = 0.d0 + matrix_1h2p = 0.d0 + call H_apply_mrpt_2p(matrix_1h2p,N_det_ref) + do istate = 1, N_states + do i = 1, N_det + do j = 1, N_det + accu(istate) += matrix_1h2p(i,j,istate) * psi_coef(i,istate) * psi_coef(j,istate) enddo enddo + print*,accu(istate) enddo - if(.False.)then - call give_1h2p_contrib(matrix_1h2p) - double precision :: accu - accu = 0.d0 - do i = 1, N_det - do j = 1, N_det - accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1) - enddo - enddo - print*, 'second order ', accu - endif - - if(.True.)then - do i = 1, N_det - do j = 1, N_det - do istate = 1, N_states - matrix_1h2p(i,j,istate) = 0.d0 - enddo - enddo - enddo - call give_1h2p_new(matrix_1h2p) - accu = 0.d0 - do i = 1, N_det - do j = 1, N_det - accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1) - enddo - enddo - endif - print*, 'third order ', accu deallocate (matrix_1h2p) end diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index 54e1a3f8..33baeb6a 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] double precision :: norm_out(N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) use bitmasks integer :: iorb @@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] orb = list_act(iorb) hole_particle = 1 spin_exc = ispin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -54,8 +54,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] 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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + one_creat(iorb,ispin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -72,7 +72,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb integer :: state_target @@ -82,7 +82,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] orb = list_act(iorb) hole_particle = -1 spin_exc = ispin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -94,8 +94,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] 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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + one_anhil(iorb,ispin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -113,7 +113,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target @@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) orb_j = list_act(jorb) hole_particle_j = 1 spin_exc_j = jspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -142,8 +142,8 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -163,7 +163,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target @@ -179,7 +179,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) orb_j = list_act(jorb) hole_particle_j = -1 spin_exc_j = jspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -192,8 +192,8 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -213,7 +213,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target double precision :: energies(n_states) @@ -227,7 +227,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 orb_j = list_act(jorb) hole_particle_j = -1 spin_exc_j = jspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -289,7 +289,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a orb_k = list_act(korb) hole_particle_k = -1 spin_exc_k = kspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -306,8 +306,8 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -330,7 +330,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -351,7 +351,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a orb_k = list_act(korb) hole_particle_k = -1 spin_exc_k = kspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -367,8 +367,8 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -391,7 +391,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -412,7 +412,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 orb_k = list_act(korb) hole_particle_k = 1 spin_exc_k = kspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -428,8 +428,8 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -452,7 +452,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -473,7 +473,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 orb_k = list_act(korb) hole_particle_k = -1 spin_exc_k = kspin - do i = 1, n_det + do i = 1, n_det_ref do j = 1, n_states psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo @@ -489,8 +489,8 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,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_ref,n_det_ref,n_det_ref,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) - three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) enddo enddo enddo @@ -515,7 +515,7 @@ END_PROVIDER integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb,i_ok integer :: state_target @@ -541,7 +541,7 @@ END_PROVIDER do state_target =1 , N_states one_anhil_one_creat_inact_virt_norm(iorb,vorb,state_target,ispin) = 0.d0 enddo - do i = 1, n_det + do i = 1, n_det_ref do j = 1, N_int psi_in_out(j,1,i) = psi_ref(j,1,i) psi_in_out(j,2,i) = psi_ref(j,2,i) @@ -616,7 +616,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: jorb,i_ok,aorb,orb_a integer :: state_target @@ -641,7 +641,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta norm = 0.d0 norm_bis = 0.d0 do ispin = 1,2 - do i = 1, n_det + do i = 1, n_det_ref do j = 1, N_int psi_in_out(j,1,i) = psi_ref(j,1,i) psi_in_out(j,2,i) = psi_ref(j,2,i) @@ -715,7 +715,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb,i_ok,aorb,orb_a integer :: state_target @@ -740,7 +740,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State norm = 0.d0 norm_bis = 0.d0 do ispin = 1,2 - do i = 1, n_det + do i = 1, n_det_ref do j = 1, N_int psi_in_out(j,1,i) = psi_ref(j,1,i) psi_in_out(j,2,i) = psi_ref(j,2,i) @@ -825,7 +825,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:) double precision, allocatable :: delta_e_det(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det_ref+1,N_det_ref+1)) + allocate (psi_in_out(N_int,2,n_det_ref),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det_ref+1,N_det_ref+1)) allocate (eigenvectors(size(H_matrix,1),N_det_ref+1)) allocate (eigenvalues(N_det_ref+1),interact_cas(N_det_ref,N_det_ref)) allocate (delta_e_det(N_det_ref,N_det_ref)) @@ -854,7 +854,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from - fock_virt_total_spin_trace(orb_v,j) enddo do ispin = 1,2 - do i = 1, n_det + do i = 1, n_det_ref do j = 1, N_int psi_in_out(j,1,i) = psi_ref(j,1,i) psi_in_out(j,2,i) = psi_ref(j,2,i) diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index f5e7bd40..eccdae0a 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -84,14 +84,21 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip do i_state = 1, N_states coef_array(i_state) = psi_ref_coef(index_i,i_state) enddo - if(dabs(hialpha).le.1.d-10)then + integer :: degree_scalar + + call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,index_i),degree_scalar,N_int) + if(dabs(hialpha).le.1.d-20)then delta_e = 1.d+20 else call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) endif + if(degree .ne. 2)then + do i_state = 1, N_states + delta_e(i_state) = 1.d+20 + enddo + endif hij_array(index_i) = hialpha call get_excitation(psi_ref(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) -! phase_array(index_i) = phase do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) enddo @@ -103,16 +110,6 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call omp_set_lock( psi_ref_bis_lock(index_i) ) do j = 1, idx_alpha(0) index_j = idx_alpha(j) -! call get_excitation(psi_ref(1,1,index_i),psi_ref(1,1,index_i),exc,degree,phase,N_int) -! if(index_j.ne.index_i)then -! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then -! print*, phase_array(index_j) , phase_array(index_i) ,phase -! call debug_det(psi_ref(1,1,index_i),N_int) -! call debug_det(psi_ref(1,1,index_j),N_int) -! call debug_det(tq(1,1,i_alpha),N_int) -! stop -! endif -! endif do i_state=1,N_states ! 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) diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 31013bb7..cb007199 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -39,155 +39,141 @@ enddo print*, '1h = ',accu - ! 1p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 1p +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 1h1p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref) - double precision :: e_corr_from_1h1p_singles(N_states) -!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) -!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 1h1p +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref) +!double precision :: e_corr_from_1h1p_singles(N_states) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 1h1p third order - if(do_third_order_1h1p)then - delta_ij_tmp = 0.d0 - call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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(3)',accu - endif +!! 1h1p third order +!if(do_third_order_1h1p)then +! delta_ij_tmp = 0.d0 +! call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) +! accu = 0.d0 +! do i_state = 1, N_states +! do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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(3)',accu +!endif - ! 2h - delta_ij_tmp = 0.d0 - call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 2h +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 2p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 2p +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 1h2p - delta_ij_tmp = 0.d0 - call give_1h2p_contrib(delta_ij_tmp) +!! 1h2p +!delta_ij_tmp = 0.d0 +!call give_1h2p_contrib(delta_ij_tmp) !!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 2h1p - delta_ij_tmp = 0.d0 - call give_2h1p_contrib(delta_ij_tmp) - !!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 2h1p +!delta_ij_tmp = 0.d0 +!call give_2h1p_contrib(delta_ij_tmp) +!!!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det_ref +! do j = 1, N_det_ref +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 - ! 2h2p - delta_ij_tmp = 0.d0 -!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref) - accu = 0.d0 - do i_state = 1, N_states - do i = 1, N_det_ref - do j = 1, N_det_ref - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_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 +!! 2h2p - double precision :: contrib_2h2p(N_states) - call give_2h2p(contrib_2h2p) - do i_state = 1, N_states - do i = 1, N_det_ref - delta_ij(i,i,i_state) += contrib_2h2p(i_state) - enddo - second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) - enddo - print*, '2h2p = ',contrib_2h2p(1) +!double precision :: contrib_2h2p(N_states) +!call give_2h2p(contrib_2h2p) +!do i_state = 1, N_states +! do i = 1, N_det_ref +! delta_ij(i,i,i_state) += contrib_2h2p(i_state) +! enddo +!second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) +!enddo +!print*, '2h2p = ',contrib_2h2p(:) ! total - accu = 0.d0 + accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det_ref -! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) - do j = i_state, N_det_ref - accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) + print*,'state ',i_state + do i = 1, N_det_ref + write(*,'(1000(F16.10,x))')delta_ij(i,:,i_state) + do j = i , N_det_ref + accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) + enddo enddo - enddo - second_order_pt_new(i_state) = accu(i_state) - print*, 'total= ',accu(i_state) + second_order_pt_new(i_state) = accu(i_state) + print*, 'total= ',accu(i_state) enddo @@ -217,7 +203,7 @@ END_PROVIDER double precision :: hij do i_state = 1, N_states do i = 1,N_det_ref - do j = i,N_det_ref + do j = 1,N_det_ref call i_h_j(psi_ref(1,1,j),psi_ref(1,1,i),N_int,hij) Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = hij & + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) ) @@ -284,9 +270,9 @@ END_PROVIDER do j = 1, N_det_ref hmatrix_tmp(j,i) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo +! print*,i,hmatrix_tmp(i,i)+nuclear_repulsion enddo - call lapack_diag(eigenvalues,eigenvectors, & - Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det_ref,N_det_ref) + call lapack_diag(eigenvalues,eigenvectors,hmatrix_tmp,N_det_ref,N_det_ref) write(*,'(A86)')'Looking for the most overlapping state within all eigenvectors of the dressed matrix' print*,'' print*,'Calculating the overlap for ...' @@ -303,7 +289,10 @@ END_PROVIDER enddo print*,'' print*,'Sorting the eigenvectors per overlap' - call dsort(overlap,iorder,n_states) + call dsort(overlap,iorder,n_det_ref) + do j = 1, N_det_ref + print*,overlap(j),iorder(j) + enddo print*,'' print*,'The most overlapping state is the ',iorder(1) print*,'with the overlap of ',dabs(overlap(1)) diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 33cb5d5b..e3b0986a 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -180,7 +180,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) double precision :: delta_e_inactive(N_states) - integer :: i_hole_inact + integer :: i_hole_inact, list_holes_inact(n_inact_orb,2) call get_excitation_degree(det_1,det_2,degree,N_int) if(degree>2)then @@ -190,8 +190,13 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) delta_e_inactive = 0.d0 + integer :: n_holes_total + n_holes_total = 0 do i = 1, n_holes_spin(1) i_hole_inact = holes_list(i,1) + n_holes_total +=1 + list_holes_inact(n_holes_total,1) = i_hole_inact + list_holes_inact(n_holes_total,2) = 1 do i_state = 1, N_states delta_e_inactive += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) enddo @@ -199,6 +204,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) do i = 1, n_holes_spin(2) i_hole_inact = holes_list(i,2) + n_holes_total +=1 + list_holes_inact(n_holes_total,1) = i_hole_inact + list_holes_inact(n_holes_total,2) = 2 do i_state = 1, N_states delta_e_inactive(i_state) += fock_core_inactive_total_spin_trace(i_hole_inact,i_state) enddo @@ -370,15 +378,41 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) enddo else if (n_holes_act == 1 .and. n_particles_act == 2) then + ! First find the particle that has been added from the inactive + ! + integer :: spin_hole_inact, spin_hole_part_act + spin_hole_inact = list_holes_inact(1,2) +! spin_hole_part_act = + if(jspin == spin_hole_inact )then + kspin = spin_hole_part_act + ispin = spin_hole_part_act + else + jspin = spin_hole_part_act + ispin = spin_hole_part_act + endif + ! by convention, you first make a movement in the cas ! first hole - ispin = hole_list_practical(1,1) - i_hole_act = hole_list_practical(2,1) + i_hole_act = hole_list_practical(2,1) + jspin = spin_hole_inact ! 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) + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 2)then + print*, '' + call debug_det(det_1,N_int) + call debug_det(det_2,N_int) + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*, s1,h1,p1 + print*, s2,h2,p2 + print*, '---' + print*, ispin, i_hole_act + print*, jspin, i_particle_act + print*, kspin, j_particle_act + pause + endif 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) diff --git a/plugins/Perturbation/pt2_new.irp.f b/plugins/MRPT_Utils/pt2_new.irp.f similarity index 100% rename from plugins/Perturbation/pt2_new.irp.f rename to plugins/MRPT_Utils/pt2_new.irp.f diff --git a/plugins/Perturbation/NEEDED_CHILDREN_MODULES b/plugins/Perturbation/NEEDED_CHILDREN_MODULES index 25b89c5f..f7999340 100644 --- a/plugins/Perturbation/NEEDED_CHILDREN_MODULES +++ b/plugins/Perturbation/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants Properties Hartree_Fock Davidson MRPT_Utils +Determinants Properties Hartree_Fock Davidson diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index b29e130f..5839c20c 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -46,36 +46,6 @@ 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) - H_pert_diag = 0.d0 - - - 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