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) print*, 'energy_cas_dyall(i)', energy_cas_dyall(i) enddo END_PROVIDER BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)] implicit none integer :: i double precision :: energies(N_states_diag) 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) energy_cas_dyall_no_exchange(i) = energies(i) print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i) enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] implicit none integer :: i,j integer :: ispin integer :: orb, hole_particle,spin_exc double precision :: norm_out(N_states_diag) 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)) use bitmasks integer :: iorb integer :: state_target double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) 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,2,i) = psi_active(j,2,i) enddo 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) one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] implicit none integer :: i,j integer :: ispin integer :: orb, hole_particle,spin_exc double precision :: norm_out(N_states_diag) 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)) integer :: iorb integer :: state_target double precision :: energies(n_states_diag) do iorb = 1,n_act_orb do ispin = 1,2 orb = list_act(iorb) 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,2,i) = psi_active(j,2,i) enddo 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) one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)] 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), 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)) integer :: iorb,jorb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)] 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), 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)) 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,2,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_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2,N_States)] implicit none integer :: i,j integer :: ispin,jspin 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), 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)) integer :: iorb,jorb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) 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) 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) 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) one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) endif enddo enddo enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin 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), 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)) integer :: iorb,jorb integer :: korb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin 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), 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)) integer :: iorb,jorb integer :: korb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin 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), 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)) integer :: iorb,jorb integer :: korb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)] implicit none integer :: i,j integer :: ispin,jspin,kspin 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), 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)) integer :: iorb,jorb integer :: korb integer :: state_target 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,2,i) = psi_active(j,2,i) enddo 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) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag) call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target) three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo enddo enddo enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt, (n_inact_orb,n_virt_orb,N_States)] &BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_norm, (n_inact_orb,n_virt_orb,N_States,2)] implicit none integer :: i,vorb,j integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v double precision :: norm_out(N_states_diag) 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)) integer :: iorb,jorb,i_ok integer :: state_target double precision :: energies(n_states_diag) 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) double precision :: thresh_norm thresh_norm = 1.d-10 do vorb = 1,n_virt_orb orb_v = list_virt(vorb) do iorb = 1, n_inact_orb orb_i = list_inact(iorb) norm = 0.d0 norm_bis = 0.d0 do ispin = 1,2 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 j = 1, N_int psi_in_out(j,1,i) = psi_det(j,1,i) psi_in_out(j,2,i) = psi_det(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok) if(i_ok.ne.1)then print*, orb_i,orb_v call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib coef = psi_coef(i,j) !* psi_coef(i,j) psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo enddo do j = 1, N_states if (dabs(norm(j,ispin)) .le. thresh_norm)then norm(j,ispin) = 0.d0 norm_no_inv(j,ispin) = norm(j,ispin) one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0 else norm_no_inv(j,ispin) = norm(j,ispin) one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin) norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) endif enddo do i = 1, N_det do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,2,i) = psi_active(j,2,i) enddo enddo do state_target = 1, N_states 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) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo enddo ! ispin do state_target = 1, N_states if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then ! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.5d0 * & ! ( energy_cas_dyall(state_target) - energies_alpha_beta(state_target,1) + & ! energy_cas_dyall(state_target) - energies_alpha_beta(state_target,2) ) ! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2) ! print*, norm_bis(state_target,1) , norm_bis(state_target,2) one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - & ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) else one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0 endif enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_States)] implicit none integer :: i,iorb,j integer :: ispin,jspin integer :: orb_i, hole_particle_i double precision :: norm_out(N_states_diag) 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)) integer :: jorb,i_ok,aorb,orb_a integer :: state_target double precision :: energies(n_states_diag) double precision :: hij double precision :: norm(N_states,2),norm_no_inv(N_states,2) double precision :: energies_alpha_beta(N_states,2) double precision :: norm_alpha_beta(N_states,2) double precision :: thresh_norm thresh_norm = 1.d-10 do aorb = 1,n_act_orb orb_a = list_act(aorb) do iorb = 1, n_inact_orb orb_i = list_inact(iorb) do state_target = 1, N_states one_anhil_inact(iorb,aorb,state_target) = 0.d0 enddo norm_alpha_beta = 0.d0 norm = 0.d0 norm_bis = 0.d0 do ispin = 1,2 do i = 1, n_det do j = 1, N_int psi_in_out(j,1,i) = psi_det(j,1,i) psi_in_out(j,2,i) = psi_det(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_a,ispin,i_ok) if(i_ok.ne.1)then do j = 1, n_states psi_in_out_coef(i,j) = 0.d0 enddo else call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib coef = psi_coef(i,j) !* psi_coef(i,j) psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo endif enddo do j = 1, N_states if (dabs(norm(j,ispin)) .le. thresh_norm)then norm(j,ispin) = 0.d0 norm_no_inv(j,ispin) = norm(j,ispin) else norm_no_inv(j,ispin) = norm(j,ispin) norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) endif enddo double precision :: norm_bis(N_states,2) do i = 1, N_det do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1)) psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1)) enddo enddo 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) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo enddo ! ispin do state_target = 1, N_states if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall(state_target) - & ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) 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 deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_States)] implicit none integer :: i,vorb,j integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v double precision :: norm_out(N_states_diag) 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)) integer :: iorb,jorb,i_ok,aorb,orb_a integer :: state_target double precision :: energies(n_states_diag) double precision :: hij double precision :: norm(N_states,2),norm_no_inv(N_states,2) double precision :: energies_alpha_beta(N_states,2) double precision :: norm_alpha_beta(N_states,2) double precision :: thresh_norm thresh_norm = 1.d-10 do aorb = 1,n_act_orb orb_a = list_act(aorb) do vorb = 1, n_virt_orb orb_v = list_virt(vorb) do state_target = 1, N_states one_creat_virt(aorb,vorb,state_target) = 0.d0 enddo norm_alpha_beta = 0.d0 norm = 0.d0 norm_bis = 0.d0 do ispin = 1,2 do i = 1, n_det do j = 1, N_int psi_in_out(j,1,i) = psi_det(j,1,i) psi_in_out(j,2,i) = psi_det(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_a,orb_v,ispin,i_ok) if(i_ok.ne.1)then do j = 1, n_states psi_in_out_coef(i,j) = 0.d0 enddo else call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib coef = psi_coef(i,j) !* psi_coef(i,j) psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo endif enddo do j = 1, N_states if (dabs(norm(j,ispin)) .le. thresh_norm)then norm(j,ispin) = 0.d0 norm_no_inv(j,ispin) = norm(j,ispin) else norm_no_inv(j,ispin) = norm(j,ispin) norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) endif enddo double precision :: norm_bis(N_states,2) do i = 1, N_det do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1)) psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1)) enddo enddo 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) ! print*, energies(state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo enddo ! ispin do state_target = 1, N_states if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - & ( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) & /( norm_bis(state_target,1) + norm_bis(state_target,2) ) else one_creat_virt(aorb,vorb,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_creat_virt(aorb,vorb,state_target) ! print*, one_anhil(aorb,1,state_target) enddo enddo enddo deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_bis, (n_inact_orb,n_virt_orb,N_det,N_States)] &BEGIN_PROVIDER [ double precision, corr_e_from_1h1p, (N_States)] implicit none integer :: i,vorb,j 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 :: 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 (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 :: hij double precision :: energies_alpha_beta(N_states,2) double precision :: accu(N_states),norm double precision :: amplitudes_alpha_beta(N_det,2) double precision :: delta_e_alpha_beta(N_det,2) corr_e_from_1h1p = 0.d0 do vorb = 1,n_virt_orb orb_v = list_virt(vorb) do iorb = 1, n_inact_orb orb_i = list_inact(iorb) ! print*, '---------------------------------------------------------------------------' do j = 1, N_states delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) & - fock_virt_total_spin_trace(orb_v,j) enddo do ispin = 1,2 do i = 1, n_det do j = 1, N_int psi_in_out(j,1,i) = psi_det(j,1,i) psi_in_out(j,2,i) = psi_det(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok) if(i_ok.ne.1)then print*, orb_i,orb_v call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif 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) interact_psi0(i) += hij * psi_coef(j,1) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,2,i) = psi_active(j,2,i) enddo call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij) diag_elem(i) = hij enddo do state_target = 1, N_states ! Building the Hamiltonian matrix H_matrix(1,1) = energy_cas_dyall(state_target) do i = 1, N_det ! interaction with psi0 H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target) H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target) ! diagonal elements H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target) ! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1) do j = i+1, N_det call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij) H_matrix(i+1,j+1) = hij !0.d0 ! H_matrix(j+1,i+1) = hij !0.d0 ! enddo enddo print*, '***' do i = 1, N_det+1 write(*,'(100(F16.10,X))')H_matrix(i,:) enddo call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1) corr_e_from_1h1p(state_target) += eigenvalues(1) - energy_cas_dyall(state_target) norm = 0.d0 do i = 1, N_det psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1) !! if(dabs(psi_coef(i,state_target)*) .gt. 1.d-8)then if(dabs(psi_in_out_coef(i,state_target)) .gt. 1.d-8)then ! if(dabs(interact_psi0(i)) .gt. 1.d-8)then delta_e_alpha_beta(i,ispin) = H_matrix(1,i+1) / psi_in_out_coef(i,state_target) ! delta_e_alpha_beta(i,ispin) = interact_psi0(i) / psi_in_out_coef(i,state_target) amplitudes_alpha_beta(i,ispin) = psi_in_out_coef(i,state_target) / psi_coef(i,state_target) else amplitudes_alpha_beta(i,ispin) = 0.d0 delta_e_alpha_beta(i,ispin) = delta_e_inact_virt(state_target) endif !! one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,ispin,state_target) = amplitudes_alpha_beta(i,ispin) norm += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target) enddo print*, 'Coef ' write(*,'(100(X,F16.10))')psi_coef(1:N_det,state_target) write(*,'(100(X,F16.10))')psi_in_out_coef(:,state_target) double precision :: coef_tmp(N_det) do i = 1, N_det coef_tmp(i) = psi_coef(i,1) * interact_psi0(i) / delta_e_alpha_beta(i,ispin) enddo write(*,'(100(X,F16.10))')coef_tmp(:) print*, 'naked interactions' write(*,'(100(X,F16.10))')interact_psi0(:) print*, '' print*, 'norm ',norm norm = 1.d0/(norm) accu(state_target) = 0.d0 do i = 1, N_det accu(state_target) += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target) * H_matrix(i+1,i+1) do j = i+1, N_det accu(state_target) += 2.d0 * psi_in_out_coef(i,state_target) * psi_in_out_coef(j,state_target) * H_matrix(i+1,j+1) enddo enddo accu(state_target) = accu(state_target) * norm print*, delta_e_inact_virt(state_target) print*, eigenvalues(1),accu(state_target),eigenvectors(1,1) print*, energy_cas_dyall(state_target) - accu(state_target), one_anhil_one_creat_inact_virt(iorb,vorb,state_target) + delta_e_inact_virt(state_target) enddo enddo ! ispin do state_target = 1, N_states do i = 1, N_det one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,state_target) = 0.5d0 * & ( delta_e_alpha_beta(i,1) + delta_e_alpha_beta(i,1)) enddo enddo print*, '***' write(*,'(100(X,F16.10))') write(*,'(100(X,F16.10))')delta_e_alpha_beta(:,2) ! write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,1,:) ! write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,2,:) print*, '---------------------------------------------------------------------------' enddo enddo deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues) print*, 'corr_e_from_1h1p,',corr_e_from_1h1p(:) END_PROVIDER subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles) implicit none double precision , intent(inout) :: matrix_1h1p(N_det,N_det,N_states) double precision , intent(out) :: e_corr_from_1h1p_singles(N_states) integer :: i,vorb,j 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 :: 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 (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 :: hij double precision :: energies_alpha_beta(N_states,2) double precision :: lamda_pt2(N_det) double precision :: accu(N_states),norm double precision :: amplitudes_alpha_beta(N_det,2) double precision :: delta_e_alpha_beta(N_det,2) double precision :: coef_array(N_states) double precision :: coef_perturb(N_det) double precision :: coef_perturb_bis(N_det) do vorb = 1,n_virt_orb orb_v = list_virt(vorb) do iorb = 1, n_inact_orb orb_i = list_inact(iorb) do j = 1, N_states delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) & - fock_virt_total_spin_trace(orb_v,j) enddo do ispin = 1,2 do i = 1, n_det do j = 1, N_int psi_in_out(j,1,i) = psi_det(j,1,i) psi_in_out(j,2,i) = psi_det(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok) if(i_ok.ne.1)then print*, orb_i,orb_v call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif 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)) interact_cas(i,j) = hij interact_psi0(i) += hij * psi_coef(j,1) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) psi_in_out(j,2,i) = psi_active(j,2,i) enddo call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij) diag_elem(i) = hij enddo do state_target = 1, N_states ! Building the Hamiltonian matrix H_matrix(1,1) = energy_cas_dyall(state_target) do i = 1, N_det ! interaction with psi0 H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target) H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target) ! diagonal elements H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target) ! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1) do j = i+1, N_det call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij) H_matrix(i+1,j+1) = hij !0.d0 ! H_matrix(j+1,i+1) = hij !0.d0 ! enddo enddo call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1) e_corr_from_1h1p_singles(state_target) += eigenvalues(1) - energy_cas_dyall(state_target) do i = 1, N_det psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1) coef_perturb(i) = 0.d0 do j = 1, N_det coef_perturb(i) += psi_coef(j,state_target) * interact_cas(i,j) *1.d0/delta_e_det(i,j) enddo coef_perturb_bis(i) = interact_psi0(i) / (eigenvalues(1) - H_matrix(i+1,i+1)) if(dabs(interact_psi0(i)) .gt. 1.d-12)then lamda_pt2(i) = psi_in_out_coef(i,state_target) / interact_psi0(i) else lamda_pt2(i) =energy_cas_dyall(state_target) - H_matrix(i+1,i+1) endif enddo if(dabs(eigenvalues(1) - energy_cas_dyall(state_target)).gt.1.d-10)then print*, '' do i = 1, N_det+1 write(*,'(100(F16.10))') H_matrix(i,:) enddo accu = 0.d0 do i = 1, N_det accu(state_target) += psi_in_out_coef(i,state_target) * interact_psi0(i) enddo print*, '' print*, 'e corr diagonal ',accu(state_target) accu = 0.d0 do i = 1, N_det accu(state_target) += coef_perturb(i) * interact_psi0(i) enddo print*, 'e corr perturb ',accu(state_target) accu = 0.d0 do i = 1, N_det accu(state_target) += coef_perturb_bis(i) * interact_psi0(i) enddo print*, 'e corr perturb EN',accu(state_target) print*, '' print*, 'coef diagonalized' write(*,'(100(F16.10,X))')psi_in_out_coef(:,state_target) print*, 'coef_perturb' write(*,'(100(F16.10,X))')coef_perturb(:) print*, 'coef_perturb EN' write(*,'(100(F16.10,X))')coef_perturb_bis(:) endif integer :: k do k = 1, N_det do i = 1, N_det matrix_1h1p(i,i,state_target) += interact_cas(k,i) * interact_cas(k,i) * lamda_pt2(k) do j = i+1, N_det matrix_1h1p(i,j,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k) matrix_1h1p(j,i,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k) enddo enddo enddo enddo enddo ! ispin enddo enddo deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues,interact_cas) deallocate(delta_e_det) end