From 4cef44732eb67945688436d38a58f8b3012b25c7 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 14 Mar 2016 16:01:55 +0100 Subject: [PATCH] Program FOBO-SCF works --- plugins/FOBOCI/EZFIO.cfg | 25 ++- plugins/FOBOCI/all_singles.irp.f | 2 +- plugins/FOBOCI/collect_all_lmct.irp.f | 4 +- plugins/FOBOCI/corr_energy_2h2p.irp.f | 64 +++++- plugins/FOBOCI/dress_simple.irp.f | 34 +-- plugins/FOBOCI/fobo_scf.irp.f | 26 ++- .../foboci_lmct_mlct_threshold_old.irp.f | 27 ++- plugins/FOBOCI/new_approach.irp.f | 6 +- plugins/FOBOCI/routines_foboci.irp.f | 211 +++++++++++++++++- 9 files changed, 357 insertions(+), 42 deletions(-) diff --git a/plugins/FOBOCI/EZFIO.cfg b/plugins/FOBOCI/EZFIO.cfg index d4a10add..88189608 100644 --- a/plugins/FOBOCI/EZFIO.cfg +++ b/plugins/FOBOCI/EZFIO.cfg @@ -1,6 +1,13 @@ -[threshold_singles] +[threshold_lmct] type: double precision -doc: threshold to select the pertinent single excitations at second order +doc: threshold to select the pertinent LMCT excitations at second order +interface: ezfio,provider,ocaml +default: 0.01 + + +[threshold_mlct] +type: double precision +doc: threshold to select the pertinent MLCT excitations at second order interface: ezfio,provider,ocaml default: 0.01 @@ -16,6 +23,20 @@ doc: if true, you do the FOBOCI calculation perturbatively interface: ezfio,provider,ocaml default: .False. + +[speed_up_convergence_foboscf] +type: logical +doc: if true, the threshold of the FOBO-SCF algorithms are increased with the iterations +interface: ezfio,provider,ocaml +default: .True. + + +[dressing_2h2p] +type: logical +doc: if true, you do dress with 2h2p excitations each FOBOCI matrix +interface: ezfio,provider,ocaml +default: .False. + [second_order_h] type: logical doc: if true, you do the FOBOCI calculation using second order intermediate Hamiltonian diff --git a/plugins/FOBOCI/all_singles.irp.f b/plugins/FOBOCI/all_singles.irp.f index c4f0b7ae..0594e56e 100644 --- a/plugins/FOBOCI/all_singles.irp.f +++ b/plugins/FOBOCI/all_singles.irp.f @@ -66,7 +66,7 @@ subroutine all_single print*,'E = ',CI_energy(i) print*,'S^2 = ',CI_eigenvectors_s2(i) enddo - do i = 1, 2 + do i = 1, max(2,N_det_generators) print*,'psi_coef = ',psi_coef(i,1) enddo deallocate(pt2,norm_pert,E_before) diff --git a/plugins/FOBOCI/collect_all_lmct.irp.f b/plugins/FOBOCI/collect_all_lmct.irp.f index ebece3ed..96eb2858 100644 --- a/plugins/FOBOCI/collect_all_lmct.irp.f +++ b/plugins/FOBOCI/collect_all_lmct.irp.f @@ -92,7 +92,7 @@ subroutine collect_lmct_mlct(hole_particle,n_couples) iorb = list_act(i) do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then n_couples +=1 hole_particle(n_couples,1) = jorb hole_particle(n_couples,2) = iorb @@ -102,7 +102,7 @@ subroutine collect_lmct_mlct(hole_particle,n_couples) enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then n_couples +=1 hole_particle(n_couples,1) = iorb hole_particle(n_couples,2) = jorb diff --git a/plugins/FOBOCI/corr_energy_2h2p.irp.f b/plugins/FOBOCI/corr_energy_2h2p.irp.f index 8c4f2fe3..ada46bf2 100644 --- a/plugins/FOBOCI/corr_energy_2h2p.irp.f +++ b/plugins/FOBOCI/corr_energy_2h2p.irp.f @@ -1,8 +1,16 @@ BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_ab_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_bb_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_a, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_b, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_double, (mo_tot_num,mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)] &BEGIN_PROVIDER [ double precision, total_corr_e_2h2p] use bitmasks + print*,'' + print*,'Providing the 2h2p correlation energy' + print*,'' implicit none integer(bit_kind) :: key_tmp(N_int,2) integer :: i,j,k,l @@ -12,9 +20,14 @@ integer :: i_ok,ispin ! Alpha - Beta correlation energy total_corr_e_2h2p = 0.d0 + corr_energy_2h2p_ab_2_orb = 0.d0 + corr_energy_2h2p_bb_2_orb = 0.d0 corr_energy_2h2p_per_orb_ab = 0.d0 corr_energy_2h2p_per_orb_aa = 0.d0 corr_energy_2h2p_per_orb_bb = 0.d0 + corr_energy_2h2p_for_1h1p_a = 0.d0 + corr_energy_2h2p_for_1h1p_b = 0.d0 + corr_energy_2h2p_for_1h1p_double = 0.d0 do i = 1, n_inact_orb ! beta i_hole = list_inact(i) do k = 1, n_virt_orb ! beta @@ -36,8 +49,24 @@ hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) contrib = hij*hij/delta_e total_corr_e_2h2p += contrib + ! Single orbital contribution corr_energy_2h2p_per_orb_ab(i_hole) += contrib corr_energy_2h2p_per_orb_ab(k_part) += contrib + ! Couple of orbital contribution for the single 1h1p + corr_energy_2h2p_for_1h1p_a(j_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_a(l_part,j_hole) += contrib + corr_energy_2h2p_for_1h1p_b(j_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_b(l_part,j_hole) += contrib + ! Couple of orbital contribution for the double 1h1p + corr_energy_2h2p_for_1h1p_double(i_hole,l_part) += contrib + corr_energy_2h2p_for_1h1p_double(l_part,i_hole) += contrib + + corr_energy_2h2p_ab_2_orb(i_hole,j_hole) += contrib + corr_energy_2h2p_ab_2_orb(j_hole,i_hole) += contrib + corr_energy_2h2p_ab_2_orb(i_hole,k_part) += contrib + corr_energy_2h2p_ab_2_orb(k_part,i_hole) += contrib + corr_energy_2h2p_ab_2_orb(k_part,l_part) += contrib + corr_energy_2h2p_ab_2_orb(l_part,k_part) += contrib enddo enddo enddo @@ -65,8 +94,12 @@ hij = hij - exc contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) total_corr_e_2h2p += contrib + ! Single orbital contribution corr_energy_2h2p_per_orb_aa(i_hole) += contrib corr_energy_2h2p_per_orb_aa(k_part) += contrib + ! Couple of orbital contribution for the single 1h1p + corr_energy_2h2p_for_1h1p_a(i_hole,k_part) += contrib + corr_energy_2h2p_for_1h1p_a(k_part,i_hole) += contrib enddo enddo enddo @@ -94,8 +127,20 @@ hij = hij - exc contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) total_corr_e_2h2p += contrib + ! Single orbital contribution corr_energy_2h2p_per_orb_bb(i_hole) += contrib corr_energy_2h2p_per_orb_bb(k_part) += contrib + corr_energy_2h2p_for_1h1p_b(i_hole,k_part) += contrib + corr_energy_2h2p_for_1h1p_b(k_part,i_hole) += contrib + + ! Two particle correlation energy + corr_energy_2h2p_bb_2_orb(i_hole,j_hole) += contrib + corr_energy_2h2p_bb_2_orb(j_hole,i_hole) += contrib + corr_energy_2h2p_bb_2_orb(i_hole,k_part) += contrib + corr_energy_2h2p_bb_2_orb(k_part,i_hole) += contrib + corr_energy_2h2p_bb_2_orb(k_part,l_part) += contrib + corr_energy_2h2p_bb_2_orb(l_part,k_part) += contrib + enddo enddo enddo @@ -103,7 +148,11 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)] + BEGIN_PROVIDER [double precision, corr_energy_2h1p_ab_bb_per_2_orb, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_a, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_b, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_double, (mo_tot_num,mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)] &BEGIN_PROVIDER [ double precision, total_corr_e_2h1p] @@ -120,6 +169,10 @@ END_PROVIDER corr_energy_2h1p_per_orb_ab = 0.d0 corr_energy_2h1p_per_orb_aa = 0.d0 corr_energy_2h1p_per_orb_bb = 0.d0 + corr_energy_2h1p_ab_bb_per_2_orb = 0.d0 + corr_energy_2h1p_for_1h1p_a = 0.d0 + corr_energy_2h1p_for_1h1p_b = 0.d0 + corr_energy_2h1p_for_1h1p_double = 0.d0 do i = 1, n_inact_orb i_hole = list_inact(i) do k = 1, n_act_orb @@ -141,6 +194,7 @@ END_PROVIDER hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map) contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) total_corr_e_2h1p += contrib + corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib corr_energy_2h1p_per_orb_ab(i_hole) += contrib corr_energy_2h1p_per_orb_ab(l_part) += contrib enddo @@ -199,6 +253,7 @@ END_PROVIDER delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) hij = hij - exc contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij)) + corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib total_corr_e_2h1p += contrib corr_energy_2h1p_per_orb_bb(i_hole) += contrib @@ -212,6 +267,7 @@ END_PROVIDER BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)] +&BEGIN_PROVIDER [double precision, corr_energy_1h2p_two_orb, (mo_tot_num,mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)] &BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)] &BEGIN_PROVIDER [ double precision, total_corr_e_1h2p] @@ -252,6 +308,8 @@ END_PROVIDER total_corr_e_1h2p += contrib corr_energy_1h2p_per_orb_ab(i_hole) += contrib corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib enddo enddo enddo @@ -282,6 +340,8 @@ END_PROVIDER total_corr_e_1h2p += contrib corr_energy_1h2p_per_orb_aa(i_hole) += contrib corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib enddo enddo enddo @@ -312,6 +372,8 @@ END_PROVIDER total_corr_e_1h2p += contrib corr_energy_1h2p_per_orb_bb(i_hole) += contrib corr_energy_1h2p_per_orb_ab(j_hole) += contrib + corr_energy_1h2p_two_orb(k_part,l_part) += contrib + corr_energy_1h2p_two_orb(l_part,k_part) += contrib enddo enddo enddo diff --git a/plugins/FOBOCI/dress_simple.irp.f b/plugins/FOBOCI/dress_simple.irp.f index 9df94140..99566a8e 100644 --- a/plugins/FOBOCI/dress_simple.irp.f +++ b/plugins/FOBOCI/dress_simple.irp.f @@ -58,24 +58,24 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen call i_h_j(det_buffer(1,1,i),det_buffer(1,1,i),Nint,haa) f = 1.d0/(E_ref-haa) - if(second_order_h)then +! if(second_order_h)then lambda_i = f - else - ! You write the new Hamiltonian matrix - do k = 1, Ndet_generators - H_matrix_tmp(k,Ndet_generators+1) = H_array(k) - H_matrix_tmp(Ndet_generators+1,k) = H_array(k) - enddo - H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa - ! Then diagonalize it - call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1) - ! Then you extract the effective denominator - accu = 0.d0 - do k = 1, Ndet_generators - accu += eigenvectors(k,1) * H_array(k) - enddo - lambda_i = eigenvectors(Ndet_generators+1,1)/accu - endif +! else +! ! You write the new Hamiltonian matrix +! do k = 1, Ndet_generators +! H_matrix_tmp(k,Ndet_generators+1) = H_array(k) +! H_matrix_tmp(Ndet_generators+1,k) = H_array(k) +! enddo +! H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa +! ! Then diagonalize it +! call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1) +! ! Then you extract the effective denominator +! accu = 0.d0 +! do k = 1, Ndet_generators +! accu += eigenvectors(k,1) * H_array(k) +! enddo +! lambda_i = eigenvectors(Ndet_generators+1,1)/accu +! endif do k=1,idx(0) contrib = H_array(idx(k)) * H_array(idx(k)) * lambda_i delta_ij_generators_(idx(k), idx(k)) += contrib diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f index 1b134733..8656b633 100644 --- a/plugins/FOBOCI/fobo_scf.irp.f +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -1,8 +1,8 @@ program foboscf implicit none + call run_prepare no_oa_or_av_opt = .True. touch no_oa_or_av_opt - call run_prepare call routine_fobo_scf call save_mos @@ -10,6 +10,8 @@ end subroutine run_prepare implicit none + no_oa_or_av_opt = .False. + touch no_oa_or_av_opt call damping_SCF call diag_inactive_virt_and_update_mos end @@ -22,6 +24,28 @@ subroutine routine_fobo_scf character*(64) :: label label = "Natural" do i = 1, 5 + print*,'*******************************************************************************' + print*,'*******************************************************************************' + print*,'FOBO-SCF Iteration ',i + print*,'*******************************************************************************' + print*,'*******************************************************************************' + if(speed_up_convergence_foboscf)then + if(i==3)then + threshold_lmct = max(threshold_lmct,0.001) + threshold_mlct = max(threshold_mlct,0.05) + soft_touch threshold_lmct threshold_mlct + endif + if(i==4)then + threshold_lmct = max(threshold_lmct,0.005) + threshold_mlct = max(threshold_mlct,0.07) + soft_touch threshold_lmct threshold_mlct + endif + if(i==5)then + threshold_lmct = max(threshold_lmct,0.01) + threshold_mlct = max(threshold_mlct,0.1) + soft_touch threshold_lmct threshold_mlct + endif + endif call FOBOCI_lmct_mlct_old_thr call save_osoci_natural_mos call damping_SCF diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index b406284f..dc6519b8 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -9,12 +9,9 @@ subroutine FOBOCI_lmct_mlct_old_thr double precision :: norm_tmp(N_states),norm_total(N_states) logical :: test_sym double precision :: thr,hij - double precision :: threshold double precision, allocatable :: dressing_matrix(:,:) logical :: verbose,is_ok verbose = .True. - threshold = threshold_singles - print*,'threshold = ',threshold thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -36,11 +33,14 @@ subroutine FOBOCI_lmct_mlct_old_thr print*,'' print*,'' print*,'DOING FIRST LMCT !!' + print*,'Threshold_lmct = ',threshold_lmct integer(bit_kind) , allocatable :: zero_bitmask(:,:) integer(bit_kind) , allocatable :: psi_singles(:,:,:) + logical :: lmct double precision, allocatable :: psi_singles_coef(:,:) allocate( zero_bitmask(N_int,2) ) do i = 1, n_inact_orb + lmct = .True. integer :: i_hole_osoci i_hole_osoci = list_inact(i) print*,'--------------------------' @@ -55,7 +55,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,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle allocate(dressing_matrix(N_det_generators,N_det_generators)) @@ -81,6 +81,9 @@ subroutine FOBOCI_lmct_mlct_old_thr call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single +! if(dressing_2h2p)then +! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct) +! endif ! ! Change the mask of the holes and particles to perform all the ! ! double excitations that starts from the active space in order @@ -148,9 +151,12 @@ subroutine FOBOCI_lmct_mlct_old_thr if(.True.)then print*,'' print*,'DOING THEN THE MLCT !!' + print*,'Threshold_mlct = ',threshold_mlct + lmct = .False. do i = 1, n_virt_orb integer :: i_particl_osoci i_particl_osoci = list_virt(i) + print*,'--------------------------' ! First set the current generators to the one of restart call set_generators_to_generators_restart @@ -172,7 +178,7 @@ 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,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle allocate(dressing_matrix(N_det_generators,N_det_generators)) @@ -187,6 +193,9 @@ subroutine FOBOCI_lmct_mlct_old_thr ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) call all_single +! if(dressing_2h2p)then +! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct) +! endif endif call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci) do k = 1, N_states @@ -221,10 +230,8 @@ subroutine FOBOCI_mlct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - double precision :: threshold logical :: verbose,is_ok verbose = .False. - threshold = 1.d-2 thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -263,7 +270,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,is_ok,verbose) + call is_a_good_candidate(threshold_mlct,is_ok,verbose) print*,'is_ok = ',is_ok is_ok =.True. if(.not.is_ok)cycle @@ -297,10 +304,8 @@ subroutine FOBOCI_lmct_old double precision :: norm_tmp,norm_total logical :: test_sym double precision :: thr - double precision :: threshold logical :: verbose,is_ok verbose = .False. - threshold = 1.d-2 thr = 1.d-12 allocate(unpaired_bitmask(N_int,2)) allocate (occ(N_int*bit_kind_size,2)) @@ -337,7 +342,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,is_ok,verbose) + call is_a_good_candidate(threshold_lmct,is_ok,verbose) print*,'is_ok = ',is_ok if(.not.is_ok)cycle ! ! so all the mono excitation on the new generators diff --git a/plugins/FOBOCI/new_approach.irp.f b/plugins/FOBOCI/new_approach.irp.f index 2e551dcd..8e2f2e53 100644 --- a/plugins/FOBOCI/new_approach.irp.f +++ b/plugins/FOBOCI/new_approach.irp.f @@ -46,7 +46,7 @@ subroutine new_approach verbose = .True. - threshold = threshold_singles + threshold = threshold_lmct print*,'threshold = ',threshold thr = 1.d-12 print*,'' @@ -623,14 +623,14 @@ subroutine new_approach print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then + if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb)) endif diff --git a/plugins/FOBOCI/routines_foboci.irp.f b/plugins/FOBOCI/routines_foboci.irp.f index 4def99e2..4aca60d7 100644 --- a/plugins/FOBOCI/routines_foboci.irp.f +++ b/plugins/FOBOCI/routines_foboci.irp.f @@ -387,14 +387,14 @@ subroutine save_osoci_natural_mos print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif @@ -519,14 +519,14 @@ subroutine set_osoci_natural_mos print*,'ACTIVE ORBITAL ',iorb do j = 1, n_inact_orb jorb = list_inact(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_lmct)then print*,'INACTIVE ' print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif enddo do j = 1, n_virt_orb jorb = list_virt(j) - if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then + if(dabs(tmp(iorb,jorb)).gt.threshold_mlct)then print*,'VIRT ' print*,'DM ',iorb,jorb,(tmp(iorb,jorb)) endif @@ -607,3 +607,206 @@ end call print_hcc end + + + subroutine dress_diag_elem_2h1p(dressing_H_mat_elem,ndet,lmct,i_hole) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet,i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 0.and.n_p==0)then ! CAS + dressing_H_mat_elem(i)+= total_corr_e_2h1p + if(lmct)then + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(i_hole) - corr_energy_2h1p_per_orb_bb(i_hole) + endif + endif + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(p1) - corr_energy_2h1p_per_orb_aa(p1) + else if (n_h == 1.and.n_p==1)then ! 1h1p +! if(degree==1)then + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i)+= - corr_energy_2h1p_per_orb_ab(h1) +! else +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) +! dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) +! endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i)+= total_corr_e_2h1p + dressing_H_mat_elem(i) += - corr_energy_2h1p_per_orb_ab(h1) + endif + enddo + + end + + subroutine dress_diag_elem_1h2p(dressing_H_mat_elem,ndet,lmct,i_hole) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet,i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 0.and.n_p==0)then ! CAS + dressing_H_mat_elem(i)+= total_corr_e_1h2p + if(.not.lmct)then + dressing_H_mat_elem(i) += - corr_energy_1h2p_per_orb_ab(i_hole) - corr_energy_1h2p_per_orb_aa(i_hole) + endif + endif + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i)+= total_corr_e_1h2p - corr_energy_1h2p_per_orb_ab(h1) + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i)+= 0.d0 + else if (n_h == 1.and.n_p==1)then ! 1h1p + if(degree==1)then + dressing_H_mat_elem(i)+= total_corr_e_1h2p + dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) + else + dressing_H_mat_elem(i) +=0.d0 + endif +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) +! dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & +! - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) +! dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) +! endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i)+= total_corr_e_1h2p + dressing_H_mat_elem(i)+= - corr_energy_1h2p_per_orb_ab(h1) - corr_energy_1h2p_per_orb_ab(h1) + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i) += 0.d0 + endif + enddo + + end + + subroutine dress_diag_elem_2h2p(dressing_H_mat_elem,ndet) + use bitmasks + double precision, intent(inout) :: dressing_H_mat_elem(Ndet) + integer, intent(in) :: ndet + implicit none + integer :: i + integer :: n_p,n_h,number_of_holes,number_of_particles + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2 + do i = 1, N_det + dressing_H_mat_elem(i)+= total_corr_e_2h2p + + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + call get_excitation(ref_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if (n_h == 1.and.n_p==0)then ! 1h + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + else if (n_h == 0.and.n_p==1)then ! 1p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + else if (n_h == 1.and.n_p==1)then ! 1h1p + if(degree==1)then + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_a(h1,p1) + corr_energy_2h2p_for_1h1p_b(h1,p1)) + else + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) + dressing_H_mat_elem(i) += 0.5d0 * (corr_energy_2h2p_for_1h1p_double(h1,p1)) + endif + else if (n_h == 2.and.n_p==1)then ! 2h1p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1) & + - corr_energy_2h2p_per_orb_ab(h2) & + - 0.5d0 * ( corr_energy_2h2p_per_orb_bb(h2) + corr_energy_2h2p_per_orb_bb(h2)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) + if(s1.ne.s2)then + dressing_H_mat_elem(i) += corr_energy_2h2p_ab_2_orb(h1,h2) + else + dressing_H_mat_elem(i) += corr_energy_2h2p_bb_2_orb(h1,h2) + endif + else if (n_h == 1.and.n_p==2)then ! 1h2p + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(h1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(h1) + corr_energy_2h2p_per_orb_bb(h1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p1) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p1) + corr_energy_2h2p_per_orb_bb(p1)) + dressing_H_mat_elem(i) += - corr_energy_2h2p_per_orb_ab(p2) & + - 0.5d0 * (corr_energy_2h2p_per_orb_aa(p2) + corr_energy_2h2p_per_orb_bb(p2)) + if(s1.ne.s2)then + dressing_H_mat_elem(i) += corr_energy_2h2p_ab_2_orb(p1,p2) + else + dressing_H_mat_elem(i) += corr_energy_2h2p_bb_2_orb(p1,p2) + endif + endif + enddo + + end + + subroutine diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole,lmct) + implicit none + double precision, allocatable :: dressing_H_mat_elem(:),energies(:) + integer, intent(in) :: i_hole + logical, intent(in) :: lmct + ! if lmct = .True. ===> LMCT + ! else ===> MLCT + integer :: i + double precision :: hij + allocate(dressing_H_mat_elem(N_det),energies(N_states_diag)) + print*,'' + print*,'dressing with the 2h2p in a CC logic' + print*,'' + do i = 1, N_det + call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij) + dressing_H_mat_elem(i) = hij + enddo + call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det) + call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,i_hole) + call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,i_hole) + call davidson_diag_hjj(psi_det,psi_coef,dressing_H_mat_elem,energies,size(psi_coef,1),N_det,N_states_diag,N_int,output_determinants) + do i = 1, 2 + print*,'psi_coef = ',psi_coef(i,1) + enddo + + + deallocate(dressing_H_mat_elem) + + + + end