From 8a91b293bf60c9c2e0af24ce6afa353ff7c31383 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 17 Nov 2016 17:03:48 +0100 Subject: [PATCH] now eigenfunction of S^2 --- config/ifort.cfg | 2 +- plugins/MRPT_Utils/energies_cas.irp.f | 192 ++--------------- plugins/MRPT_Utils/excitations_cas.irp.f | 259 +++++++++++++++++++---- plugins/MRPT_Utils/new_way.irp.f | 120 ----------- plugins/MRPT_Utils/print_1h2p.irp.f | 12 +- plugins/MRPT_Utils/psi_active_prov.irp.f | 35 +-- plugins/MRPT_Utils/special_hij.irp.f | 183 ++++++++++++++++ src/Bitmask/bitmask_cas_routines.irp.f | 21 ++ 8 files changed, 477 insertions(+), 347 deletions(-) create mode 100644 plugins/MRPT_Utils/special_hij.irp.f diff --git a/config/ifort.cfg b/config/ifort.cfg index 4cf7829e..843e887b 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index c1ce50e7..f09c30cb 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -520,7 +520,7 @@ END_PROVIDER integer :: iorb,jorb,i_ok integer :: state_target double precision :: energies(n_states) - double precision :: hij + double precision :: hij,hij_test 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) @@ -552,11 +552,18 @@ END_PROVIDER call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif +! call i_H_j_no_k_operators_from_act(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij_test) call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) +! if(i==1.and.dabs(hij)>1.d-8)then +! if(dabs(hij)>1.d-8)then +! print*, ispin,vorb,iorb +! print*, i,hij,hij_test +! pause +! endif 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 + psi_in_out_coef(i,j) = coef * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo enddo @@ -582,22 +589,18 @@ END_PROVIDER 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 + 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,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) + 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) 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) - & +! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - & + one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(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 @@ -688,24 +691,20 @@ 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,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) 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) - & +! one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall(state_target) - & + one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall_no_exchange(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 @@ -735,7 +734,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State double precision :: thresh_norm - thresh_norm = 1.d-10 + thresh_norm = 1.d-20 do aorb = 1,n_act_orb orb_a = list_act(aorb) @@ -791,15 +790,16 @@ 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,state_target) -! print*, energies(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) + 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) 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) - & +! one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - & + one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall_no_exchange(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 @@ -815,154 +815,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State 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_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),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) - 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) diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f index 10cfe7c0..6028d4fa 100644 --- a/plugins/MRPT_Utils/excitations_cas.irp.f +++ b/plugins/MRPT_Utils/excitations_cas.irp.f @@ -212,52 +212,97 @@ double precision function diag_H_mat_elem_no_elec_check(det_in,Nint) core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb) - mo_bielec_integral_jj_exchange(jorb,iorb) enddo enddo -! print*,'core_act = ',core_act -! print*,'alpha_alpha = ',alpha_alpha -! print*,'alpha_beta = ',alpha_beta -! print*,'beta_beta = ',beta_beta -! print*,'mono_elec = ',mono_elec - -! do i = 1, n_core_inact_orb -! iorb = list_core_inact(i) -! diag_H_mat_elem_no_elec_check += 2.d0 * fock_core_inactive_total_spin_trace(iorb,1) -! enddo - - -!!!!!!!!!!!! -return -!!!!!!!!!!!! - - - ! alpha - alpha - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - diag_H_mat_elem_no_elec_check += 1.d0 * mo_mono_elec_integral(iorb,iorb) - do j = i+1, n_core_inact_orb - jorb = list_core_inact(j) - diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb) - 1.d0 * mo_bielec_integral_jj_exchange(jorb,iorb) - enddo - enddo - - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - diag_H_mat_elem_no_elec_check += 1.d0 * mo_mono_elec_integral(iorb,iorb) - do j = i+1, n_core_inact_orb - jorb = list_core_inact(j) - diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb) - 1.d0 * mo_bielec_integral_jj_exchange(jorb,iorb) - enddo - enddo - - do i = 1, n_core_inact_orb - iorb = list_core_inact(i) - do j = 1, n_core_inact_orb - jorb = list_core_inact(j) - diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb) - enddo - enddo end + + + +double precision function diag_H_mat_elem_no_elec_check_no_spin(det_in,Nint) + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + + integer :: i, j, iorb, jorb + integer :: occ(Nint*bit_kind_size,2) + integer :: elec_num_tab_local(2) + + double precision :: core_act + double precision :: alpha_alpha + double precision :: alpha_beta + double precision :: beta_beta + double precision :: mono_elec + core_act = 0.d0 + alpha_alpha = 0.d0 + alpha_beta = 0.d0 + beta_beta = 0.d0 + mono_elec = 0.d0 + + diag_H_mat_elem_no_elec_check_no_spin = 0.d0 + call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), N_int) + call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), N_int) + ! alpha - alpha + do i = 1, elec_num_tab_local(1) + iorb = occ(i,1) + diag_H_mat_elem_no_elec_check_no_spin += mo_mono_elec_integral(iorb,iorb) + mono_elec += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(1) + jorb = occ(j,1) + diag_H_mat_elem_no_elec_check_no_spin += mo_bielec_integral_jj(jorb,iorb) + alpha_alpha += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + + ! beta - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + diag_H_mat_elem_no_elec_check_no_spin += mo_mono_elec_integral(iorb,iorb) + mono_elec += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(2) + jorb = occ(j,2) + diag_H_mat_elem_no_elec_check_no_spin += mo_bielec_integral_jj(jorb,iorb) + beta_beta += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + + + ! alpha - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + do j = 1, elec_num_tab_local(1) + jorb = occ(j,1) + diag_H_mat_elem_no_elec_check_no_spin += mo_bielec_integral_jj(jorb,iorb) + alpha_beta += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + + + ! alpha - core-act + do i = 1, elec_num_tab_local(1) + iorb = occ(i,1) + do j = 1, n_core_inact_orb + jorb = list_core_inact(j) + diag_H_mat_elem_no_elec_check_no_spin += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + + ! beta - core-act + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + do j = 1, n_core_inact_orb + jorb = list_core_inact(j) + diag_H_mat_elem_no_elec_check_no_spin += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + +end + + subroutine i_H_j_dyall(key_i,key_j,Nint,hij) use bitmasks implicit none @@ -389,6 +434,133 @@ subroutine i_H_j_dyall(key_i,key_j,Nint,hij) end +subroutine i_H_j_dyall_no_spin(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem_no_elec_check, phase,phase_2 + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,1,2) .and. exc(1,1,2) == exc(1,2,1) )then + hij = 0.d0 + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, n_occ_ab(1) + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, n_occ_ab(2) + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, n_occ_ab(1) + hij = hij + mipi(occ(k,1)) !- miip(occ(k,1)) + enddo + do k = 1, n_occ_ab(2) + hij = hij + mipi(occ(k,2)) + enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, n_occ_ab(2) + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, n_occ_ab(1) + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, n_occ_ab(1) + hij = hij + mipi(occ(k,1)) + enddo + do k = 1, n_occ_ab(2) + hij = hij + mipi(occ(k,2)) !- miip(occ(k,2)) + enddo + + endif + hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) ) + + case (0) + double precision :: diag_H_mat_elem_no_elec_check_no_spin + hij = diag_H_mat_elem_no_elec_check_no_spin(key_i,Nint) + end select +end + + + subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target) use bitmasks implicit none @@ -414,6 +586,7 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe do j = 1, ndet if(psi_coef_tmp(j)==0.d0)cycle call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) +! call i_H_j_dyall_no_spin(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij enddo enddo diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index fa5812e1..3624b7d3 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -393,126 +393,6 @@ subroutine give_1h2p_contrib(matrix_1h2p) end -subroutine give_1h1p_contrib(matrix_1h1p) - use bitmasks - implicit none - double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*) - integer :: i,j,r,a,b - integer :: iorb, jorb, rorb, aorb, borb - integer :: ispin,jspin - integer :: idet,jdet - integer :: inint - integer :: elec_num_tab_local(2),acu_elec - integer(bit_kind) :: det_tmp(N_int,2) - integer :: exc(0:2,2,2) - integer :: accu_elec - double precision :: get_mo_bielec_integral - double precision :: active_int(n_act_orb,2) - double precision :: hij,phase - integer :: degree(N_det) - integer :: idx(0:N_det) - integer :: istate - double precision :: hja,delta_e_inact_virt(N_states) - integer :: kspin,degree_scalar -!matrix_1h1p = 0.d0 - - elec_num_tab_local = 0 - do inint = 1, N_int - elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) - elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) - enddo - do i = 1, n_inact_orb ! First inactive - iorb = list_inact(i) - do r = 1, n_virt_orb ! First virtual - rorb = list_virt(r) - do j = 1, N_states - delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(iorb,j) & - - fock_virt_total_spin_trace(rorb,j) - enddo - do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations - do jdet = 1, idx(0) - do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) - do inint = 1, N_int - det_tmp(inint,1) = psi_det(inint,1,idet) - det_tmp(inint,2) = psi_det(inint,2,idet) - enddo - ! Do the excitation inactive -- > virtual - double precision :: himono,delta_e(N_states),coef_mono(N_states) - call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin - call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin - call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono) - - do state_target = 1, N_states -! delta_e(state_target) = one_anhil_one_creat_inact_virt(i,r,state_target) + delta_e_inact_virt(state_target) - delta_e(state_target) = one_anhil_one_creat_inact_virt_bis(i,r,idet,state_target) - coef_mono(state_target) = himono / delta_e(state_target) - enddo - if(idx(jdet).ne.idet)then - call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) - if (exc(0,1,1) == 1) then - ! Mono alpha - aorb = (exc(1,2,1)) !!! a^{\dagger}_a - borb = (exc(1,1,1)) !!! a_{b} - jspin = 1 - else - ! Mono beta - aorb = (exc(1,2,2)) !!! a^{\dagger}_a - borb = (exc(1,1,2)) !!! a_{b} - jspin = 2 - endif - - call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int) - if(degree_scalar .ne. 2)then - print*, 'pb !!!' - print*, degree_scalar - call debug_det(psi_det(1,1,idx(jdet)),N_int) - call debug_det(det_tmp,N_int) - stop - endif - call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) - if(ispin == jspin )then - hij = -get_mo_bielec_integral(iorb,aorb,rorb,borb,mo_integrals_map) & - + get_mo_bielec_integral(iorb,aorb,borb,rorb,mo_integrals_map) - else - hij = get_mo_bielec_integral(iorb,borb,rorb,aorb,mo_integrals_map) - endif - hij = hij * phase - double precision :: hij_test - integer :: state_target - call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test) - if(dabs(hij - hij_test).gt.1.d-10)then - print*, 'ahah pb !!' - print*, 'hij .ne. hij_test' - print*, hij,hij_test - call debug_det(psi_det(1,1,idx(jdet)),N_int) - call debug_det(det_tmp,N_int) - print*, ispin, jspin - print*,iorb,borb,rorb,aorb - print*, phase - call i_H_j_verbose(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test) - stop - endif - do state_target = 1, N_states - matrix_1h1p(idx(jdet),idet,state_target) += hij* coef_mono(state_target) - enddo - else - do state_target = 1, N_states - matrix_1h1p(idet,idet,state_target) += himono * coef_mono(state_target) - enddo - endif - enddo - enddo - - - - enddo - enddo - enddo -end - subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) use bitmasks implicit none diff --git a/plugins/MRPT_Utils/print_1h2p.irp.f b/plugins/MRPT_Utils/print_1h2p.irp.f index d10e1fb5..03851d8a 100644 --- a/plugins/MRPT_Utils/print_1h2p.irp.f +++ b/plugins/MRPT_Utils/print_1h2p.irp.f @@ -2,7 +2,17 @@ program print_1h2p implicit none read_wf = .True. touch read_wf - call routine + call routine_2 +end + +subroutine routine_2 + implicit none + integer :: i,j + do i =1, n_inact_orb + write(*,'(100(F16.10,X))')one_anhil_one_creat_inact_virt(i,:,1) + enddo + + end subroutine routine diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 67501727..33cb5d5b 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -293,27 +293,38 @@ 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 + 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 + endif else if (n_holes_act == 1 .and. n_particles_act == 0) then ispin = hole_list_practical(1,1) i_hole_act = hole_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_act_reverse(h1) + i_part = list_virt_reverse(p1) + do i_state = 1, N_states + delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state) + enddo + else if (degree == 2)then do i_state = 1, N_states delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) enddo + endif else if (n_holes_act == 1 .and. n_particles_act == 1) then ! first hole @@ -415,7 +426,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) i_hole = list_inact_reverse(h1) i_part = list_virt_reverse(p1) do i_state = 1, N_states -! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) + delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) enddo endif else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then diff --git a/plugins/MRPT_Utils/special_hij.irp.f b/plugins/MRPT_Utils/special_hij.irp.f new file mode 100644 index 00000000..597d8ee3 --- /dev/null +++ b/plugins/MRPT_Utils/special_hij.irp.f @@ -0,0 +1,183 @@ + + +subroutine i_H_j_no_k_operators_from_act(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral, phase + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size), miip_other(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + has_mipi = .False. + logical :: is_i_in_active + double precision :: accu_a, accu_b, accu_core + accu_a = 0.d0 + accu_b = 0.d0 + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + if(.not.is_i_in_active(i))then + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + else +! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + accu_a += miip(i) + endif + enddo + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + if(.not.is_i_in_active(i))then + miip_other(i) = 0.d0 + else +! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + accu_b += miip(i) + endif + enddo +! print*, accu_a,accu_b,accu_a + accu_b + accu_a = 0.d0 + accu_b = 0.d0 + accu_core = mo_mono_elec_integral(m,p) + + do k = 1, elec_alpha_num + hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) + accu_Core += mipi(occ(k,1)) + if(is_i_in_active(occ(k,1)))then + accu_a += miip(occ(k,1)) + else + accu_Core -= miip(occ(k,1)) + endif + enddo +! print*, hij,accu_core + do k = 1, elec_beta_num + hij = hij + mipi(occ(k,2)) - miip_other(occ(k,2)) + accu_Core += mipi(occ(k,2)) + if(is_i_in_active(occ(k,2)))then + accu_b += miip_other(occ(k,2)) + else + accu_Core -= miip_other(occ(k,2)) + endif + enddo +! print*, hij,accu_core,accu_core - accu_a - accu_b +! print*, accu_a,accu_b,accu_a + accu_b +! pause + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + if(.not.is_i_in_active(i))then + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + else + miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + endif + enddo + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + if(.not.is_i_in_active(i))then + miip_other(i) = 0.d0 + else + miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + endif + enddo + do k = 1, elec_alpha_num + hij = hij + mipi(occ(k,1)) - miip_other(occ(k,1)) + enddo + do k = 1, elec_beta_num + hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hij = phase*(hij + mo_mono_elec_integral(m,p)) + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 87a02d10..5c170632 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -560,3 +560,24 @@ logical function is_i_in_virtual(i) endif end + +logical function is_i_in_active(i) + implicit none + integer,intent(in) :: i + integer(bit_kind) :: key(N_int) + integer :: k,j + integer :: accu + is_i_in_active = .False. + key= 0_bit_kind + k = ishft(i-1,-bit_kind_shift)+1 + j = i-ishft(k-1,bit_kind_shift)-1 + key(k) = ibset(key(k),j) + accu = 0 + do k = 1, N_int + accu += popcnt(iand(key(k),cas_bitmask(k,1,1))) + enddo + if(accu .ne. 0)then + is_i_in_active= .True. + endif + +end