diff --git a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f new file mode 100644 index 00000000..c65e89d8 --- /dev/null +++ b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f @@ -0,0 +1,60 @@ +program MRPT_Utils + implicit none + read_wf = .True. + touch read_wf +! call routine +! call routine_2 + call routine_3 +end + + +subroutine routine_3 + implicit none + provide one_creation + +end + +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) + 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) + 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 new file mode 100644 index 00000000..dda6e11e --- /dev/null +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -0,0 +1,457 @@ +BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] + implicit none + integer :: i + double precision :: energies(N_states_diag) + do i = 1, N_states + call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i) + energy_cas_dyall(i) = energies(i) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, one_creation, (n_act_orb,2)] + implicit none + integer :: i,j + integer :: ispin + integer :: orb, hole_particle,spin_exc + double precision :: norm_out(N_states_diag) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb + 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 + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_creation(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_anhilation, (n_act_orb,2)] + implicit none + integer :: i,j + integer :: ispin + integer :: orb, hole_particle,spin_exc + double precision :: norm_out(N_states_diag) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb + 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 + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_anhilation(iorb,ispin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_creation, (n_act_orb,n_act_orb,2,2)] + implicit none + integer :: i,j + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,jorb + integer :: state_target + state_target = 1 + double precision :: energies(n_states_diag) + do iorb = 1,n_act_orb + do ispin = 1,2 + orb_i = list_act(iorb) + hole_particle_i = 1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = 1 + spin_exc_j = jspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_anhilation, (n_act_orb,n_act_orb,2,2)] + implicit none + integer :: i,j + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,jorb + integer :: state_target + state_target = 1 + double precision :: energies(n_states_diag) + do iorb = 1,n_act_orb + do ispin = 1,2 + orb_i = list_act(iorb) + hole_particle_i = -1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = -1 + spin_exc_j = jspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_anhilation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_anhilation_one_creation, (n_act_orb,n_act_orb,2,2)] + implicit none + integer :: i,j + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,jorb + integer :: state_target + state_target = 1 + double precision :: energies(n_states_diag) + do iorb = 1,n_act_orb + do ispin = 1,2 + orb_i = list_act(iorb) + hole_particle_i = 1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = -1 + spin_exc_j = jspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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) + one_anhilation_one_creation(iorb,jorb,ispin,jspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_anhilation_one_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] + implicit none + integer :: i,j + integer :: ispin,jspin,kspin + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,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 + orb_i = list_act(iorb) + hole_particle_i = 1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = -1 + spin_exc_j = jspin + do korb = 1, n_act_orb + do kspin = 1,2 + orb_k = list_act(korb) + hole_particle_k = -1 + spin_exc_k = kspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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) + two_anhilation_one_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_creation_one_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] + implicit none + integer :: i,j + integer :: ispin,jspin,kspin + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,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 + orb_i = list_act(iorb) + hole_particle_i = 1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = 1 + spin_exc_j = jspin + do korb = 1, n_act_orb + do kspin = 1,2 + orb_k = list_act(korb) + hole_particle_k = -1 + spin_exc_k = kspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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) + two_creation_one_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, three_creation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] + implicit none + integer :: i,j + integer :: ispin,jspin,kspin + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,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 + orb_i = list_act(iorb) + hole_particle_i = 1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = 1 + spin_exc_j = jspin + do korb = 1, n_act_orb + do kspin = 1,2 + orb_k = list_act(korb) + hole_particle_k = 1 + spin_exc_k = kspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_creation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, three_anhilation, (n_act_orb,n_act_orb,n_act_orb,2,2,2)] + implicit none + integer :: i,j + integer :: ispin,jspin,kspin + 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) + integer(bit_kind) :: psi_in_out(N_int,2,n_det) + double precision :: psi_in_out_coef(n_det,N_states_diag) + use bitmasks + + integer :: iorb,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 + orb_i = list_act(iorb) + hole_particle_i = -1 + spin_exc_i = ispin + do jorb = 1, n_act_orb + do jspin = 1,2 + orb_j = list_act(jorb) + hole_particle_j = -1 + spin_exc_j = jspin + do korb = 1, n_act_orb + do kspin = 1,2 + orb_k = list_act(korb) + hole_particle_k = -1 + spin_exc_k = kspin + do i = 1, n_det + do j = 1, n_states_diag + psi_in_out_coef(i,j) = psi_coef(i,j) + enddo + do j = 1, N_int + psi_in_out(j,1,i) = psi_active(j,1,i) + psi_in_out(j,1,i) = psi_active(j,2,i) + 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_anhilation(iorb,jorb,korb,ispin,jspin,kspin) = energy_cas_dyall(state_target) - energies(state_target) + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f b/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f new file mode 100644 index 00000000..6ecffcce --- /dev/null +++ b/plugins/MRPT_Utils/energies_cas_spin_averaged.irp.f @@ -0,0 +1,190 @@ + +BEGIN_PROVIDER [ double precision, one_creation_spin_averaged, (n_act_orb)] + implicit none + integer :: i + do i = 1, n_act_orb + one_creation_spin_averaged(i) = one_creation(i,1) + one_creation(i,2) + one_creation_spin_averaged(i) = 0.5d0 * one_creation_spin_averaged(i) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, one_anhilation_spin_averaged, (n_act_orb)] + implicit none + integer :: i + do i = 1, n_act_orb + one_anhilation_spin_averaged(i) = one_anhilation(i,1) + one_anhilation(i,2) + one_anhilation_spin_averaged(i) = 0.5d0 * one_anhilation_spin_averaged(i) + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_creation_spin_averaged, (n_act_orb,n_act_orb)] + implicit none + integer :: i,j + integer :: ispin,jspin + double precision :: counting + do i = 1, n_act_orb + do j = 1, n_act_orb + two_creation_spin_averaged(j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + two_creation_spin_averaged(j,i) += two_creation(j,i,ispin,jspin) + counting += 1.d0 + enddo + enddo + two_creation_spin_averaged(j,i) = two_creation_spin_averaged(j,i) / counting + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_anhilation_spin_averaged, (n_act_orb,n_act_orb)] + implicit none + integer :: i,j + integer :: ispin,jspin + double precision :: counting + do i = 1, n_act_orb + do j = 1, n_act_orb + two_anhilation_spin_averaged(j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + two_anhilation_spin_averaged(j,i) += two_anhilation(j,i,ispin,jspin) + counting += 1.d0 + enddo + enddo + two_anhilation_spin_averaged(j,i) = two_anhilation_spin_averaged(j,i) / counting + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, one_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb)] + implicit none + integer :: i,j + integer :: ispin,jspin + double precision :: counting + do i = 1, n_act_orb + do j = 1, n_act_orb + one_anhilation_one_creation_spin_averaged(j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + one_anhilation_one_creation_spin_averaged(j,i) += one_anhilation_one_creation(j,i,jspin,ispin) + counting += 1.d0 + enddo + enddo + one_anhilation_one_creation_spin_averaged(j,i) = one_anhilation_one_creation_spin_averaged(j,i) / counting + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, two_anhilation_one_creation_spin_averaged, (n_act_orb,n_act_orb,n_act_orb)] + 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_anhilation_one_creation_spin_averaged(k,j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + do kspin = 1,2 + two_anhilation_one_creation_spin_averaged(k,j,i) += two_anhilation_one_creation(k,j,i,kspin,jspin,ispin) + counting += 1.d0 + enddo + enddo + two_anhilation_one_creation_spin_averaged(k,j,i) = two_anhilation_one_creation_spin_averaged(k,j,i) / counting + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, two_creation_one_anhilation_spin_averaged, (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_creation_one_anhilation_spin_averaged(k,j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + do kspin = 1,2 + two_creation_one_anhilation_spin_averaged(k,j,i) += two_creation_one_anhilation(k,j,i,kspin,jspin,ispin) + counting += 1.d0 + enddo + enddo + two_creation_one_anhilation_spin_averaged(k,j,i) = two_creation_one_anhilation_spin_averaged(k,j,i) / counting + enddo + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, three_creation_spin_averaged, (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_creation_spin_averaged(k,j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + do kspin = 1,2 + three_creation_spin_averaged(k,j,i) += three_creation(k,j,i,kspin,jspin,ispin) + counting += 1.d0 + enddo + enddo + three_creation_spin_averaged(k,j,i) = three_creation_spin_averaged(k,j,i) / counting + enddo + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, three_anhilation_spin_averaged, (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_anhilation_spin_averaged(k,j,i) = 0.d0 + counting = 0.d0 + do ispin = 1, 2 + do jspin = 1,2 + do kspin = 1,2 + three_anhilation_spin_averaged(k,j,i) += three_anhilation(k,j,i,kspin,jspin,ispin) + counting += 1.d0 + enddo + enddo + three_anhilation_spin_averaged(k,j,i) = three_anhilation_spin_averaged(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 new file mode 100644 index 00000000..c7e21d0c --- /dev/null +++ b/plugins/MRPT_Utils/fock_like_operators.irp.f @@ -0,0 +1,170 @@ + BEGIN_PROVIDER [double precision, fock_core_inactive, (mo_tot_num)] + BEGIN_DOC +! inactive part of the fock operator with contributions only from the inactive + END_DOC + implicit none + integer :: i,j + double precision :: accu + + integer :: j_inact_core_orb,i_inact_core_orb + do i = 1, n_core_inact_orb + accu = 0.d0 + i_inact_core_orb = list_core_inact(i) + do j = 1, n_core_inact_orb +! do j = 1, elec_alpha_num +! j_inact_core_orb = j + j_inact_core_orb = list_core_inact(j) + accu += 2.d0 * mo_bielec_integral_jj(i_inact_core_orb,j_inact_core_orb) & + - mo_bielec_integral_jj_exchange(i_inact_core_orb,j_inact_core_orb) + enddo + fock_core_inactive(i_inact_core_orb) = accu + mo_mono_elec_integral(i_inact_core_orb,i_inact_core_orb) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, fock_virt_from_core_inact, (mo_tot_num)] + BEGIN_DOC +! fock operator for the virtuals that comes from the doubly occupied orbitals + END_DOC + implicit none + integer :: i,j + double precision :: accu + + integer :: j_inact_core_orb,i_virt_orb + do i = 1, n_virt_orb + accu = 0.d0 + i_virt_orb = list_virt(i) + do j = 1, n_core_inact_orb +! do j = 1, elec_alpha_num +! j_inact_core_orb = j + j_inact_core_orb = list_core_inact(j) + accu += 2.d0 * mo_bielec_integral_jj(i_virt_orb,j_inact_core_orb) & + - mo_bielec_integral_jj_exchange(i_virt_orb,j_inact_core_orb) + enddo + fock_virt_from_core_inact(i_virt_orb) = accu + mo_mono_elec_integral(i_virt_orb,i_virt_orb) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, fock_core_inactive_from_act, (mo_tot_num,2)] + BEGIN_DOC +! inactive part of the fock operator with contributions only from the active + END_DOC + implicit none + integer :: i,j + double precision :: accu_coulomb,accu_exchange(2) + double precision :: na,nb,ntot + double precision :: coulomb, exchange + integer :: j_act_orb,i_inact_core_orb + + 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 + 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_DOC +! virtual part of the fock operator with contributions only from the active + END_DOC + implicit none + integer :: i,j + double precision :: accu_coulomb,accu_exchange(2) + double precision :: na,nb,ntot + double precision :: coulomb, exchange + integer :: j_act_orb,i_virt_orb + + 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 + 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_averaged, (mo_tot_num)] + 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_averaged(i_inact_core_orb) = 0.5d0 * (fock_core_inactive_total(i_inact_core_orb,1) + fock_core_inactive_total(i_inact_core_orb,2)) + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, fock_virt_total, (mo_tot_num,2)] +&BEGIN_PROVIDER [double precision, fock_virt_total_spin_averaged, (mo_tot_num)] + BEGIN_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) + fock_virt_total(i_virt_orb,2) = fock_virt_from_core_inact(i_virt_orb) + fock_virt_from_act(i_virt_orb,2) + fock_virt_total_spin_averaged(i_virt_orb) = 0.5d0 * ( fock_virt_total(i_virt_orb,1) + fock_virt_total(i_virt_orb,2) ) + enddo + END_PROVIDER + + + + + + BEGIN_PROVIDER [double precision, fock_operator_active_from_core_inact, (n_act_orb,n_act_orb)] + BEGIN_DOC +! active part of the fock operator with contributions only from the inactive + END_DOC + implicit none + integer :: i,j,k + double precision :: accu + double precision :: get_mo_bielec_integral,coulomb, exchange + PROVIDE mo_bielec_integrals_in_map + do i = 1, n_act_orb + do j = 1, n_act_orb + accu = 0.d0 + do k = 1, n_core_inact_orb + coulomb = get_mo_bielec_integral(k,i,k,j,mo_integrals_map) + exchange = get_mo_bielec_integral(k,i,i,k,mo_integrals_map) + accu += 2.d0 * coulomb - exchange + enddo + fock_operator_active_from_core_inact(i,j) = accu + enddo + enddo + + END_PROVIDER + + + + diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f new file mode 100644 index 00000000..e5848814 --- /dev/null +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -0,0 +1,279 @@ + + use bitmasks +BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] + BEGIN_DOC +! active part of psi + END_DOC + implicit none + use bitmasks + integer :: i,j,k,l + do i = 1, N_det + do j = 1, N_int + psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1)) + psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1)) + enddo + enddo +END_PROVIDER + + +subroutine give_holes_and_particles_in_active_space(det_1,det_2,n_holes_spin,n_particles_spin,n_holes,n_particles,& + holes_active_list,particles_active_list) + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_1(N_int,2) + integer(bit_kind),intent(in ) :: det_2(N_int,2) + integer, intent(out) :: n_holes_spin(2),n_particles_spin(2) + integer, intent(out) :: n_holes,n_particles + integer, intent(out) :: holes_active_list(2 * n_act_orb,2) + integer, intent(out) :: particles_active_list(2 * n_act_orb,2) + integer :: i + integer(bit_kind) :: holes(N_int,2) + integer(bit_kind) :: particles(N_int,2) + integer(bit_kind) :: det_tmp_2(N_int,2),det_tmp_1(N_int,2) + BEGIN_DOC +! returns the holes and particles operators WITHIN THE ACTIVE SPACE +! that connect det_1 and det_2. By definition, the holes/particles +! are such that one starts from det_1 and goes to det_2 +! +! n_holes is the total number of holes +! n_particles is the total number of particles +! n_holes_spin is the number of number of holes per spin (1=alpha, 2=beta) +! n_particles_spin is the number of number of particles per spin (1=alpha, 2=beta) +! holes_active_list is the index of the holes per spin, that ranges from 1 to n_act_orb +! particles_active_list is the index of the particles per spin, that ranges from 1 to n_act_orb + END_DOC + + call give_active_part_determinant(det_1,det_tmp_1) + call give_active_part_determinant(det_2,det_tmp_2) + do i = 1, N_int + holes(i,1) = iand(det_tmp_1(i,1),xor(det_tmp_1(i,1),det_tmp_2(i,1))) + holes(i,2) = iand(det_tmp_1(i,2),xor(det_tmp_1(i,2),det_tmp_2(i,2))) + particles(i,1) = iand(det_tmp_2(i,1),xor(det_tmp_1(i,1),det_tmp_2(i,1))) + particles(i,2) = iand(det_tmp_2(i,2),xor(det_tmp_1(i,2),det_tmp_2(i,2))) + enddo + + integer :: holes_list(N_int*bit_kind_size,2) + holes_list = 0 + call bitstring_to_list(holes(1,1), holes_list(1,1), n_holes_spin(1), N_int) + call bitstring_to_list(holes(1,2), holes_list(1,2), n_holes_spin(2), N_int) + + n_holes = 0 + do i = 1, n_holes_spin(1) + n_holes +=1 + holes_active_list(i,1) = list_act_reverse(holes_list(i,1)) + enddo + do i = 1, n_holes_spin(2) + n_holes +=1 + holes_active_list(i,2) = list_act_reverse(holes_list(i,2)) + enddo + + + integer :: particles_list(N_int*bit_kind_size,2) + particles_list = 0 + call bitstring_to_list(particles(1,1), particles_list(1,1), n_particles_spin(1), N_int) + call bitstring_to_list(particles(1,2), particles_list(1,2), n_particles_spin(2), N_int) + n_particles = 0 + do i = 1, n_particles_spin(1) + n_particles += 1 + particles_active_list(i,1) = list_act_reverse(particles_list(i,1)) + enddo + do i = 1, n_particles_spin(2) + n_particles += 1 + particles_active_list(i,2) = list_act_reverse(particles_list(i,2)) + enddo + +end + +subroutine give_holes_in_inactive_space(det_1,n_holes_spin,n_holes,holes_list) + BEGIN_DOC +! returns the holes operators WITHIN THE INACTIVE SPACE +! that has lead to det_1. +! +! n_holes is the total number of holes +! n_holes_spin is the number of number of holes per spin (1=alpha, 2=beta) +! holes_inactive_list is the index of the holes per spin, that ranges from 1 to mo_tot_num + END_DOC + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_1(N_int,2) + integer, intent(out) :: n_holes_spin(2) + integer, intent(out) :: n_holes + integer, intent(out) :: holes_list(N_int*bit_kind_size,2) + integer :: i + integer(bit_kind) :: holes(N_int,2) + integer(bit_kind) :: det_tmp_1(N_int,2) + + call give_core_inactive_part_determinant(det_1,det_tmp_1) + + do i = 1, N_int + holes(i,1) = iand(reunion_of_core_inact_bitmask(i,1),xor(det_tmp_1(i,1),reunion_of_core_inact_bitmask(i,1))) + holes(i,2) = iand(reunion_of_core_inact_bitmask(i,2),xor(det_tmp_1(i,2),reunion_of_core_inact_bitmask(i,2))) + enddo + holes_list = 0 + call bitstring_to_list(holes(1,1), holes_list(1,1), n_holes_spin(1), N_int) + call bitstring_to_list(holes(1,2), holes_list(1,2), n_holes_spin(2), N_int) + n_holes = n_holes_spin(1) + n_holes_spin(2) + +end + +subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,particles_list) + BEGIN_DOC +! returns the holes operators WITHIN THE VIRTUAL SPACE +! that has lead to det_1. +! +! n_particles is the total number of particles +! n_particles_spin is the number of number of particles per spin (1=alpha, 2=beta) +! particles_inactive_list is the index of the particles per spin, that ranges from 1 to mo_tot_num + END_DOC + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_1(N_int,2) + integer, intent(out) :: n_particles_spin(2) + integer, intent(out) :: n_particles + integer, intent(out) :: particles_list(N_int*bit_kind_size,2) + integer :: i + integer(bit_kind) :: det_tmp_1(N_int,2) + integer(bit_kind) :: particles(N_int,2) + + call give_virt_part_determinant(det_1,det_tmp_1) + + do i = 1, N_int + particles(i,1) = iand(virt_bitmask(i,1),xor(det_tmp_1(i,1),virt_bitmask(i,1))) + particles(i,2) = iand(virt_bitmask(i,2),xor(det_tmp_1(i,2),virt_bitmask(i,2))) + enddo + + particles_list = 0 + call bitstring_to_list(particles(1,1), particles_list(1,1), n_particles_spin(1), N_int) + call bitstring_to_list(particles(1,2), particles_list(1,2), n_particles_spin(2), N_int) + n_particles = n_particles_spin(1) + n_particles_spin(2) + + +end + +subroutine get_delta_e_dyall(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_averaged(i_hole_inact) + enddo + + do i = 1, n_holes_spin(2) + i_hole_inact = holes_list(i,2) + delta_e_inactive += fock_core_inactive_total_spin_averaged(i_hole_inact) + 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_averaged(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 + icount = 0 + do i = 1, n_holes_spin_act(1) + icount += 1 + holes_active_list_spin_traced(icount) = holes_active_list(i,1) + enddo + do i = 1, n_holes_spin_act(2) + icount += 1 + holes_active_list_spin_traced(icount) = holes_active_list(i,2) + enddo + if(icount .ne. n_holes) then + print * , 'pb in holes_active_list_spin_traced !!' + stop + endif + + icount = 0 + do i = 1, n_particles_spin_act(1) + icount += 1 + particles_active_list_spin_traced(icount) = particles_active_list(i,1) + enddo + do i = 1, n_particles_spin_act(2) + icount += 1 + particles_active_list_spin_traced(icount) = particles_active_list(i,2) + enddo + if(icount .ne. n_particles) then + 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 + + + if (n_holes_act == 1 .and. n_particles_act == 0) then + i_hole_act = holes_active_list_spin_traced(1) + delta_e_act += one_creation_spin_averaged(i_hole_act) + + else if (n_holes_act == 0 .and. n_particles_act == 1) then + i_particle_act = particles_active_list_spin_traced(1) + delta_e_act += one_anhilation_spin_averaged(i_particle_act) + + else if (n_holes_act == 1 .and. n_particles_act == 1) then + i_hole_act = holes_active_list_spin_traced(1) + i_particle_act = particles_active_list_spin_traced(1) + delta_e_act += one_anhilation_one_creation_spin_averaged(i_hole_act,i_particle_act) + + else if (n_holes_act == 2 .and. n_particles_act == 1) then + i_hole_act = holes_active_list_spin_traced(1) + j_hole_act = holes_active_list_spin_traced(2) + i_particle_act = particles_active_list_spin_traced(1) + delta_e_act += two_anhilation_one_creation_spin_averaged(i_hole_act,j_hole_act,i_particle_act) + + else if (n_holes_act == 1 .and. n_particles_act == 2) then + i_hole_act = holes_active_list_spin_traced(1) + i_particle_act = particles_active_list_spin_traced(1) + j_particle_act = particles_active_list_spin_traced(2) + delta_e_act += two_creation_one_anhilation_spin_averaged(i_hole_act,i_particle_act,j_particle_act) + + else if (n_holes_act == 3 .and. n_particles_act == 0) then + i_hole_act = holes_active_list_spin_traced(1) + j_hole_act = holes_active_list_spin_traced(2) + k_hole_act = holes_active_list_spin_traced(3) + delta_e_act += three_anhilation_spin_averaged(i_hole_act,j_hole_act,k_hole_act) + + else if (n_holes_act == 0 .and. n_particles_act == 3) then + i_particle_act = particles_active_list_spin_traced(1) + j_particle_act = particles_active_list_spin_traced(2) + k_particle_act = particles_active_list_spin_traced(3) + delta_e_act += three_creation_spin_averaged(i_particle_act,j_particle_act,k_particle_act) + + endif + + delta_e_final = delta_e_act + delta_e_virt + delta_e_inactive + +end diff --git a/plugins/MRPT_Utils/utils_bitmask.irp.f b/plugins/MRPT_Utils/utils_bitmask.irp.f new file mode 100644 index 00000000..1b262eb6 --- /dev/null +++ b/plugins/MRPT_Utils/utils_bitmask.irp.f @@ -0,0 +1,36 @@ + +subroutine give_active_part_determinant(det_in,det_out) + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_in(N_int,2) + integer(bit_kind),intent(out) :: det_out(N_int,2) + integer :: i + do i = 1,N_int + det_out(i,1) = iand(det_in(i,1),cas_bitmask(i,1,1)) + det_out(i,2) = iand(det_in(i,2),cas_bitmask(i,1,1)) + enddo +end + +subroutine give_core_inactive_part_determinant(det_in,det_out) + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_in(N_int,2) + integer(bit_kind),intent(out) :: det_out(N_int,2) + integer :: i + do i = 1,N_int + det_out(i,1) = iand(det_in(i,1),reunion_of_core_inact_bitmask(i,1)) + det_out(i,2) = iand(det_in(i,2),reunion_of_core_inact_bitmask(i,1)) + enddo +end + +subroutine give_virt_part_determinant(det_in,det_out) + implicit none + use bitmasks + integer(bit_kind),intent(in) :: det_in(N_int,2) + integer(bit_kind),intent(out) :: det_out(N_int,2) + integer :: i + do i = 1,N_int + det_out(i,1) = iand(det_in(i,1),virt_bitmask(i,1)) + det_out(i,2) = iand(det_in(i,2),virt_bitmask(i,1)) + enddo +end diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index f23d51e9..c52ed837 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -431,12 +431,40 @@ END_PROVIDER END_PROVIDER + + BEGIN_PROVIDER [ integer, list_core_inact, (n_core_inact_orb)] +&BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_tot_num)] + + implicit none + integer :: occ_inact(N_int*bit_kind_size) + integer :: itest,i + occ_inact = 0 + + call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), occ_inact(1), itest, N_int) + + list_core_inact_reverse = 0 + do i = 1, n_core_inact_orb + list_core_inact(i) = occ_inact(i) + list_core_inact_reverse(occ_inact(i)) = i + enddo + + END_PROVIDER + + BEGIN_PROVIDER [ integer, n_core_inact_orb ] + implicit none + integer :: i + n_core_inact_orb = 0 + do i = 1, N_int + n_core_inact_orb += popcnt(reunion_of_core_inact_bitmask(i,1)) + enddo + ENd_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] implicit none BEGIN_DOC ! Reunion of the core and inactive and virtual bitmasks END_DOC - integer :: i,j + integer :: i do i = 1, N_int reunion_of_core_inact_bitmask(i,1) = ior(core_bitmask(i,1),inact_bitmask(i,1)) reunion_of_core_inact_bitmask(i,2) = ior(core_bitmask(i,2),inact_bitmask(i,2)) @@ -474,6 +502,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), inact_virt_bitmask, (N_int,2)] +&BEGIN_PROVIDER [ integer(bit_kind), core_inact_virt_bitmask, (N_int,2)] implicit none BEGIN_DOC ! Reunion of the inactive and virtual bitmasks @@ -482,6 +511,8 @@ END_PROVIDER do i = 1, N_int inact_virt_bitmask(i,1) = ior(inact_bitmask(i,1),virt_bitmask(i,1)) inact_virt_bitmask(i,2) = ior(inact_bitmask(i,2),virt_bitmask(i,2)) + core_inact_virt_bitmask(i,1) = ior(core_bitmask(i,1),inact_virt_bitmask(i,1)) + core_inact_virt_bitmask(i,2) = ior(core_bitmask(i,2),inact_virt_bitmask(i,2)) enddo END_PROVIDER diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f index b7233beb..a487cc23 100644 --- a/src/Determinants/create_excitations.irp.f +++ b/src/Determinants/create_excitations.irp.f @@ -45,3 +45,16 @@ subroutine set_bit_to_integer(i_physical,key,Nint) j = i_physical-ishft(k-1,bit_kind_shift)-1 key(k) = ibset(key(k),j) end + + +subroutine clear_bit_to_integer(i_physical,key,Nint) + use bitmasks + implicit none + integer, intent(in) :: i_physical,Nint + integer(bit_kind), intent(inout) :: key(Nint) + integer :: k,j,i + k = ishft(i_physical-1,-bit_kind_shift)+1 + j = i_physical-ishft(k-1,bit_kind_shift)-1 + key(k) = ibclr(key(k),j) +end +