diff --git a/config/ifort.cfg b/config/ifort.cfg index a738a83c..da414912 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/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index 021aa422..8a51c4fe 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -89,11 +89,11 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen end -subroutine is_a_good_candidate(threshold,is_ok,verbose) +subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop) use bitmasks implicit none double precision, intent(in) :: threshold - logical, intent(out) :: is_ok + logical, intent(out) :: is_ok,exit_loop logical, intent(in) :: verbose integer :: l,k,m @@ -111,7 +111,7 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose) enddo enddo !call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input) - call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose) + call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) if(do_it_perturbative)then if(is_ok)then N_det = N_det_generators @@ -135,14 +135,14 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose) end -subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose) +subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop) use bitmasks implicit none integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators) integer, intent(in) :: Ndet_generators double precision, intent(in) :: threshold logical, intent(in) :: verbose - logical, intent(out) :: is_ok + logical, intent(out) :: is_ok,exit_loop double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states) double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators) @@ -151,6 +151,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average logical :: is_a_ref_det(Ndet_generators) + exit_loop = .False. is_a_ref_det = .False. do i = 1, N_det_generators @@ -191,6 +192,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then is_ok = .False. + exit_loop = .True. return endif endif diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f index 8656b633..0b0902b0 100644 --- a/plugins/FOBOCI/fobo_scf.irp.f +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -1,6 +1,6 @@ program foboscf implicit none - call run_prepare +!call run_prepare no_oa_or_av_opt = .True. touch no_oa_or_av_opt call routine_fobo_scf diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index e81b3fc1..a072a918 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -38,6 +38,7 @@ subroutine FOBOCI_lmct_mlct_old_thr integer(bit_kind) , allocatable :: psi_singles(:,:,:) logical :: lmct double precision, allocatable :: psi_singles_coef(:,:) + logical :: exit_loop allocate( zero_bitmask(N_int,2) ) do i = 1, n_inact_orb lmct = .True. @@ -55,7 +56,7 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'Passed set generators' call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) - call is_a_good_candidate(threshold_lmct,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop) print*,'is_ok = ',is_ok if(.not.is_ok)cycle allocate(dressing_matrix(N_det_generators,N_det_generators)) @@ -161,9 +162,9 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'--------------------------' ! First set the current generators to the one of restart + call check_symetry(i_particl_osoci,thr,test_sym) call set_generators_to_generators_restart call set_psi_det_to_generators - call check_symetry(i_particl_osoci,thr,test_sym) if(.not.test_sym)cycle print*,'i_particl_osoci= ',i_particl_osoci ! Initialize the bitmask to the restart ones @@ -180,9 +181,15 @@ subroutine FOBOCI_lmct_mlct_old_thr call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) !! ! so all the mono excitation on the new generators - call is_a_good_candidate(threshold_mlct,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop) print*,'is_ok = ',is_ok - if(.not.is_ok)cycle + if(.not. is_ok)then + if(exit_loop)then + exit + else + cycle + endif + endif allocate(dressing_matrix(N_det_generators,N_det_generators)) if(.not.do_it_perturbative)then dressing_matrix = 0.d0 @@ -234,7 +241,7 @@ subroutine FOBOCI_mlct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - logical :: verbose,is_ok + logical :: verbose,is_ok,exit_loop verbose = .False. thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) @@ -274,7 +281,7 @@ subroutine FOBOCI_mlct_old call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) ! ! so all the mono excitation on the new generators - call is_a_good_candidate(threshold_mlct,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop) print*,'is_ok = ',is_ok is_ok =.True. if(.not.is_ok)cycle @@ -308,7 +315,7 @@ subroutine FOBOCI_lmct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - logical :: verbose,is_ok + logical :: verbose,is_ok,exit_loop verbose = .False. thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) @@ -346,7 +353,7 @@ subroutine FOBOCI_lmct_old call set_generators_to_psi_det call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) - call is_a_good_candidate(threshold_lmct,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop) print*,'is_ok = ',is_ok if(.not.is_ok)cycle ! ! so all the mono excitation on the new generators diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 4aca60d7..3ecd7977 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -212,8 +212,8 @@ subroutine update_density_matrix_osoci integer :: iorb,jorb do i = 1, mo_tot_num do j = 1, mo_tot_num - one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha(i,j) - one_body_dm_mo_alpha_generators_restart(i,j)) - one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta(i,j) - one_body_dm_mo_beta_generators_restart(i,j)) + one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j)) + one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j)) enddo enddo @@ -588,14 +588,14 @@ end integer :: i double precision :: accu_tot,accu_sd print*,'touched the one_body_dm_mo_beta' - one_body_dm_mo_alpha = one_body_dm_mo_alpha_osoci - one_body_dm_mo_beta = one_body_dm_mo_beta_osoci + one_body_dm_mo_alpha_average = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_average = one_body_dm_mo_beta_osoci touch one_body_dm_mo_alpha one_body_dm_mo_beta accu_tot = 0.d0 accu_sd = 0.d0 do i = 1, mo_tot_num - accu_tot += one_body_dm_mo_alpha(i,i) + one_body_dm_mo_beta(i,i) - accu_sd += one_body_dm_mo_alpha(i,i) - one_body_dm_mo_beta(i,i) + accu_tot += one_body_dm_mo_alpha_average(i,i) + one_body_dm_mo_beta_average(i,i) + accu_sd += one_body_dm_mo_alpha_average(i,i) - one_body_dm_mo_beta_average(i,i) enddo print*,'accu_tot = ',accu_tot print*,'accu_sdt = ',accu_sd diff --git a/plugins/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f index fb52a719..ac399ce7 100644 --- a/plugins/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -10,14 +10,28 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] 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) :: psi_in_out(N_int,2,N_det) - double precision :: psi_in_out_coef(N_det,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 @@ -45,6 +59,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -54,9 +69,10 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] 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) + 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 @@ -83,6 +99,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -93,9 +110,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) 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) + 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 @@ -131,6 +149,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -141,9 +160,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) 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) + 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 @@ -178,6 +198,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -188,10 +209,11 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 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) + 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) @@ -205,23 +227,28 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 orb_j = list_act(jorb) hole_particle_j = -1 spin_exc_j = jspin - do i = 1, n_det - do j = 1, n_states_diag - psi_in_out_coef(i,j) = psi_coef(i,j) + 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 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) + 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 - 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) - 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) - enddo enddo enddo enddo @@ -229,6 +256,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 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 @@ -237,9 +265,10 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a 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) + 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 @@ -286,6 +315,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -297,9 +327,10 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a 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) + 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 @@ -345,6 +376,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -356,9 +388,10 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 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) + 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 @@ -404,6 +437,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 enddo enddo enddo + deallocate(psi_in_out,psi_in_out_coef) END_PROVIDER @@ -415,9 +449,10 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 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) + 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 @@ -463,5 +498,617 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 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 diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f index 213abf7b..805070f7 100644 --- a/plugins/MRPT_Utils/excitations_cas.irp.f +++ b/plugins/MRPT_Utils/excitations_cas.irp.f @@ -382,8 +382,6 @@ subroutine i_H_j_dyall(key_i,key_j,Nint,hij) endif hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) ) -! hij = phase*(mo_mono_elec_integral(m,p) ) ! + fock_operator_active_from_core_inact(m,p) ) -! hij = 0.d0 case (0) hij = diag_H_mat_elem_no_elec_check(key_i,Nint) @@ -422,3 +420,289 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe energies(state_target) = accu deallocate(psi_coef_tmp) end + + +double precision function coulomb_value_no_check(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 + + coulomb_value_no_check = 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) + do j = i+1, elec_num_tab_local(1) + jorb = occ(j,1) + coulomb_value_no_check += mo_bielec_integral_jj_anti(jorb,iorb) + alpha_alpha += mo_bielec_integral_jj_anti(jorb,iorb) + enddo + enddo + + ! beta - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + do j = i+1, elec_num_tab_local(2) + jorb = occ(j,2) + coulomb_value_no_check += mo_bielec_integral_jj_anti(jorb,iorb) + beta_beta += mo_bielec_integral_jj_anti(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) + coulomb_value_no_check += mo_bielec_integral_jj(jorb,iorb) + alpha_beta += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + + +end + +subroutine i_H_j_dyall_no_exchange(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_schwartz + 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_no_exchange, phase,phase_2 + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(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,2,2) .and. exc(1,2,1) == exc(1,1,2))then + hij = 0.d0 + else + hij = phase*get_mo_bielec_integral_schwartz( & + 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_schwartz( & + 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_schwartz( & + 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_schwartz(m,i,p,i,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_schwartz(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)) + 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_schwartz(m,i,p,i,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_schwartz(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)) + enddo + + endif + hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) ) + + case (0) + hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint) + end select +end + + +double precision function diag_H_mat_elem_no_elec_check_no_exchange(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_exchange(2) + core_act_exchange = 0.d0 + diag_H_mat_elem_no_elec_check_no_exchange = 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_exchange += 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_exchange += 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_exchange += 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_exchange += 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_exchange += 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_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + core_act_exchange(1) += - mo_bielec_integral_jj_exchange(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_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb) + core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb) + enddo + enddo + +end + +subroutine u0_H_dyall_u0_no_exchange(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target) + use bitmasks + implicit none + integer, intent(in) :: N_states_in,ndet,dim_psi_in,dim_psi_coef,state_target + integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in) + double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states_in) + double precision, intent(out) :: energies(N_states_in) + + integer :: i,j + double precision :: hij,accu + energies = 0.d0 + accu = 0.d0 + double precision, allocatable :: psi_coef_tmp(:) + allocate(psi_coef_tmp(ndet)) + + do i = 1, ndet + psi_coef_tmp(i) = psi_in_coef(i,state_target) + enddo + + double precision :: hij_bis + do i = 1, ndet + if(psi_coef_tmp(i)==0.d0)cycle + do j = 1, ndet + if(psi_coef_tmp(j)==0.d0)cycle + call i_H_j_dyall_no_exchange(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 + energies(state_target) = accu + deallocate(psi_coef_tmp) +end diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 0c8ec98c..275af0e4 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -79,8 +79,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip phase_array =0.d0 do i = 1,idx_alpha(0) index_i = idx_alpha(i) - call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) + double precision :: coef_array(N_states) + do i_state = 1, N_states + coef_array(i_state) = psi_coef(index_i,i_state) + enddo + call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) hij_array(index_i) = hialpha call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) ! phase_array(index_i) = phase diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index e298ae67..c1da3670 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -57,6 +57,9 @@ ! 1h1p delta_ij_tmp = 0.d0 call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) + double precision :: e_corr_from_1h1p_singles(N_states) +!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) +!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) accu = 0.d0 do i_state = 1, N_states do i = 1, N_det @@ -69,6 +72,23 @@ enddo print*, '1h1p = ',accu + ! 1h1p third order + delta_ij_tmp = 0.d0 + call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) +!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) +!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) + accu = 0.d0 + do i_state = 1, N_states + do i = 1, N_det + do j = 1, N_det + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) + enddo + enddo + second_order_pt_new_1h1p(i_state) = accu(i_state) + enddo + print*, '1h1p(3)',accu + ! 2h delta_ij_tmp = 0.d0 call H_apply_mrpt_2h(delta_ij_tmp,N_det) @@ -101,6 +121,7 @@ ! 1h2p delta_ij_tmp = 0.d0 +!call give_1h2p_contrib(delta_ij_tmp) call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states @@ -116,6 +137,7 @@ ! 2h1p delta_ij_tmp = 0.d0 +!call give_2h1p_contrib(delta_ij_tmp) call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) accu = 0.d0 do i_state = 1, N_states @@ -159,7 +181,7 @@ accu = 0.d0 do i_state = 1, N_states do i = 1, N_det - write(*,'(1000(F16.10,x))')delta_ij(i,:,:) +! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) do j = i_state, N_det accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) enddo diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index 8cf6ec5f..09016ab0 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -391,3 +391,568 @@ 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_schwartz + 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_schwartz(iorb,aorb,rorb,borb,mo_integrals_map) & + + get_mo_bielec_integral_schwartz(iorb,aorb,borb,rorb,mo_integrals_map) + else + hij = get_mo_bielec_integral_schwartz(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 + double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*) + integer :: i,j,r,a,b + integer :: iorb, jorb, rorb, aorb, borb,s,sorb + integer :: ispin,jspin + integer :: idet,jdet + integer :: inint + integer :: elec_num_tab_local(2),acu_elec + integer(bit_kind) :: det_tmp(N_int,2),det_tmp_bis(N_int,2) + integer(bit_kind) :: det_pert(N_int,2,n_inact_orb,n_virt_orb,2) + double precision :: coef_det_pert(n_inact_orb,n_virt_orb,2,N_states,2) + double precision :: delta_e_det_pert(n_inact_orb,n_virt_orb,2,N_states) + double precision :: hij_det_pert(n_inact_orb,n_virt_orb,2,N_states) + integer :: exc(0:2,2,2) + integer :: accu_elec + double precision :: get_mo_bielec_integral_schwartz + 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 + double precision :: himono,delta_e(N_states),coef_mono(N_states) + integer :: state_target + do idet = 1, N_det + call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + 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 ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + do state_target = 1, N_states + coef_det_pert(i,r,ispin,state_target,1) = 0.d0 + coef_det_pert(i,r,ispin,state_target,2) = 0.d0 + enddo + 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 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 + 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 inint = 1, N_int + det_pert(inint,1,i,r,ispin) = det_tmp(inint,1) + det_pert(inint,2,i,r,ispin) = det_tmp(inint,2) + enddo + do state_target = 1, N_states + delta_e_det_pert(i,r,ispin,state_target) = one_anhil_one_creat_inact_virt(i,r,state_target) + delta_e_inact_virt(state_target) + coef_det_pert(i,r,ispin,state_target,1) = himono / delta_e_det_pert(i,r,ispin,state_target) + enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements + !!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations + enddo ! ispin + enddo ! rorb + enddo ! iorb + + 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 ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + do inint = 1, N_int + det_tmp(inint,1) = det_pert(inint,1,i,r,ispin) + det_tmp(inint,2) = det_pert(inint,2,i,r,ispin) + enddo + do j = 1, n_inact_orb ! First inactive + jorb = list_inact(j) + do s = 1, n_virt_orb ! First virtual + sorb = list_virt(s) + do jspin = 1, 2 ! spin of the couple a-a^dagger (i,r) + if(i==j.and.r==s.and.ispin==jspin)cycle + do inint = 1, N_int + det_tmp_bis(inint,1) = det_pert(inint,1,j,s,jspin) + det_tmp_bis(inint,2) = det_pert(inint,2,j,s,jspin) + enddo + call i_H_j(det_tmp_bis,det_tmp,N_int,himono) + do state_target = 1, N_states + coef_det_pert(i,r,ispin,state_target,2) += & + coef_det_pert(j,s,jspin,state_target,1) * himono / delta_e_det_pert(i,r,ispin,state_target) + enddo + enddo + enddo + enddo + enddo ! ispin + enddo ! rorb + enddo ! iorb + 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 ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + do state_target = 1, N_states + coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2) + enddo + + do inint = 1, N_int + det_tmp(inint,1) = det_pert(inint,1,i,r,ispin) + det_tmp(inint,2) = det_pert(inint,2,i,r,ispin) + enddo + do jdet = 1, idx(0) +! + 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 + 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) + double precision :: hij_test + hij_test = 0.d0 + call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test) + do state_target = 1, N_states + matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2) + enddo + else + hij_test = 0.d0 + call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij_test) + do state_target = 1, N_states + matrix_1h1p(idet,idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2) + enddo + endif + enddo + enddo + enddo + enddo + + enddo ! idet +end + + +subroutine give_1p_sec_order_singles_contrib(matrix_1p) + use bitmasks + implicit none + double precision , intent(inout) :: matrix_1p(N_det,N_det,*) + integer :: i,j,r,a,b + integer :: iorb, jorb, rorb, aorb, borb,s,sorb + integer :: ispin,jspin + integer :: idet,jdet + integer :: inint + integer :: elec_num_tab_local(2),acu_elec + integer(bit_kind) :: det_tmp(N_int,2),det_tmp_bis(N_int,2) + integer(bit_kind) :: det_pert(N_int,2,n_act_orb,n_virt_orb,2) + double precision :: coef_det_pert(n_act_orb,n_virt_orb,2,N_states,2) + double precision :: delta_e_det_pert(n_act_orb,n_virt_orb,2,N_states) + double precision :: hij_det_pert(n_act_orb,n_virt_orb,2) + integer :: exc(0:2,2,2) + integer :: accu_elec + double precision :: get_mo_bielec_integral_schwartz + double precision :: hij,phase + integer :: degree(N_det) + integer :: idx(0:N_det) + integer :: istate + double precision :: hja,delta_e_act_virt(N_states) + integer :: kspin,degree_scalar +!matrix_1p = 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 + double precision :: himono,delta_e(N_states),coef_mono(N_states) + integer :: state_target + do idet = 1, N_det + call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + do i = 1, n_act_orb ! First active + iorb = list_act(i) + do r = 1, n_virt_orb ! First virtual + rorb = list_virt(r) + do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + do state_target = 1, N_states + coef_det_pert(i,r,ispin,state_target,1) = 0.d0 + coef_det_pert(i,r,ispin,state_target,2) = 0.d0 + enddo + do j = 1, N_states + delta_e_act_virt(j) = - fock_virt_total_spin_trace(rorb,j) + enddo + 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 active -- > virtual + call do_mono_excitation(det_tmp,iorb,rorb,ispin,i_ok) + integer :: i_ok + if(i_ok .ne.1)then + do state_target = 1, N_states + coef_det_pert(i,r,ispin,state_target,1) = -1.d+10 + coef_det_pert(i,r,ispin,state_target,2) = -1.d+10 + hij_det_pert(i,r,ispin) = 0.d0 + enddo + do inint = 1, N_int + det_pert(inint,1,i,r,ispin) = 0_bit_kind + det_pert(inint,2,i,r,ispin) = 0_bit_kind + enddo + cycle + endif + call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono) + do inint = 1, N_int + det_pert(inint,1,i,r,ispin) = det_tmp(inint,1) + det_pert(inint,2,i,r,ispin) = det_tmp(inint,2) + enddo + do state_target = 1, N_states + delta_e_det_pert(i,r,ispin,state_target) = one_creat_virt(i,r,state_target) + delta_e_act_virt(state_target) + coef_det_pert(i,r,ispin,state_target,1) = himono / delta_e_det_pert(i,r,ispin,state_target) + hij_det_pert(i,r,ispin) = himono + enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements + !!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations + enddo ! ispin + enddo ! rorb + enddo ! iorb + +! do i = 1, n_act_orb ! First active +! do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) +! if(coef_det_pert(i,1,ispin,1,1) == -1.d+10)cycle +! iorb = list_act(i) +! do r = 1, n_virt_orb ! First virtual +! rorb = list_virt(r) +! do inint = 1, N_int +! det_tmp(inint,1) = det_pert(inint,1,i,r,ispin) +! det_tmp(inint,2) = det_pert(inint,2,i,r,ispin) +! enddo +! do j = 1, n_act_orb ! First active +! do jspin = 1, 2 ! spin of the couple a-a^dagger (i,r) +! if(coef_det_pert(j,1,jspin,1,1) == -1.d+10)cycle +! jorb = list_act(j) +! do s = 1, n_virt_orb ! First virtual +! sorb = list_virt(s) +! if(i==j.and.r==s.and.ispin==jspin)cycle +! do inint = 1, N_int +! det_tmp_bis(inint,1) = det_pert(inint,1,j,s,jspin) +! det_tmp_bis(inint,2) = det_pert(inint,2,j,s,jspin) +! enddo +! call i_H_j(det_tmp_bis,det_tmp,N_int,himono) +! do state_target = 1, N_states +! coef_det_pert(i,r,ispin,state_target,2) += & +! coef_det_pert(j,s,jspin,state_target,1) * himono / delta_e_det_pert(i,r,ispin,state_target) +! enddo +! enddo +! enddo +! enddo +! enddo ! ispin +! enddo ! rorb +! enddo ! iorb + + do i = 1, n_act_orb ! First active + do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + if(coef_det_pert(i,1,ispin,1,1) == -1.d+10)cycle + iorb = list_act(i) + do r = 1, n_virt_orb ! First virtual + rorb = list_virt(r) +! do state_target = 1, N_states +! coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2) +! enddo + do inint = 1, N_int + det_tmp(inint,1) = det_pert(inint,1,i,r,ispin) + det_tmp(inint,2) = det_pert(inint,2,i,r,ispin) + enddo + do jdet = 1,N_det + double precision :: coef_array(N_states),hij_test + call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono) + call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,coef_array,hij_test,delta_e) + do state_target = 1, N_states +! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1) + matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target) + enddo + enddo + enddo + enddo + enddo + + enddo ! idet +end + + + +subroutine give_1h1p_only_doubles_spin_cross(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_schwartz + 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(bit_kind) :: pert_det(N_int,2,n_act_orb,n_act_orb,2) + double precision :: pert_det_coef(n_act_orb,n_act_orb,2,N_states) + integer :: kspin,degree_scalar + integer :: other_spin(2) + other_spin(1) = 2 + other_spin(2) = 1 + double precision :: hidouble,delta_e(N_states) +!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_double_alpha_beta(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations + do ispin = 1, 2 + jspin = other_spin(ispin) + do a = 1, n_act_orb + aorb = list_act(a) + do b = 1, n_act_orb + borb = list_act(b) + 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 (i-->a)(ispin) + (b-->r)(other_spin(ispin)) + integer :: i_ok,corb,dorb + call do_mono_excitation(det_tmp,iorb,aorb,ispin,i_ok) + if(i_ok .ne. 1)then + do state_target = 1, N_states + pert_det_coef(a,b,ispin,state_target) = -100000.d0 + enddo + do inint = 1, N_int + pert_det(inint,1,a,b,ispin) = 0_bit_kind + pert_det(inint,2,a,b,ispin) = 0_bit_kind + enddo + cycle + endif + call do_mono_excitation(det_tmp,borb,rorb,jspin,i_ok) + if(i_ok .ne. 1)then + do state_target = 1, N_states + pert_det_coef(a,b,ispin,state_target) = -100000.d0 + enddo + do inint = 1, N_int + pert_det(inint,1,a,b,ispin) = 0_bit_kind + pert_det(inint,2,a,b,ispin) = 0_bit_kind + enddo + cycle + endif + do inint = 1, N_int + pert_det(inint,1,a,b,ispin) = det_tmp(inint,1) + pert_det(inint,2,a,b,ispin) = det_tmp(inint,2) + enddo + + call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hidouble) + do state_target = 1, N_states + delta_e(state_target) = one_anhil_one_creat(a,b,ispin,jspin,state_target) + delta_e_inact_virt(state_target) + pert_det_coef(a,b,ispin,state_target) = hidouble / delta_e(state_target) + matrix_1h1p(idet,idet,state_target) += hidouble * pert_det_coef(a,b,ispin,state_target) + enddo + enddo + enddo + enddo + do jdet = 1, idx(0) + if(idx(jdet).ne.idet)then + call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + integer :: c,d,state_target + integer(bit_kind) :: det_tmp_bis(N_int,2) + ! excitation from I --> J + ! (a->c) (alpha) + (b->d) (beta) + aorb = exc(1,1,1) + corb = exc(1,2,1) + c = list_act_reverse(corb) + borb = exc(1,1,2) + dorb = exc(1,2,2) + d = list_act_reverse(dorb) + ispin = 1 + jspin = 2 + do inint = 1, N_int + det_tmp(inint,1) = pert_det(inint,1,c,d,1) + det_tmp(inint,2) = pert_det(inint,2,c,d,1) + det_tmp_bis(inint,1) = pert_det(inint,1,c,d,2) + det_tmp_bis(inint,2) = pert_det(inint,2,c,d,2) + enddo + double precision :: hjdouble_1,hjdouble_2 + call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1) + call i_H_j(psi_det(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2) + do state_target = 1, N_states + matrix_1h1p(idx(jdet),idet,state_target) += (pert_det_coef(c,d,1,state_target) * hjdouble_1 + pert_det_coef(c,d,2,state_target) * hjdouble_2 ) + enddo + endif + enddo + + + + enddo + enddo + enddo + + + + + +end + + diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 5a60d093..b4c7e6f4 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -14,8 +14,6 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] 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 - -! call debug_det(psi_active(1,1,i),N_int) enddo END_PROVIDER @@ -154,7 +152,7 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti end -subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) +subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final) BEGIN_DOC ! routine that returns the delta_e with the Moller Plesset and Dyall operators ! @@ -172,6 +170,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) use bitmasks double precision, intent(out) :: delta_e_final(N_states) integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) + double precision, intent(in) :: coef_array(N_states),hij integer :: i,j,k,l integer :: i_state @@ -292,23 +291,52 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) ! delta_e_act += one_creat_spin_trace(i_particle_act ) ispin = particle_list_practical(1,1) i_particle_act = particle_list_practical(2,1) - do i_state = 1, N_states - delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) - enddo + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 1)then + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + i_hole = list_inact_reverse(h1) + i_part = list_act_reverse(p1) + do i_state = 1, N_states + delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state) + enddo + else if (degree == 2)then + do i_state = 1, N_states + delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state) + enddo + endif else if (n_holes_act == 1 .and. n_particles_act == 0) then ! i_hole_act = holes_active_list_spin_traced(1) ! delta_e_act += one_anhil_spin_trace(i_hole_act ) ispin = hole_list_practical(1,1) i_hole_act = hole_list_practical(2,1) - do i_state = 1, N_states - delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state) - enddo + 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 + if(isnan(one_creat_virt(i_hole,i_part,i_state)))then + print*, i_hole,i_part,i_state + call debug_det(det_1,N_int) + call debug_det(det_2,N_int) + stop + endif + 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 ! i_hole_act = holes_active_list_spin_traced(1) ! i_particle_act = particles_active_list_spin_traced(1) ! delta_e_act += one_anhil_one_creat_spin_trace(i_hole_act,i_particle_act) + ! first hole ispin = hole_list_practical(1,1) i_hole_act = hole_list_practical(2,1) @@ -422,6 +450,34 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) do i_state = 1, N_states delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state) enddo + + else if (n_holes_act .eq. 0 .and. n_particles_act .eq.0)then + integer :: degree + integer(bit_kind) :: det_1_active(N_int,2) + integer :: h1,h2,p1,p2,s1,s2 + integer :: exc(0:2,2,2) + integer :: i_hole, i_part + double precision :: phase + call get_excitation_degree(det_1,det_2,degree,N_int) + if(degree == 1)then +! call debug_det(det_1,N_int) + call get_excitation(det_1,det_2,exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + i_hole = list_inact_reverse(h1) + i_part = list_virt_reverse(p1) + do i_state = 1, N_states +! if(one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1).gt.1.d-10)then +! print*, hij, one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) +! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) & +! * coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) +! print*, coef_array(i_state)* hij*coef_array(i_state)* hij,one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1), & +! coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1) +! else + delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) +! endif + enddo + endif + else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then @@ -438,3 +494,4 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) enddo end + diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 2621a9c6..c284e01e 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -192,6 +192,11 @@ subroutine pt2_moller_plesset ($arguments) endif do i =1,N_st H_pert_diag(i) = h +! if(dabs(i_H_psi_array(i)).gt.1.d-8)then +! print*, i_H_psi_array(i) +! call debug_det(det_pert,N_int) +! print*, h1,p1,h2,p2,s1,s2 +! endif c_pert(i) = i_H_psi_array(i) *delta_e e_2_pert(i) = c_pert(i) * i_H_psi_array(i) enddo diff --git a/plugins/Properties/hyperfine_constants.irp.f b/plugins/Properties/hyperfine_constants.irp.f index e31b3ba4..de14da9f 100644 --- a/plugins/Properties/hyperfine_constants.irp.f +++ b/plugins/Properties/hyperfine_constants.irp.f @@ -102,6 +102,11 @@ END_PROVIDER conversion_factor_gauss_hcc(3) = 619.9027742370165d0 conversion_factor_cm_1_hcc(3) = 579.4924475562677d0 + ! boron + conversion_factor_mhz_hcc(5) = 1434.3655101868d0 + conversion_factor_gauss_hcc(5) = 511.817264334d0 + conversion_factor_cm_1_hcc(5) = 478.4528336953d0 + ! carbon conversion_factor_mhz_hcc(6) = 1124.18303629792945d0 conversion_factor_gauss_hcc(6) = 401.136570647523058d0 diff --git a/plugins/Selectors_full/zmq.irp.f b/plugins/Selectors_full/zmq.irp.f index 952e5c06..9f6f616c 100644 --- a/plugins/Selectors_full/zmq.irp.f +++ b/plugins/Selectors_full/zmq.irp.f @@ -98,7 +98,14 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id) if (N_det_selectors_read > 0) then N_det_selectors = N_det_selectors_read endif - SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators + SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators psi_coef_generators psi_det_generators +! n_det_generators +! n_det_selectors +! psi_coef +! psi_coef_generators +! psi_det +! psi_det_generators + end diff --git a/src/Bitmask/bitmask_cas_routines.irp.f b/src/Bitmask/bitmask_cas_routines.irp.f index 5cd09aa2..961be5df 100644 --- a/src/Bitmask/bitmask_cas_routines.irp.f +++ b/src/Bitmask/bitmask_cas_routines.irp.f @@ -541,3 +541,24 @@ use bitmasks end +logical function is_i_in_virtual(i) +use bitmasks + implicit none + integer,intent(in) :: i + integer(bit_kind) :: key(N_int) + integer :: k,j + integer :: accu + is_i_in_virtual = .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),virt_bitmask(k,1))) + enddo + if(accu .ne. 0)then + is_i_in_virtual = .True. + endif + +end diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index c52ed837..964c4ed8 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -37,7 +37,7 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask_4, (N_int,4) ] enddo END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ] + BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ] implicit none integer :: i do i=1,N_int @@ -473,17 +473,33 @@ END_PROVIDER BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] +&BEGIN_PROVIDER [ integer, n_core_inact_act_orb ] implicit none BEGIN_DOC ! Reunion of the core, inactive and active bitmasks END_DOC integer :: i,j + n_core_inact_act_orb = 0 do i = 1, N_int reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1)) reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1)) + n_core_inact_act_orb +=popcnt(reunion_of_core_inact_act_bitmask(i,1)) enddo END_PROVIDER + BEGIN_PROVIDER [ integer, list_core_inact_act, (n_core_inact_act_orb)] +&BEGIN_PROVIDER [ integer, list_core_inact_act_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_act_bitmask(1,1), occ_inact(1), itest, N_int) + list_inact_reverse = 0 + do i = 1, n_core_inact_act_orb + list_core_inact_act(i) = occ_inact(i) + list_core_inact_act_reverse(occ_inact(i)) = i + enddo +END_PROVIDER @@ -543,8 +559,8 @@ END_PROVIDER integer :: i,j n_core_orb = 0 do i = 1, N_int - core_bitmask(i,1) = xor(closed_shell_ref_bitmask(i,1),reunion_of_cas_inact_bitmask(i,1)) - core_bitmask(i,2) = xor(closed_shell_ref_bitmask(i,2),reunion_of_cas_inact_bitmask(i,2)) + core_bitmask(i,1) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,1),virt_bitmask(i,1))) + core_bitmask(i,2) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1))) n_core_orb += popcnt(core_bitmask(i,1)) enddo print*,'n_core_orb = ',n_core_orb diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f index 6af49681..71301dbc 100644 --- a/src/Determinants/create_excitations.irp.f +++ b/src/Determinants/create_excitations.irp.f @@ -37,6 +37,72 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) endif end +subroutine do_spin_flip(key_in,i_flip,ispin,i_ok) + implicit none + BEGIN_DOC + ! flip the spin ispin in the orbital i_flip + ! on key_in + ! ispin = 1 == alpha + ! ispin = 2 == beta + ! i_ok = 1 == the flip is possible + ! i_ok = -1 == the flip is not possible + END_DOC + integer, intent(in) :: i_flip,ispin + integer(bit_kind), intent(inout) :: key_in(N_int,2) + integer, intent(out) :: i_ok + integer :: k,j,i + integer(bit_kind) :: key_tmp(N_int,2) + i_ok = -1 + key_tmp = 0_bit_kind + k = ishft(i_flip-1,-bit_kind_shift)+1 + j = i_flip-ishft(k-1,bit_kind_shift)-1 + key_tmp(k,1) = ibset(key_tmp(k,1),j) + integer :: other_spin(2) + other_spin(1) = 2 + other_spin(2) = 1 + if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then + ! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip" + key_in(k,ispin) = ibclr(key_in(k,ispin),j) ! destroy the electron ispin in the orbital i_flip + key_in(k,other_spin(ispin)) = ibset(key_in(k,other_spin(ispin)),j) ! create an electron of spin other_spin in the same orbital + i_ok = 1 + else + return + endif + + + +end + +logical function is_spin_flip_possible(key_in,i_flip,ispin) + implicit none + BEGIN_DOC + ! returns .True. if the spin-flip of spin ispin in the orbital i_flip is possible + ! on key_in + END_DOC + integer, intent(in) :: i_flip,ispin + integer(bit_kind), intent(in) :: key_in(N_int,2) + integer :: k,j,i + integer(bit_kind) :: key_tmp(N_int,2) + is_spin_flip_possible = .False. + key_tmp = 0_bit_kind + k = ishft(i_flip-1,-bit_kind_shift)+1 + j = i_flip-ishft(k-1,bit_kind_shift)-1 + key_tmp(k,1) = ibset(key_tmp(k,1),j) + integer :: other_spin(2) + other_spin(1) = 2 + other_spin(2) = 1 + if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then + ! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip" + is_spin_flip_possible = .True. + return + else + return + endif + + + +end + subroutine set_bit_to_integer(i_physical,key,Nint) use bitmasks implicit none diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index 9810b219..344e0160 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -334,7 +334,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nma accu_precision_of_diag = 0.d0 do i = 1, nstates do j = i+1, nstates - if( ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then + if( ( dabs(s2(i,i) - s2(j,j)) .le.0.5d0 ) ) then s2(i,j) = 0.d0 s2(j,i) = 0.d0 endif diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index cb96957a..7f0e7e57 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1389,6 +1389,55 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s l = l+1 endif enddo + else + + print*, 'get_excitation_degree_vector_mono_or_exchange not yet implemented for N_int > 1 ...' + stop + + endif + idx(0) = l-1 +end + + + + +subroutine get_excitation_degree_vector_double_alpha_beta(key1,key2,degree,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Applies get_excitation_degree to an array of determinants and return only the mono excitations + ! and the connections through exchange integrals + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: degree(sze) + integer, intent(out) :: idx(0:sze) + integer(bit_kind) :: key_tmp(Nint,2) + + integer :: i,l,d,m + integer :: degree_alpha, degree_beta + + ASSERT (Nint > 0) + ASSERT (sze > 0) + + l=1 + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (d .ne.4)cycle + key_tmp(1,1) = xor(key1(1,1,i),key2(1,1)) + key_tmp(1,2) = xor(key1(1,2,i),key2(1,2)) + degree_alpha = popcnt(key_tmp(1,1)) + degree_beta = popcnt(key_tmp(1,2)) + if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + enddo else if (Nint==2) then !DIR$ LOOP COUNT (1000) @@ -1397,29 +1446,17 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s popcnt(xor( key1(1,2,i), key2(1,2))) + & popcnt(xor( key1(2,1,i), key2(2,1))) + & popcnt(xor( key1(2,2,i), key2(2,2))) - exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,2),key2(1,2)))) + & - popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,2),key2(2,2)))) - exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + & - popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) - if (d > 4)cycle - if (d ==4)then - if(exchange_1 .eq. 0 ) then + if (d .ne.4)cycle + key_tmp(1,1) = xor(key1(1,1,i),key2(1,1)) + key_tmp(1,2) = xor(key1(1,2,i),key2(1,2)) + key_tmp(2,1) = xor(key1(2,1,i),key2(2,1)) + key_tmp(2,2) = xor(key1(2,2,i),key2(2,2)) + degree_alpha = popcnt(key_tmp(1,1)) + popcnt(key_tmp(2,1)) + degree_beta = popcnt(key_tmp(1,2)) + popcnt(key_tmp(2,2)) + if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin degree(l) = ishft(d,-1) idx(l) = i l = l+1 - else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - else - cycle - endif -! pause - else - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - endif enddo else if (Nint==3) then @@ -1432,31 +1469,19 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s popcnt(xor( key1(2,2,i), key2(2,2))) + & popcnt(xor( key1(3,1,i), key2(3,1))) + & popcnt(xor( key1(3,2,i), key2(3,2))) - exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + & - popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,1),key2(2,2)))) + & - popcnt(xor(ior(key1(3,1,i),key1(3,2,i)),ior(key2(3,1),key2(3,2)))) - exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + & - popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) + & - popcnt(ior(xor(key1(3,1,i),key2(3,1)),xor(key1(3,2,i),key2(3,2)))) - if (d > 4)cycle - if (d ==4)then - if(exchange_1 .eq. 0 ) then + if (d .ne.4)cycle + key_tmp(1,1) = xor(key1(1,1,i),key2(1,1)) + key_tmp(1,2) = xor(key1(1,2,i),key2(1,2)) + key_tmp(2,1) = xor(key1(2,1,i),key2(2,1)) + key_tmp(2,2) = xor(key1(2,2,i),key2(2,2)) + key_tmp(3,1) = xor(key1(3,1,i),key2(3,1)) + key_tmp(3,2) = xor(key1(3,2,i),key2(3,2)) + degree_alpha = popcnt(key_tmp(1,1)) + popcnt(key_tmp(2,1)) + popcnt(key_tmp(3,1)) + degree_beta = popcnt(key_tmp(1,2)) + popcnt(key_tmp(2,2)) + popcnt(key_tmp(3,2)) + if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin degree(l) = ishft(d,-1) idx(l) = i l = l+1 - else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - else - cycle - endif -! pause - else - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - endif enddo else @@ -1464,39 +1489,26 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s !DIR$ LOOP COUNT (1000) do i=1,sze d = 0 - exchange_1 = 0 !DIR$ LOOP COUNT MIN(4) do m=1,Nint d = d + popcnt(xor( key1(m,1,i), key2(m,1))) & + popcnt(xor( key1(m,2,i), key2(m,2))) - exchange_1 = popcnt(xor(ior(key1(m,1,i),key1(m,2,i)),ior(key2(m,1),key2(m,2)))) - exchange_2 = popcnt(ior(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2)))) + key_tmp(m,1) = xor(key1(m,1,i),key2(m,1)) + key_tmp(m,2) = xor(key1(m,2,i),key2(m,2)) + degree_alpha = popcnt(key_tmp(m,1)) + degree_beta = popcnt(key_tmp(m,2)) enddo - if (d > 4)cycle - if (d ==4)then - if(exchange_1 .eq. 0 ) then + if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin degree(l) = ishft(d,-1) idx(l) = i l = l+1 - else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - else - cycle - endif -! pause - else - degree(l) = ishft(d,-1) - idx(l) = i - l = l+1 - endif enddo endif idx(0) = l-1 end + subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degree,Nint,sze,idx) use bitmasks implicit none diff --git a/src/Integrals_Bielec/EZFIO.cfg b/src/Integrals_Bielec/EZFIO.cfg index feed02c1..04d58236 100644 --- a/src/Integrals_Bielec/EZFIO.cfg +++ b/src/Integrals_Bielec/EZFIO.cfg @@ -10,7 +10,7 @@ type: logical doc: If True, do not compute the bielectronic integrals when 4 indices are virtual interface: ezfio,provider,ocaml default: False -ezfio_name: None +ezfio_name: no_vvvv_integrals [disk_access_mo_integrals] type: Disk_access diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 3557772d..dae73a01 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -22,6 +22,7 @@ end BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] implicit none integer(bit_kind) :: mask_ijkl(N_int,4) + integer(bit_kind) :: mask_ijk(N_int,3) BEGIN_DOC ! If True, the map of MO bielectronic integrals is provided @@ -39,44 +40,128 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] if(no_vvvv_integrals)then integer :: i,j,k,l +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 4 + ! + print*, '' + print*, '' do i = 1,N_int mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) - mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3) - mask_ijkl(i,4) = core_inact_act_bitmask_4(i,4) - enddo - call add_integrals_to_map(mask_ijkl) - ! (core+inact+act) ^ 3 (virt) ^1 - do i = 1,N_int - mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) - mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3) - mask_ijkl(i,4) = virt_bitmask(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) enddo call add_integrals_to_map(mask_ijkl) + call set_integrals_exchange_jj_into_map + call set_integrals_jj_into_map +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 2 (virt) ^2 + ! = J_iv + print*, '' + print*, '' do i = 1,N_int mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) - mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) - mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + + ! (core+inact+act) ^ 2 (virt) ^2 + ! = (iv|iv) + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! + print*, '' + print*, '' + do i = 1,N_int + mask_ijk(i,1) = virt_bitmask(i,1) + mask_ijk(i,2) = virt_bitmask(i,1) + mask_ijk(i,3) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_three_indices(mask_ijk) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 3 (virt) ^1 + ! + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) mask_ijkl(i,4) = virt_bitmask(i,1) enddo call add_integrals_to_map(mask_ijkl) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! ! (core+inact+act) ^ 1 (virt) ^3 + ! + print*, '' + print*, '' do i = 1,N_int mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) mask_ijkl(i,2) = virt_bitmask(i,1) mask_ijkl(i,3) = virt_bitmask(i,1) mask_ijkl(i,4) = virt_bitmask(i,1) enddo - call add_integrals_to_map(mask_ijkl) - + call add_integrals_to_map_no_exit_34(mask_ijkl) + else call add_integrals_to_map(full_ijkl_bitmask_4) endif END_PROVIDER +subroutine set_integrals_jj_into_map + use bitmasks + implicit none + integer :: i,j,n_integrals,i0,j0 + double precision :: buffer_value(mo_tot_num * mo_tot_num) + integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num) + n_integrals = 0 + do j0 = 1, n_virt_orb + j = list_virt(j0) + do i0 = j0, n_virt_orb + i = list_virt(i0) + n_integrals += 1 +! mo_bielec_integral_jj_exchange(i,j) = mo_bielec_integral_vv_exchange_from_ao(i,j) + call mo_bielec_integrals_index(i,j,i,j,buffer_i(n_integrals)) + buffer_value(n_integrals) = mo_bielec_integral_vv_from_ao(i,j) + enddo + enddo + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + call map_unique(mo_integrals_map) +end + +subroutine set_integrals_exchange_jj_into_map + use bitmasks + implicit none + integer :: i,j,n_integrals,i0,j0 + double precision :: buffer_value(mo_tot_num * mo_tot_num) + integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num) + n_integrals = 0 + do j0 = 1, n_virt_orb + j = list_virt(j0) + do i0 = j0+1, n_virt_orb + i = list_virt(i0) + n_integrals += 1 + call mo_bielec_integrals_index(i,j,j,i,buffer_i(n_integrals)) + buffer_value(n_integrals) = mo_bielec_integral_vv_exchange_from_ao(i,j) + enddo + enddo + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + call map_unique(mo_integrals_map) + +end + subroutine add_integrals_to_map(mask_ijkl) use bitmasks implicit none @@ -115,6 +200,672 @@ subroutine add_integrals_to_map(mask_ijkl) !Get list of MOs for i,j,k and l !------------------------------- + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + character*(2048) :: output(1) + print*, 'i' + call bitstring_to_str( output(1), mask_ijkl(1,1), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,1)) + enddo + if(j==0)then + return + endif + + print*, 'j' + call bitstring_to_str( output(1), mask_ijkl(1,2), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,2)) + enddo + if(j==0)then + return + endif + + print*, 'k' + call bitstring_to_str( output(1), mask_ijkl(1,3), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,3)) + enddo + if(j==0)then + return + endif + + print*, 'l' + call bitstring_to_str( output(1), mask_ijkl(1,4), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijkl(i,4)) + enddo + if(j==0)then + return + endif + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + double precision :: accu_bis + accu_bis = 0.d0 + + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num,accu_bis) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 +!$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num +!IRP_IF COARRAY +! if (mod(l1-this_image(),num_images()) /= 0 ) then +! cycle +! endif +!IRP_ENDIF + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_l + l = list_ijkl(l0,4) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + do j0 = 1, n_j +! print*, 'l :: j0',l + j = list_ijkl(j0,2) +! print*, 'j :: 2',j + if (j > l) then +! print*, 'j>l' +! print*, j,l + exit + endif + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) +! print*, 'l :: k0',l +! print*, 'k :: 3',j + i1 = ishft((k*k-k),-1) + if (i1<=j1) then + continue + else +! print*, 'k>l' +! print*, k,l + exit + endif + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) +! print*, 'l :: i0',l +! print*, 'i :: 1',i + if (i>k) then +! print*, 'i>k' +! print*, i,k + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) +! i1+=1 + enddo + +! do i = 1, min(k,j1-i1+list_ijkl(1,1)) +! do i = 1, min(k,j1-i1+list_ijkl(1,1)-1) + do i0 = 1, n_i + i = list_ijkl(i0,1) + if(i> min(k,j1-i1+list_ijkl(1,1)-1))then +! if (i>k) then !min(k,j1-i1) + exit + endif +! print*, i,j,k,l +! print*, k,j1,i1,j1-i1 + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif +! print*, i,j,k,l + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if(i==12.and.k==12 .and.j==12.and.l==12)then +! print*, i,j,k,l,buffer_i(n_integrals) +! accu_bis += buffer_value(n_integrals) +! print*, buffer_value(n_integrals),accu_bis +! endif + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + integer :: index_needed + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL +!IRP_IF COARRAY +! print*, 'Communicating the map' +! call communicate_mo_integrals() +!IRP_ENDIF + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + integer(map_size_kind) :: map_idx + map_idx = ishft(106,map_shift) +! call get_cache_map_verbose(mo_integrals_map,map_idx) + + if (write_mo_integrals) then + call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') + call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + endif + +end + + +subroutine add_integrals_to_map_three_indices(mask_ijk) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijk(N_int,3) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijk(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijk(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijk(1,3), list_ijkl(1,3), n_k, N_int ) + character*(2048) :: output(1) + print*, 'i' + call bitstring_to_str( output(1), mask_ijk(1,1), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,1)) + enddo + if(j==0)then + return + endif + + print*, 'j' + call bitstring_to_str( output(1), mask_ijk(1,2), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,2)) + enddo + if(j==0)then + return + endif + + print*, 'k' + call bitstring_to_str( output(1), mask_ijk(1,3), N_int ) + print *, trim(output(1)) + j = 0 + do i = 1, N_int + j += popcnt(mask_ijk(i,3)) + enddo + if(j==0)then + return + endif + + size_buffer = min(ao_num*ao_num*ao_num,16000000) + print*, 'Providing the molecular integrals ' + print*, 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + double precision :: accu_bis + accu_bis = 0.d0 + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & + !$OMP wall_0,thread_num,accu_bis) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + n_integrals = 0 + wall_0 = wall_1 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + thread_num = 0 +!$ thread_num = omp_get_thread_num() + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num +!IRP_IF COARRAY +! if (mod(l1-this_image(),num_images()) /= 0 ) then +! cycle +! endif +!IRP_ENDIF + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_j + l = list_ijkl(l0,2) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + j0 = l0 + j = list_ijkl(j0,2) + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) +! if (i1<=j1) then +! continue +! else +! exit +! endif + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + enddo + +! do i = 1, min(k,j1-i1+list_ijkl(1,1)) +! do i = 1, min(k,j1-i1+list_ijkl(1,1)-1) + do i0 = 1, n_i + i = list_ijkl(i0,1) +! if(i> min(k,j1-i1+list_ijkl(1,1)-1))then + if (i==k) then !min(k,j1-i1) + cycle + endif +! print*, i,j,k,l +! print*, k,j1,i1,j1-i1 + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif +! print*, i,j,k,l + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if(i==12.and.k==12 .and.j==12.and.l==12)then +! print*, i,j,k,l,buffer_i(n_integrals) +! accu_bis += buffer_value(n_integrals) +! print*, buffer_value(n_integrals),accu_bis +! endif + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(l1)/float(ao_num), '% in ', & + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + endif + endif + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + integer :: index_needed + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL +!IRP_IF COARRAY +! print*, 'Communicating the map' +! call communicate_mo_integrals() +!IRP_ENDIF + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + print*,'Molecular integrals provided:' + print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + print*,' Number of MO integrals: ', mo_map_size + print*,' cpu time :',cpu_2 - cpu_1, 's' + print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + integer(map_size_kind) :: map_idx + map_idx = ishft(106,map_shift) +! call get_cache_map_verbose(mo_integrals_map,map_idx) + + if (write_mo_integrals) then + call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') + call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") + endif + + + +end + + +subroutine add_integrals_to_map_no_exit_34(mask_ijkl) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + PROVIDE ao_bielec_integrals_in_map mo_coef + + !Get list of MOs for i,j,k and l + !------------------------------- + allocate(list_ijkl(mo_tot_num,4)) call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) @@ -268,11 +1019,6 @@ subroutine add_integrals_to_map(mask_ijkl) do k0 = 1, n_k k = list_ijkl(k0,3) i1 = ishft((k*k-k),-1) - if (i1<=j1) then - continue - else - exit - endif bielec_tmp_1 = 0.d0 do i0 = 1, n_i i = list_ijkl(i0,1) @@ -282,7 +1028,12 @@ subroutine add_integrals_to_map(mask_ijkl) bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) enddo - do i = 1, min(k,j1-i1) + do i0 = 1, n_i + i = list_ijkl(i0,1) + if(i> k)then + exit + endif + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then cycle endif @@ -345,8 +1096,6 @@ end - - BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_from_ao, (mo_tot_num_align,mo_tot_num) ] &BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] &BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] @@ -478,6 +1227,155 @@ end mo_bielec_integral_jj_anti_from_ao = mo_bielec_integral_jj_from_ao - mo_bielec_integral_jj_exchange_from_ao +END_PROVIDER + + BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_vv_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! mo_bielec_integral_vv_from_ao(i,j) = J_ij + ! mo_bielec_integral_vv_exchange_from_ao(i,j) = J_ij + ! mo_bielec_integral_vv_anti_from_ao(i,j) = J_ij - K_ij + ! but only for the virtual orbitals + END_DOC + + integer :: i,j,p,q,r,s + integer :: i0,j0 + double precision :: c + real(integral_kind) :: integral + integer :: n, pp + real(integral_kind), allocatable :: int_value(:) + integer, allocatable :: int_idx(:) + + double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) + + if (.not.do_direct_integrals) then + PROVIDE ao_bielec_integrals_in_map mo_coef + endif + + mo_bielec_integral_vv_from_ao = 0.d0 + mo_bielec_integral_vv_exchange_from_ao = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i0,j0,i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & + !$OMP iqrs, iqsr,iqri,iqis) & + !$OMP SHARED(n_virt_orb,mo_tot_num,list_virt,mo_coef_transp,mo_tot_num_align,ao_num,& + !$OMP ao_integrals_threshold,do_direct_integrals) & + !$OMP REDUCTION(+:mo_bielec_integral_vv_from_ao,mo_bielec_integral_vv_exchange_from_ao) + + allocate( int_value(ao_num), int_idx(ao_num), & + iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& + iqsr(mo_tot_num_align,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqrs(i,j) = 0.d0 + iqsr(i,j) = 0.d0 + enddo + enddo + + if (do_direct_integrals) then + double precision :: ao_bielec_integral + do r=1,ao_num + call compute_ao_bielec_integrals(q,r,s,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call compute_ao_bielec_integrals(q,s,r,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i =list_virt(i0) + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + + else + + do r=1,ao_num + call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i =list_virt(i0) + iqrs(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral = int_value(pp) + if (abs(integral) > ao_integrals_threshold) then + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqsr(i,r) += mo_coef_transp(i,p) * integral + enddo + endif + enddo + enddo + endif + iqis = 0.d0 + iqri = 0.d0 + do r=1,ao_num + !DIR$ VECTOR ALIGNED + do i0=1,n_virt_orb + i = list_virt(i0) + iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) + iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) + enddo + enddo + do i0=1,n_virt_orb + i= list_virt(i0) + !DIR$ VECTOR ALIGNED + do j0=1,n_virt_orb + j = list_virt(j0) + c = mo_coef_transp(j,q)*mo_coef_transp(j,s) + mo_bielec_integral_vv_from_ao(j,i) += c * iqis(i) + mo_bielec_integral_vv_exchange_from_ao(j,i) += c * iqri(i) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs,iqsr,int_value,int_idx) + !$OMP END PARALLEL + + mo_bielec_integral_vv_anti_from_ao = mo_bielec_integral_vv_from_ao - mo_bielec_integral_vv_exchange_from_ao +! print*, '**********' +! do i0 =1, n_virt_orb +! i = list_virt(i0) +! print*, mo_bielec_integral_vv_from_ao(i,i) +! enddo +! print*, '**********' + + END_PROVIDER @@ -495,16 +1393,48 @@ END_PROVIDER double precision :: get_mo_bielec_integral PROVIDE mo_bielec_integrals_in_map - mo_bielec_integral_jj = 0.d0 mo_bielec_integral_jj_exchange = 0.d0 - do j=1,mo_tot_num - do i=1,mo_tot_num - mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) - mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) - mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + + if(.not.no_vvvv_integrals)then + do j=1,mo_tot_num + do i=1,mo_tot_num + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + enddo + else + integer :: j0,i0 + do j0=1,n_core_inact_act_orb + j = list_core_inact_act(j0) + do i0=1,n_core_inact_act_orb + i = list_core_inact_act(i0) + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) enddo - enddo + enddo + do j0 = 1, n_virt_orb + j = list_virt(j0) + do i0 = 1, n_virt_orb + i = list_virt(i0) + mo_bielec_integral_jj(i,j) = mo_bielec_integral_vv_from_ao(i,j) + mo_bielec_integral_jj_exchange(i,j) = mo_bielec_integral_vv_exchange_from_ao(i,j) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + do i0=1,n_core_inact_act_orb + i = list_core_inact_act(i0) + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + mo_bielec_integral_jj(j,i) = mo_bielec_integral_jj(i,j) + mo_bielec_integral_jj_exchange(j,i) = mo_bielec_integral_jj_exchange(i,j) + mo_bielec_integral_jj_anti(j,i) = mo_bielec_integral_jj_anti(i,j) + enddo + enddo + + endif END_PROVIDER @@ -518,7 +1448,7 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to do i=1,mo_tot_num do k=1,mo_tot_num - mo_bielec_integral_schwartz(k,i) = dsqrt(mo_bielec_integral_jj(k,i)) + mo_bielec_integral_schwartz(k,i) = 1.d10 enddo enddo diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index 47adc83e..af2af34f 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -484,6 +484,7 @@ subroutine map_get(map, key, value) integer(map_size_kind) :: idx_cache integer(cache_map_size_kind) :: idx + ! index in tha pointers array idx_cache = ishft(key,map_shift) !DIR$ FORCEINLINE call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx) @@ -853,3 +854,25 @@ subroutine get_cache_map(map,map_idx,keys,values,n_elements) enddo end + +subroutine get_cache_map_verbose(map,map_idx) + use map_module + implicit none + type (map_type), intent(in) :: map + integer(map_size_kind), intent(in) :: map_idx + integer(cache_map_size_kind) :: n_elements + integer(key_kind) :: keys(2**16) + double precision :: values(2**16) + integer(cache_map_size_kind) :: i + integer(key_kind) :: shift + + shift = ishft(map_idx,-map_shift) + + n_elements = map%map(map_idx)%n_elements + do i=1,n_elements + keys(i) = map%map(map_idx)%key(i) + shift + values(i) = map%map(map_idx)%value(i) + print*, ',key,values',keys(i),values(i) + enddo + +end