diff --git a/plugins/MRPT/MRPT_Utils.main.irp.f b/plugins/MRPT/jm_mrpt2.irp.f similarity index 100% rename from plugins/MRPT/MRPT_Utils.main.irp.f rename to plugins/MRPT/jm_mrpt2.irp.f diff --git a/plugins/MRPT/mrpt.irp.f b/plugins/MRPT/mrpt.irp.f deleted file mode 100644 index 8c8d746d..00000000 --- a/plugins/MRPT/mrpt.irp.f +++ /dev/null @@ -1,38 +0,0 @@ -program MRPT - implicit none - BEGIN_DOC -! TODO - END_DOC - print *, ' _/ ' - print *, ' -:\_?, _Jm####La ' - print *, 'J"(:" > _]#AZ#Z#UUZ##, ' - print *, '_,::./ %(|i%12XmX1*1XL _?, ' - print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( ' - print *, ' .:< ]J=mQD?WXn|,)nr" ' - print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" ' - print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 ' - print *, ' "23Z#1S2oo2XXSnnnoSo2>v" ' - print *, ' miX#L -~`""!!1}oSoe|i7 ' - print *, ' 4cn#m, v221=|v[ ' - print *, ' ]hI3Zma,;..__wXSe=+vo ' - print *, ' ]Zov*XSUXXZXZXSe||vo2 ' - print *, ' ]Z#>=|< ' - print *, ' -ziiiii||||||+||==+> ' - print *, ' -%|+++||=|=+|=|==/ ' - print *, ' -a>====+|====-:- ' - print *, ' "~,- -- /- ' - print *, ' -. )> ' - print *, ' .~ +- ' - print *, ' . .... : . ' - print *, ' -------~ ' - print *, '' -end diff --git a/plugins/MRPT_Utils/density_matrix_based.irp.f b/plugins/MRPT_Utils/density_matrix_based.irp.f new file mode 100644 index 00000000..b2f3b8cf --- /dev/null +++ b/plugins/MRPT_Utils/density_matrix_based.irp.f @@ -0,0 +1,193 @@ +subroutine contrib_1h2p_dm_based(accu) + implicit none + integer :: i_i,i_r,i_v,i_a,i_b + integer :: i,r,v,a,b + integer :: ispin,jspin + integer :: istate + double precision, intent(out) :: accu(N_states) + double precision :: active_int(n_act_orb,2) + double precision :: delta_e(n_act_orb,2,N_states) + double precision :: get_mo_bielec_integral + accu = 0.d0 +!do i_i = 1, 1 + do i_i = 1, n_inact_orb + i = list_inact(i_i) +! do i_r = 1, 1 + do i_r = 1, n_virt_orb + r = list_virt(i_r) +! do i_v = 1, 1 + do i_v = 1, n_virt_orb + v = list_virt(i_v) + do i_a = 1, n_act_orb + a = list_act(i_a) + active_int(i_a,1) = get_mo_bielec_integral(i,a,r,v,mo_integrals_map) ! direct + active_int(i_a,2) = get_mo_bielec_integral(i,a,v,r,mo_integrals_map) ! exchange + do istate = 1, N_states + do jspin=1, 2 + delta_e(i_a,jspin,istate) = one_anhil(i_a,jspin,istate) & + - fock_virt_total_spin_trace(r,istate) & + - fock_virt_total_spin_trace(v,istate) & + + fock_core_inactive_total_spin_trace(i,istate) + delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate) + enddo + enddo + enddo + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb +! do i_b = i_a, i_a + b = list_act(i_b) + do ispin = 1, 2 ! spin of (i --> r) + do jspin = 1, 2 ! spin of (a --> v) + if(ispin == jspin .and. r.le.v)cycle ! condition not to double count + do istate = 1, N_states + if(ispin == jspin)then + accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_mo_spin_index(a,b,istate,ispin) & + * (active_int(i_b,1) - active_int(i_b,2)) & + * delta_e(i_a,jspin,istate) + else + accu(istate) += active_int(i_a,1) * one_body_dm_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & + * active_int(i_b,1) + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + +end + +subroutine contrib_2h1p_dm_based(accu) + implicit none + integer :: i_i,i_j,i_v,i_a,i_b + integer :: i,j,v,a,b + integer :: ispin,jspin + integer :: istate + double precision, intent(out) :: accu(N_states) + double precision :: active_int(n_act_orb,2) + double precision :: delta_e(n_act_orb,2,N_states) + double precision :: get_mo_bielec_integral + accu = 0.d0 + do i_i = 1, n_inact_orb + i = list_inact(i_i) + do i_j = 1, n_inact_orb + j = list_inact(i_j) + do i_v = 1, n_virt_orb + v = list_virt(i_v) + do i_a = 1, n_act_orb + a = list_act(i_a) + active_int(i_a,1) = get_mo_bielec_integral(i,j,v,a,mo_integrals_map) ! direct + active_int(i_a,2) = get_mo_bielec_integral(i,j,a,v,mo_integrals_map) ! exchange + do istate = 1, N_states + do jspin=1, 2 +! delta_e(i_a,jspin,istate) = +! + delta_e(i_a,jspin,istate) = one_creat(i_a,jspin,istate) - fock_virt_total_spin_trace(v,istate) & + + fock_core_inactive_total_spin_trace(i,istate) & + + fock_core_inactive_total_spin_trace(j,istate) + delta_e(i_a,jspin,istate) = 1.d0/delta_e(i_a,jspin,istate) + enddo + enddo + enddo + do i_a = 1, n_act_orb + a = list_act(i_a) + do i_b = 1, n_act_orb +! do i_b = i_a, i_a + b = list_act(i_b) + do ispin = 1, 2 ! spin of (i --> v) + do jspin = 1, 2 ! spin of (j --> a) + if(ispin == jspin .and. i.le.j)cycle ! condition not to double count + do istate = 1, N_states + if(ispin == jspin)then + accu(istate) += (active_int(i_a,1) - active_int(i_a,2)) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) & + * (active_int(i_b,1) - active_int(i_b,2)) & + * delta_e(i_a,jspin,istate) + else + accu(istate) += active_int(i_a,1) * one_body_dm_dagger_mo_spin_index(a,b,istate,ispin) * delta_e(i_a,ispin,istate) & + * active_int(i_b,1) + endif + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + +end + + +!subroutine contrib_2p_dm_based(accu) +!implicit none +!integer :: i_r,i_v,i_a,i_b,i_c,i_d +!integer :: r,v,a,b,c,d +!integer :: ispin,jspin +!integer :: istate +!double precision, intent(out) :: accu(N_states) +!double precision :: active_int(n_act_orb,n_act_orb,2) +!double precision :: delta_e(n_act_orb,n_act_orb,2,2,N_states) +!double precision :: get_mo_bielec_integral +!accu = 0.d0 +!do i_r = 1, n_virt_orb +! r = list_virt(i_r) +! do i_v = 1, n_virt_orb +! v = list_virt(i_v) +! do i_a = 1, n_act_orb +! a = list_act(i_a) +! do i_b = 1, n_act_orb +! b = list_act(i_b) +! active_int(i_a,i_b,1) = get_mo_bielec_integral(a,b,r,v,mo_integrals_map) ! direct +! active_int(i_a,i_b,2) = get_mo_bielec_integral(a,b,v,r,mo_integrals_map) ! direct +! do istate = 1, N_states +! do jspin=1, 2 ! spin of i_a +! do ispin = 1, 2 ! spin of i_b +! delta_e(i_a,i_b,jspin,ispin,istate) = two_anhil(i_a,i_b,jspin,ispin,istate) & +! - fock_virt_total_spin_trace(r,istate) & +! - fock_virt_total_spin_trace(v,istate) +! delta_e(i_a,i_b,jspin,ispin,istate) = 1.d0/delta_e(i_a,i_b,jspin,ispin,istate) +! enddo +! enddo +! enddo +! enddo +! enddo +! ! diagonal terms +! do i_a = 1, n_act_orb +! a = list_act(i_a) +! do i_b = 1, n_act_orb +! b = list_act(i_b) +! do ispin = 1, 2 ! spin of (a --> r) +! do jspin = 1, 2 ! spin of (b --> v) +! if(ispin == jspin .and. r.le.v)cycle ! condition not to double count +! if(ispin == jspin .and. a.le.b)cycle ! condition not to double count +! do istate = 1, N_states +! if(ispin == jspin)then +! double precision :: contrib_spin +! if(ispin == 1)then +! contrib_spin = two_body_dm_aa_diag_act(i_a,i_b) +! else +! contrib_spin = two_body_dm_bb_diag_act(i_a,i_b) +! endif +! accu(istate) += (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) * contrib_spin & +! * (active_int(i_a,i_b,1) - active_int(i_a,i_b,2)) & +! * delta_e(i_a,i_b,ispin,jspin,istate) +! else +! accu(istate) += 0.5d0 * active_int(i_a,i_b,1) * two_body_dm_ab_diag_act(i_a,i_b) * delta_e(i_a,i_b,ispin,jspin,istate) & +! * active_int(i_a,i_b,1) +! endif +! enddo +! enddo +! enddo +! enddo +! enddo +! enddo +! enddo + + +!end + diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index dd79edbe..30429053 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -1,9 +1,9 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] implicit none integer :: i - double precision :: energies(N_states_diag) + double precision :: energies(N_states) do i = 1, N_states - call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i) + call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i) energy_cas_dyall(i) = energies(i) print*, 'energy_cas_dyall(i)', energy_cas_dyall(i) enddo @@ -13,9 +13,9 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)] implicit none integer :: i - double precision :: energies(N_states_diag) + double precision :: energies(N_states) do i = 1, N_states - call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i) + call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i) energy_cas_dyall_no_exchange(i) = energies(i) print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i) enddo @@ -28,22 +28,22 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] integer :: i,j integer :: ispin integer :: orb, hole_particle,spin_exc - double precision :: norm_out(N_states_diag) + 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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) use bitmasks integer :: iorb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) hole_particle = 1 spin_exc = ispin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -53,8 +53,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] enddo do state_target = 1,N_states call apply_exc_to_psi(orb,hole_particle,spin_exc, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -68,22 +68,22 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] integer :: i,j integer :: ispin integer :: orb, hole_particle,spin_exc - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) hole_particle = -1 spin_exc = ispin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -93,8 +93,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] enddo do state_target = 1, N_states call apply_exc_to_psi(orb,hole_particle,spin_exc, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -109,15 +109,15 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) integer :: ispin,jspin integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -129,7 +129,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) hole_particle_j = 1 spin_exc_j = jspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -139,10 +139,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) enddo do state_target = 1 , N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -159,16 +159,16 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) integer :: ispin,jspin integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: state_target state_target = 1 - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -180,7 +180,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) hole_particle_j = -1 spin_exc_j = jspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -189,10 +189,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) 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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -208,15 +208,15 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 integer :: ispin,jspin integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -228,7 +228,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 hole_particle_j = -1 spin_exc_j = jspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -238,14 +238,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) if(orb_i == orb_j .and. ispin .ne. jspin)then - call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) else - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) endif enddo @@ -264,16 +264,16 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j integer :: orb_k, hole_particle_k,spin_exc_k - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: korb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -290,7 +290,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a hole_particle_k = -1 spin_exc_k = kspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -301,12 +301,12 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a 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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -326,16 +326,16 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j integer :: orb_k, hole_particle_k,spin_exc_k - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: korb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -352,7 +352,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a hole_particle_k = -1 spin_exc_k = kspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -362,12 +362,12 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a enddo do state_target = 1, N_states call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -387,16 +387,16 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j integer :: orb_k, hole_particle_k,spin_exc_k - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: korb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -413,7 +413,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 hole_particle_k = 1 spin_exc_k = kspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -423,12 +423,12 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -448,16 +448,16 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 integer :: orb_i, hole_particle_i,spin_exc_i integer :: orb_j, hole_particle_j,spin_exc_j integer :: orb_k, hole_particle_k,spin_exc_k - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb integer :: korb integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) do iorb = 1,n_act_orb do ispin = 1,2 orb_i = list_act(iorb) @@ -474,7 +474,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 hole_particle_k = -1 spin_exc_k = kspin do i = 1, n_det - do j = 1, n_states_diag + do j = 1, n_states psi_in_out_coef(i,j) = psi_coef(i,j) enddo do j = 1, N_int @@ -484,12 +484,12 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,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 u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -511,15 +511,15 @@ END_PROVIDER integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb,i_ok integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) double precision :: hij double precision :: norm(N_states,2),norm_no_inv(N_states,2),norm_bis(N_states,2) double precision :: energies_alpha_beta(N_states,2) @@ -585,7 +585,7 @@ END_PROVIDER energies_alpha_beta(state_target, ispin) = - mo_bielec_integral_jj_exchange(orb_i,orb_v) ! energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo @@ -616,15 +616,15 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta integer :: i,iorb,j integer :: ispin,jspin integer :: orb_i, hole_particle_i - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: jorb,i_ok,aorb,orb_a integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) double precision :: hij double precision :: norm(N_states,2),norm_no_inv(N_states,2) double precision :: energies_alpha_beta(N_states,2) @@ -688,7 +688,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo @@ -701,11 +701,6 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta else one_anhil_inact(iorb,aorb,state_target) = 0.d0 endif -! print*, '********' -! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2) -! print*, norm_bis(state_target,1) , norm_bis(state_target,2) -! print*, one_anhil_inact(iorb,aorb,state_target) -! print*, one_creat(aorb,1,state_target) enddo enddo enddo @@ -719,15 +714,15 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v - double precision :: norm_out(N_states_diag) + double precision :: norm_out(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,N_states_diag)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) integer :: iorb,jorb,i_ok,aorb,orb_a integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) double precision :: hij double precision :: norm(N_states,2),norm_no_inv(N_states,2) double precision :: energies_alpha_beta(N_states,2) @@ -791,7 +786,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) ! print*, energies(state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif @@ -825,19 +820,19 @@ END_PROVIDER integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v - double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det) + double precision :: norm_out(N_states),diag_elem(N_det),interact_psi0(N_det) double precision :: delta_e_inact_virt(N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1)) allocate (eigenvectors(size(H_matrix,1),N_det+1)) allocate (eigenvalues(N_det+1)) integer :: iorb,jorb,i_ok integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) double precision :: hij double precision :: energies_alpha_beta(N_states,2) @@ -973,21 +968,21 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v - double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det) + double precision :: norm_out(N_states),diag_elem(N_det),interact_psi0(N_det) double precision :: delta_e_inact_virt(N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:) double precision, allocatable :: delta_e_det(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1)) allocate (eigenvectors(size(H_matrix,1),N_det+1)) allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det)) allocate (delta_e_det(N_det,N_det)) integer :: iorb,jorb,i_ok integer :: state_target - double precision :: energies(n_states_diag) + double precision :: energies(n_states) double precision :: hij double precision :: energies_alpha_beta(N_states,2) double precision :: lamda_pt2(N_det) @@ -1023,7 +1018,7 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from interact_psi0(i) = 0.d0 do j = 1 , N_det call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij) - call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j)) + call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),delta_e_det(i,j)) interact_cas(i,j) = hij interact_psi0(i) += hij * psi_coef(j,1) enddo diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 275af0e4..67a12af2 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -84,7 +84,8 @@ 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_coef(index_i,i_state) enddo - call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) + call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) +! call get_delta_e_dyall_general_mp(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) hij_array(index_i) = hialpha call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) ! phase_array(index_i) = phase diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index d7b1f0f6..b829dd40 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -58,8 +58,6 @@ delta_ij_tmp = 0.d0 call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) 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 @@ -121,7 +119,7 @@ ! 1h2p delta_ij_tmp = 0.d0 -!call give_1h2p_contrib(delta_ij_tmp) + call give_1h2p_contrib(delta_ij_tmp) call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states @@ -137,7 +135,7 @@ ! 2h1p delta_ij_tmp = 0.d0 -!call give_2h1p_contrib(delta_ij_tmp) + call give_2h1p_contrib(delta_ij_tmp) call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states @@ -223,9 +221,9 @@ END_PROVIDER enddo END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states_diag) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ] + BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states) ] + &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states) ] + &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states) ] BEGIN_DOC ! Eigenvectors/values of the CI matrix END_DOC @@ -244,14 +242,14 @@ END_PROVIDER double precision, allocatable :: e_array(:) integer, allocatable :: iorder(:) - ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors - do j=1,min(N_states_diag,N_det) + ! Guess values for the "N_states" states of the CI_dressed_pt2_new_eigenvectors + do j=1,min(N_states,N_det) do i=1,N_det CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j) enddo enddo - do j=N_det+1,N_states_diag + do j=N_det+1,N_states do i=1,N_det CI_dressed_pt2_new_eigenvectors(i,j) = 0.d0 enddo @@ -267,8 +265,8 @@ END_PROVIDER allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) allocate (eigenvalues(N_det)) call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) - CI_electronic_energy(:) = 0.d0 + Hmatrix_dressed_pt2_new_symmetrized(1,1,1),size(H_matrix_all_dets,1),N_det) + CI_electronic_dressed_pt2_new_energy(:) = 0.d0 if (s2_eig) then i_state = 0 allocate (s2_eigvalues(N_det)) @@ -291,54 +289,54 @@ END_PROVIDER ! Fill the first "i_state" states that have a correct S^2 value do j = 1, i_state do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) enddo - CI_electronic_energy(j) = eigenvalues(index_good_state_array(j)) - CI_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j)) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) enddo i_other_state = 0 do j = 1, N_det if(good_state_array(j))cycle i_other_state +=1 - if(i_state+i_other_state.gt.n_states_diag)then + if(i_state+i_other_state.gt.n_states)then exit endif do i=1,N_det - CI_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) enddo - CI_electronic_energy(i_state+i_other_state) = eigenvalues(j) - CI_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) enddo else print*,'' print*,'!!!!!!!! WARNING !!!!!!!!!' print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' + print*,' and the ',N_states,'states requested' print*,' We did not find any state with S^2 values close to ',expected_s2 print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_eigenvectors' + print*,' as the CI_dressed_pt2_new_eigenvectors' print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' print*,'' - do j=1,min(N_states_diag,N_det) + do j=1,min(N_states,N_det) do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo - CI_electronic_energy(j) = eigenvalues(j) - CI_eigenvectors_s2(j) = s2_eigvalues(j) + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) + CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(j) enddo endif deallocate(index_good_state_array,good_state_array) deallocate(s2_eigvalues) else - call u_0_S2_u_0(CI_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(eigenvectors,1)) - ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states_diag) + call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& + min(N_det,N_states),size(eigenvectors,1)) + ! Select the "N_states" states of lowest energy + do j=1,min(N_det,N_states) do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) + CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo - CI_electronic_energy(j) = eigenvalues(j) + CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) enddo endif deallocate(eigenvectors,eigenvalues) @@ -348,7 +346,7 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ] +BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states) ] implicit none BEGIN_DOC ! N_states lowest eigenvalues of the CI matrix @@ -357,11 +355,11 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag) ] integer :: j character*(8) :: st call write_time(output_determinants) - do j=1,N_states_diag + do j=1,N_states CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st)) - call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) + call write_double(output_determinants,CI_dressed_pt2_new_eigenvectors_s2(j),'S^2 of state '//trim(st)) enddo END_PROVIDER diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index fa5812e1..9388fe5b 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -804,7 +804,7 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) do jdet = 1,N_det double precision :: coef_array(N_states),hij_test call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono) - call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,coef_array,hij_test,delta_e) + call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,delta_e) do state_target = 1, N_states ! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1) matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target) diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 794742b4..981fb8c2 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -152,7 +152,7 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti end -subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) +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 ! @@ -170,7 +170,6 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) use bitmasks double precision, intent(out) :: delta_e_final(N_states) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) - double precision, intent(in) :: coef_array(N_states),hij integer :: i,j,k,l integer :: i_state @@ -293,20 +292,9 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) if (n_holes_act == 0 .and. n_particles_act == 1) then ispin = particle_list_practical(1,1) i_particle_act = particle_list_practical(2,1) -! call get_excitation_degree(det_1,det_2,degree,N_int) -! if(degree == 1)then -! call get_excitation(det_1,det_2,exc,degree,phase,N_int) -! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) -! i_hole = list_inact_reverse(h1) -! i_part = list_act_reverse(p1) -! do i_state = 1, N_states -! delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) -! enddo -! else if (degree == 2)then do i_state = 1, N_states delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) enddo -! endif else if (n_holes_act == 1 .and. n_particles_act == 0) then ispin = hole_list_practical(1,1) @@ -433,3 +421,159 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) end + + +subroutine get_delta_e_dyall_general_mp(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 as the sum of energies of excitations a la MP + ! + END_DOC + implicit none + use bitmasks + 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(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) + 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) + 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(N_states) + 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) + 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) + do i_state = 1, N_states + delta_e_virt += fock_virt_total_spin_trace(i_part_virt,i_state) + enddo + 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(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) + 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 ! spin + hole_list_practical(2,icountbis) = holes_active_list(i,1) ! index of active orb + 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 + + do i = 1, n_holes_act + ispin = hole_list_practical(1,i) + i_hole_act = hole_list_practical(2,i) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) + enddo + enddo + + do i = 1, n_particles_act + ispin = particle_list_practical(1,i) + i_particle_act = particle_list_practical(2,i) + do i_state = 1, N_states + delta_e_act(i_state) += one_creat(i_particle_act, ispin,i_state) + enddo + enddo + + 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/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index e4e94b7f..e0764d96 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -15,6 +15,57 @@ enddo END_PROVIDER + BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] + implicit none + integer :: i,j,ispin,istate + ispin = 1 + do istate = 1, N_states + do j = 1, mo_tot_num + do i = 1, mo_tot_num + one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_alpha(i,j,istate) + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do j = 1, mo_tot_num + do i = 1, mo_tot_num + one_body_dm_mo_spin_index(i,j,istate,ispin) = one_body_dm_mo_beta(i,j,istate) + enddo + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [ double precision, one_body_dm_dagger_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] + implicit none + integer :: i,j,ispin,istate + ispin = 1 + do istate = 1, N_states + do j = 1, mo_tot_num + one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_alpha(j,j,istate) + do i = j+1, mo_tot_num + one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate) + one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_alpha(i,j,istate) + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do j = 1, mo_tot_num + one_body_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_body_dm_mo_beta(j,j,istate) + do i = j+1, mo_tot_num + one_body_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_body_dm_mo_beta(i,j,istate) + one_body_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_body_dm_mo_beta(i,j,istate) + enddo + enddo + 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 @@ -285,6 +336,8 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha, (ao_num_align,ao_num) ] &BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,ao_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_ao_alpha_no_align, (ao_num,ao_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta_no_align, (ao_num,ao_num) ] BEGIN_DOC ! one body density matrix on the AO basis : rho_AO(alpha) , rho_AO(beta) END_DOC @@ -303,11 +356,16 @@ END_PROVIDER ! 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 - enddo enddo enddo enddo + do i = 1, ao_num + do j = 1, ao_num + one_body_dm_ao_alpha_no_align(j,i) = one_body_dm_ao_alpha(j,i) + one_body_dm_ao_beta_no_align(j,i) = one_body_dm_ao_beta(j,i) + enddo + enddo END_PROVIDER