From 2e45413f44ae3cf7d1aea0167efd4cfca915b406 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 19 Jan 2023 17:59:10 +0100 Subject: [PATCH 01/14] added optimization for Slater_tc in two-e elements --- src/bi_ort_ints/total_twoe_pot.irp.f | 26 ++++ src/determinants/slater_rules.irp.f | 8 +- src/tc_bi_ortho/slater_tc_opt.irp.f | 208 +++++++++++++++++++++++++ src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 119 +------------- 4 files changed, 246 insertions(+), 115 deletions(-) create mode 100644 src/tc_bi_ortho/slater_tc_opt.irp.f diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index 31cf0624..fef43f93 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -199,3 +199,29 @@ END_PROVIDER ! --- + + BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! mo_bi_ortho_tc_two_e_jj(i,j) = J_ij = + ! mo_bi_ortho_tc_two_e_jj_exchange(i,j) = K_ij = + ! mo_bi_ortho_tc_two_e_jj_anti(i,j) = J_ij - K_ij + END_DOC + + integer :: i,j + double precision :: get_two_e_integral + + mo_bi_ortho_tc_two_e_jj = 0.d0 + mo_bi_ortho_tc_two_e_jj_exchange = 0.d0 + + do i=1,mo_num + do j=1,mo_num + mo_bi_ortho_tc_two_e_jj(i,j) = mo_bi_ortho_tc_two_e(j,i,j,i) + mo_bi_ortho_tc_two_e_jj_exchange(i,j) = mo_bi_ortho_tc_two_e(i,j,j,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j) + enddo + enddo + +END_PROVIDER diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index b9710fd1..78607b7c 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -1790,12 +1790,12 @@ double precision function diag_H_mat_elem(det_in,Nint) integer :: tmp(2) !DIR$ FORCEINLINE call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) - ASSERT (tmp(1) == nexc(1)) - ASSERT (tmp(2) == nexc(2)) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta !DIR$ FORCEINLINE call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) - ASSERT (tmp(1) == nexc(1)) - ASSERT (tmp(2) == nexc(2)) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta det_tmp = ref_bitmask do ispin=1,2 diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f new file mode 100644 index 00000000..0374cb81 --- /dev/null +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -0,0 +1,208 @@ + BEGIN_PROVIDER [ double precision, ref_tc_energy_tot] +&BEGIN_PROVIDER [ double precision, ref_tc_energy_1e] +&BEGIN_PROVIDER [ double precision, ref_tc_energy_2e] +&BEGIN_PROVIDER [ double precision, ref_tc_energy_3e] + implicit none + BEGIN_DOC +! Various component of the TC energy for the reference "HF" Slater determinant + END_DOC + double precision :: hmono, htwoe, htot, hthree + call diag_htilde_mu_mat_bi_ortho(N_int,HF_bitmask , hmono, htwoe, htot) + ref_tc_energy_1e = hmono + ref_tc_energy_2e = htwoe + if(three_body_h_tc)then + call diag_htilde_three_body_ints_bi_ort(N_int, HF_bitmask, hthree) + ref_tc_energy_3e = hthree + else + ref_tc_energy_3e = 0.d0 + endif + ref_tc_energy_tot = ref_tc_energy_1e + ref_tc_energy_2e + ref_tc_energy_3e + END_PROVIDER + +subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, htot) + implicit none + BEGIN_DOC + ! Computes $\langle i|H|i \rangle$. + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + double precision, intent(out) :: hmono,htwoe,htot,hthree + + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: i, nexc(2), ispin + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer(bit_kind) :: det_tmp(Nint,2) + integer :: na, nb + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + + nexc(1) = 0 + nexc(2) = 0 + do i=1,Nint + hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),det_in(i,1)) + particle(i,2) = iand(hole(i,2),det_in(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + + if (nexc(1)+nexc(2) == 0) then + htot = ref_tc_energy_tot + return + endif + + !call debug_det(det_in,Nint) + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta + + + det_tmp = ref_bitmask + hmono = ref_tc_energy_1e + htwoe = ref_tc_energy_2e + hthree= ref_tc_energy_3e + + do ispin=1,2 + na = elec_num_tab(ispin) + nb = elec_num_tab(iand(ispin,1)+1) + do i=1,nexc(ispin) + !DIR$ FORCEINLINE + call ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, hmono,htwoe,hthree, Nint,na,nb) + !DIR$ FORCEINLINE + call a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, hmono,htwoe,hthree, Nint,na,nb) + enddo + enddo + htot = hmono+htwoe+hthree +end + +subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes one- and two-body energy corresponding + ! + ! to the ADDITION of an electron in an orbital 'iorb' of spin 'ispin' + ! + ! onto a determinant 'key'. + ! + ! in output, the determinant key is changed by the ADDITION of that electron + ! + ! and the quantities hmono,htwoe,hthree are INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hmono,htwoe,hthree + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + + if (iorb < 1) then + print *, irp_here, ': iorb < 1' + print *, iorb, mo_num + stop -1 + endif + if (iorb > mo_num) then + print *, irp_here, ': iorb > mo_num' + print *, iorb, mo_num + stop -1 + endif + + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + ASSERT (tmp(1) == elec_alpha_num) + ASSERT (tmp(2) == elec_beta_num) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb) + + ! Same spin + do i=1,na + htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) + enddo + na = na+1 +end + +subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes one- and two-body energy corresponding + ! + ! to the REMOVAL of an electron in an orbital 'iorb' of spin 'ispin' + ! + ! onto a determinant 'key'. + ! + ! in output, the determinant key is changed by the REMOVAL of that electron + ! + ! and the quantities hmono,htwoe,hthree are INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hmono,htwoe,hthree + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + integer :: tmp(2) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k>0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + na = na-1 + + hmono = hmono - mo_bi_ortho_tc_one_e(iorb,iorb) + + ! Same spin + do i=1,na + htwoe= htwoe- mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) + enddo + +end diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 2d71b6b2..094c9bbc 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,121 +11,18 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - ! call routine_2 - call test_rout + call test_slater_tc_opt end -subroutine test_rout +subroutine test_slater_tc_opt implicit none - integer :: i,j,ii,jj - use bitmasks ! you need to include the bitmasks_module.f90 features - integer(bit_kind), allocatable :: det_i(:,:) - allocate(det_i(N_int,2)) - det_i(:,:)= psi_det(:,:,1) - call debug_det(det_i,N_int) - integer, allocatable :: occ(:,:) - integer :: n_occ_ab(2) - allocate(occ(N_int*bit_kind_size,2)) - call bitstring_to_list_ab(det_i, occ, n_occ_ab, N_int) - double precision :: hmono, htwoe, htot - call diag_htilde_mu_mat_bi_ortho(N_int, det_i, hmono, htwoe, htot) - print*,'hmono, htwoe, htot' - print*, hmono, htwoe, htot - print*,'alpha electrons orbital occupancy' - do i = 1, n_occ_ab(1) ! browsing the alpha electrons - j = occ(i,1) - print*,j,mo_bi_ortho_tc_one_e(j,j) - enddo - print*,'beta electrons orbital occupancy' - do i = 1, n_occ_ab(2) ! browsing the beta electrons - j = occ(i,2) - print*,j,mo_bi_ortho_tc_one_e(j,j) - enddo - print*,'alpha beta' - do i = 1, n_occ_ab(1) - ii = occ(i,1) - do j = 1, n_occ_ab(2) - jj = occ(j,2) - print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii) - enddo - enddo - print*,'alpha alpha' - do i = 1, n_occ_ab(1) - ii = occ(i,1) - do j = 1, n_occ_ab(1) - jj = occ(j,1) - print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii) - enddo - enddo - - print*,'beta beta' - do i = 1, n_occ_ab(2) - ii = occ(i,2) - do j = 1, n_occ_ab(2) - jj = occ(j,2) - print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii) - enddo - enddo - - -end - -subroutine routine_2 - implicit none - integer :: i - double precision :: bi_ortho_mo_ints - print*,'H matrix' + integer :: i,j + double precision :: hmono, htwoe, htot, hthree + double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot do i = 1, N_det - write(*,'(1000(F16.5,X))')htilde_matrix_elmt_bi_ortho(:,i) - enddo - i = 1 - double precision :: phase - integer :: degree,h1, p1, h2, p2, s1, s2, exc(0:2,2,2) - call get_excitation_degree(ref_bitmask, psi_det(1,1,i), degree, N_int) - if(degree==2)then - call get_double_excitation(ref_bitmask, psi_det(1,1,i), exc, phase, N_int) - call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) - print*,'h1,h2,p1,p2' - print*, h1,h2,p1,p2 - print*,mo_bi_ortho_tc_two_e(p1,p2,h1,h2),mo_bi_ortho_tc_two_e(h1,h2,p1,p2) - endif - - - print*,'coef' - do i = 1, ao_num - print*,i,mo_l_coef(i,8),mo_r_coef(i,8) - enddo -! print*,'mdlqfmlqgmqglj' -! print*,'mo_bi_ortho_tc_two_e()',mo_bi_ortho_tc_two_e(2,2,3,3) -! print*,'bi_ortho_mo_ints ',bi_ortho_mo_ints(2,2,3,3) - print*,'Overlap' - do i = 1, mo_num - write(*,'(100(F16.10,X))')overlap_bi_ortho(:,i) + call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) + call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) + print*,htot,hnewtot,dabs(htot-hnewtot) enddo end - -subroutine routine - implicit none - double precision :: hmono,htwoe,hthree,htot - integer(bit_kind), allocatable :: key1(:,:) - integer(bit_kind), allocatable :: key2(:,:) - allocate(key1(N_int,2),key2(N_int,2)) - use bitmasks - key1 = ref_bitmask - call htilde_mu_mat_bi_ortho(key1,key1, N_int, hmono,htwoe,hthree,htot) - key2 = key1 - integer :: h,p,i_ok - h = 1 - p = 8 - call do_single_excitation(key2,h,p,1,i_ok) - call debug_det(key2,N_int) - call htilde_mu_mat_bi_ortho(key2,key1, N_int, hmono,htwoe,hthree,htot) -! print*,'fock_matrix_tc_mo_alpha(p,h) = ',fock_matrix_tc_mo_alpha(p,h) - print*,'htot = ',htot - print*,'hmono = ',hmono - print*,'htwoe = ',htwoe - double precision :: bi_ortho_mo_ints - print*,'bi_ortho_mo_ints(1,p,1,h)',bi_ortho_mo_ints(1,p,1,h) - -end From 1651242fba9b0450003c32ea6c868b28f77f8443 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 19 Jan 2023 18:45:43 +0100 Subject: [PATCH 02/14] beginning to optimize the single excitations on tc --- src/tc_bi_ortho/slater_tc_opt.irp.f | 1 + src/tc_bi_ortho/slater_tc_opt_single.irp.f | 202 +++++++++++++++++++++ 2 files changed, 203 insertions(+) create mode 100644 src/tc_bi_ortho/slater_tc_opt_single.irp.f diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index 0374cb81..50886fe3 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -206,3 +206,4 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) enddo end + diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f new file mode 100644 index 00000000..a69f5d2e --- /dev/null +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -0,0 +1,202 @@ + + +subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe, htot) + BEGIN_DOC + ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: hmono, htwoe, htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int, phase + double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13 + integer :: other_spin(2) + integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2) + + other_spin(1) = 2 + other_spin(2) = 1 + + hmono = 0.d0 + htwoe= 0.d0 + htot = 0.d0 + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.ne.1)then + return + endif + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + + call get_single_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) + call get_single_excitation_from_fock_tc(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,hthree,htot) +end + + BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)] +&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)] + implicit none + BEGIN_DOC + ! tc_2e_3idx_coulomb_integrals(j,k,i) = + ! + ! tc_2e_3idx_exchange_integrals(j,k,i) = + END_DOC + integer :: i,j,k,l + double precision :: get_two_e_integral + double precision :: integral + + do k = 1, mo_num + do i = 1, mo_num + do j = 1, mo_num + tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) + tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] + implicit none + BEGIN_DOC +! Closed-shell part of the Fock operator for the TC operator + END_DOC + integer :: h0,p0,h,p,k0,k,i + integer :: n_occ_ab(2) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab_virt(2) + integer :: occ_virt(N_int*bit_kind_size,2) + integer(bit_kind) :: key_test(N_int) + integer(bit_kind) :: key_virt(N_int,2) + + call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) + do i = 1, N_int + key_virt(i,1) = full_ijkl_bitmask(i) + key_virt(i,2) = full_ijkl_bitmask(i) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + ! virt ---> virt single excitations + do h0 = 1, n_occ_ab_virt(1) + h=occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + ! docc ---> docc single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + +END_PROVIDER + + +subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,htwoe,hthree,htot) + use bitmasks + implicit none + integer,intent(in) :: h,p,spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) + double precision, intent(out) :: hmono,htwoe,hthree,htot + integer(bit_kind) :: differences(N_int,2) + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: partcl(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_partcl(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i + double precision :: buffer_c(mo_num),buffer_x(mo_num) + do i=1, mo_num + buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) + buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) + enddo + do i = 1, N_int + differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),key_i(i,1)) + partcl(i,2) = iand(differences(i,2),key_i(i,2)) + enddo + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + hmono = mo_bi_ortho_tc_one_e(p,h) + htwoe = fock_op_2_e_tc_closed_shell(p,h) + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + htwoe -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + htwoe -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + htwoe += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + htwoe += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + htwoe += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + htwoe -= buffer_x(i) + enddo + hthree = 0.d0 + htwoe = htwoe * phase + hmono = hmono * phase + htot = htwoe + hmono + hthree + +end + From 4ee080215068c055603c24989f78ba6ceb25a43b Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 19 Jan 2023 19:29:26 +0100 Subject: [PATCH 03/14] two body part up to single excitations work with fock --- src/bi_ort_ints/total_twoe_pot.irp.f | 23 +++ src/determinants/single_excitations.irp.f | 2 +- src/tc_bi_ortho/slater_tc.irp.f | 3 + src/tc_bi_ortho/slater_tc_opt_single.irp.f | 214 ++++++++++++--------- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 20 +- 5 files changed, 165 insertions(+), 97 deletions(-) diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index fef43f93..e74c6d2a 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -224,4 +224,27 @@ END_PROVIDER enddo enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)] +&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)] + implicit none + BEGIN_DOC + ! tc_2e_3idx_coulomb_integrals(j,k,i) = + ! + ! tc_2e_3idx_exchange_integrals(j,k,i) = + END_DOC + integer :: i,j,k,l + double precision :: get_two_e_integral + double precision :: integral + + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) + tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) + enddo + enddo + enddo + END_PROVIDER diff --git a/src/determinants/single_excitations.irp.f b/src/determinants/single_excitations.irp.f index ccfeaa2e..1c25e314 100644 --- a/src/determinants/single_excitations.irp.f +++ b/src/determinants/single_excitations.irp.f @@ -28,7 +28,7 @@ BEGIN_PROVIDER [double precision, fock_operator_closed_shell_ref_bitmask, (mo_nu integer :: occ_virt(N_int*bit_kind_size,2) integer(bit_kind) :: key_test(N_int) integer(bit_kind) :: key_virt(N_int,2) - + fock_operator_closed_shell_ref_bitmask = 0.d0 call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) do i = 1, N_int key_virt(i,1) = full_ijkl_bitmask(i) diff --git a/src/tc_bi_ortho/slater_tc.irp.f b/src/tc_bi_ortho/slater_tc.irp.f index 33b738ba..2c0ae2ca 100644 --- a/src/tc_bi_ortho/slater_tc.irp.f +++ b/src/tc_bi_ortho/slater_tc.irp.f @@ -324,6 +324,9 @@ subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) +! if(h1==14.and.p1==2)then +! print*,'h1,p1 old = ',h1,p1 +! endif hmono = mo_bi_ortho_tc_one_e(p1,h1) * phase diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f index a69f5d2e..df930136 100644 --- a/src/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -1,6 +1,6 @@ -subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe, htot) +subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe, hthree, htot) BEGIN_DOC ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS !! @@ -14,7 +14,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe implicit none integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) - double precision, intent(out) :: hmono, htwoe, htot + double precision, intent(out) :: hmono, htwoe, hthree, htot integer :: occ(Nint*bit_kind_size,2) integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: degree,exc(0:2,2,2) @@ -27,109 +27,21 @@ subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe other_spin(1) = 2 other_spin(2) = 1 - hmono = 0.d0 - htwoe= 0.d0 - htot = 0.d0 + hmono = 0.d0 + htwoe = 0.d0 + hthree = 0.d0 + htot = 0.d0 call get_excitation_degree(key_i, key_j, degree, Nint) if(degree.ne.1)then return endif - call bitstring_to_list_ab(key_i, occ, Ne, Nint) + call bitstring_to_list_ab(key_i, occ, Ne, Nint) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) call get_single_excitation_from_fock_tc(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,hthree,htot) end - BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)] -&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)] - implicit none - BEGIN_DOC - ! tc_2e_3idx_coulomb_integrals(j,k,i) = - ! - ! tc_2e_3idx_exchange_integrals(j,k,i) = - END_DOC - integer :: i,j,k,l - double precision :: get_two_e_integral - double precision :: integral - - do k = 1, mo_num - do i = 1, mo_num - do j = 1, mo_num - tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) - tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) - enddo - enddo - enddo - -END_PROVIDER - - -BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] - implicit none - BEGIN_DOC -! Closed-shell part of the Fock operator for the TC operator - END_DOC - integer :: h0,p0,h,p,k0,k,i - integer :: n_occ_ab(2) - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab_virt(2) - integer :: occ_virt(N_int*bit_kind_size,2) - integer(bit_kind) :: key_test(N_int) - integer(bit_kind) :: key_virt(N_int,2) - - call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) - do i = 1, N_int - key_virt(i,1) = full_ijkl_bitmask(i) - key_virt(i,2) = full_ijkl_bitmask(i) - key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) - key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) - enddo - call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) - ! docc ---> virt single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu - enddo - enddo - - ! virt ---> virt single excitations - do h0 = 1, n_occ_ab_virt(1) - h=occ_virt(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu - enddo - enddo - - ! docc ---> docc single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab(1) - p = occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - array_exchange(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu - enddo - enddo - -END_PROVIDER - subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,htwoe,hthree,htot) use bitmasks @@ -200,3 +112,115 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h end + +BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] + implicit none + BEGIN_DOC +! Closed-shell part of the Fock operator for the TC operator + END_DOC + integer :: h0,p0,h,p,k0,k,i + integer :: n_occ_ab(2) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab_virt(2) + integer :: occ_virt(N_int*bit_kind_size,2) + integer(bit_kind) :: key_test(N_int) + integer(bit_kind) :: key_virt(N_int,2) + double precision :: accu + + fock_op_2_e_tc_closed_shell = -1000.d0 + call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) + do i = 1, N_int + key_virt(i,1) = full_ijkl_bitmask(i) + key_virt(i,2) = full_ijkl_bitmask(i) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab(1) + p=occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + ! virt ---> virt single excitations + do h0 = 1, n_occ_ab_virt(1) + h=occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p=occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + + ! docc ---> docc single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + do h0 = 1, n_occ_ab(1) + h = occ(h0,1) + do p0 = 1, n_occ_ab(1) + p=occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + +! do i = 1, mo_num +! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i) +! enddo + +END_PROVIDER + diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 094c9bbc..069a1d53 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -19,10 +19,28 @@ subroutine test_slater_tc_opt integer :: i,j double precision :: hmono, htwoe, htot, hthree double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot + double precision :: accu ,i_count + accu = 0.d0 + i_count = 0.d0 do i = 1, N_det +! do i = 14,14 call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) - print*,htot,hnewtot,dabs(htot-hnewtot) + do j = 1, N_det +! do j = 1, 1 + if(i==j)cycle + call single_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, htot) + call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) + if(dabs(htot).gt.1.d-10)then +! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then + print*,j,i + i_count += 1.D0 + print*,htot,hnewtot,dabs(htot-hnewtot) + accu += dabs(htot-hnewtot) +! endif + endif + enddo enddo + print*,'accu = ',accu/i_count end From f0178d09a2572990111211fbe4b894bde9fa05ee Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 19 Jan 2023 22:34:11 +0100 Subject: [PATCH 04/14] diagonal matrix elements work with 3-e a la fock --- src/bi_ort_ints/three_body_ijm.irp.f | 2 +- src/tc_bi_ortho/slater_tc_3e.irp.f | 13 ----- src/tc_bi_ortho/slater_tc_opt.irp.f | 78 ++++++++++++++++++++++++-- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 18 ++++-- 4 files changed, 88 insertions(+), 23 deletions(-) diff --git a/src/bi_ort_ints/three_body_ijm.irp.f b/src/bi_ort_ints/three_body_ijm.irp.f index 0e42264b..4d21cb93 100644 --- a/src/bi_ort_ints/three_body_ijm.irp.f +++ b/src/bi_ort_ints/three_body_ijm.irp.f @@ -60,7 +60,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num ! ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation ! - ! three_e_3_idx_direct_bi_ort(m,j,i) = + ! three_e_3_idx_cycle_1_bi_ort(m,j,i) = ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index a56a432f..0d5f8542 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -49,8 +49,6 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) if(Ne(1)+Ne(2).ge.3)then !! ! alpha/alpha/beta three-body - double precision :: accu - accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) do j = i+1, Ne(1) @@ -62,14 +60,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR hthree += direct_int - exchange_int - accu += direct_int - exchange_int enddo enddo enddo - !print*,'aab = ',accu ! beta/beta/alpha three-body - accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) do j = i+1, Ne(2) @@ -79,14 +74,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) hthree += direct_int - exchange_int - accu += direct_int - exchange_int enddo enddo enddo - !print*,'abb = ',accu ! alpha/alpha/alpha three-body - accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) ! 1 do j = i+1, Ne(1) @@ -95,14 +87,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,1) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS - accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo - !print*,'aaa = ',accu ! beta/beta/beta three-body - accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) ! 1 do j = i+1, Ne(2) @@ -111,11 +100,9 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,2) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS - accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo - !print*,'bbb = ',accu endif end diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index 50886fe3..4048f481 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -55,6 +55,9 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, enddo if (nexc(1)+nexc(2) == 0) then + hmono = ref_tc_energy_1e + htwoe = ref_tc_energy_2e + hthree= ref_tc_energy_3e htot = ref_tc_energy_tot return endif @@ -75,7 +78,6 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, hmono = ref_tc_energy_1e htwoe = ref_tc_energy_2e hthree= ref_tc_energy_3e - do ispin=1,2 na = elec_num_tab(ispin) nb = elec_num_tab(iand(ispin,1)+1) @@ -110,7 +112,9 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer :: occ(Nint*bit_kind_size,2) integer :: other_spin - integer :: k,l,i + integer :: k,l,i,jj,mm,j,m + double precision :: three_e_diag_parrallel_spin, direct_int, exchange_int + if (iorb < 1) then print *, irp_here, ': iorb < 1' @@ -151,6 +155,39 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) do i=1,nb htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo + + if(three_body_h_tc)then + !!!!! 3-e part + !! same-spin/same-spin + do j = 1, na + jj = occ(j,ispin) + do m = j+1, na + mm = occ(m,ispin) + hthree += three_e_diag_parrallel_spin(mm,jj,iorb) + enddo + enddo + !! same-spin/oposite-spin + do j = 1, na + jj = occ(j,ispin) + do m = 1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo + enddo + !! oposite-spin/opposite-spin + do j = 1, nb + jj = occ(j,other_spin) + do m = j+1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo + enddo + endif + na = na+1 end @@ -172,10 +209,11 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer, intent(inout) :: na, nb integer(bit_kind), intent(inout) :: key(Nint,2) double precision, intent(inout) :: hmono,htwoe,hthree - + + double precision :: direct_int, exchange_int, three_e_diag_parrallel_spin integer :: occ(Nint*bit_kind_size,2) integer :: other_spin - integer :: k,l,i + integer :: k,l,i,jj,mm,j,m integer :: tmp(2) ASSERT (iorb > 0) @@ -205,5 +243,37 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo + if(three_body_h_tc)then + !!!!! 3-e part + !! same-spin/same-spin + do j = 1, na + jj = occ(j,ispin) + do m = j+1, na + mm = occ(m,ispin) + hthree -= three_e_diag_parrallel_spin(mm,jj,iorb) + enddo + enddo + !! same-spin/oposite-spin + do j = 1, na + jj = occ(j,ispin) + do m = 1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree -= (direct_int - exchange_int) + enddo + enddo + !! oposite-spin/opposite-spin + do j = 1, nb + jj = occ(j,other_spin) + do m = j+1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree -= (direct_int - exchange_int) + enddo + enddo + endif + end diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 069a1d53..dda4bd00 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -19,28 +19,36 @@ subroutine test_slater_tc_opt integer :: i,j double precision :: hmono, htwoe, htot, hthree double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot - double precision :: accu ,i_count + double precision :: accu_d ,i_count, accu accu = 0.d0 + accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det ! do i = 14,14 call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) + call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) +! print*,hthree,hnewthnewree +! print*,htot,hnewtot,dabs(hnewtot-htot) + accu_d += dabs(htot-hnewtot) +! if(dabs(htot-hnewtot).gt.1.d-8)then + print*,i + print*,htot,hnewtot,dabs(htot-hnewtot) +! endif do j = 1, N_det -! do j = 1, 1 if(i==j)cycle call single_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, htot) call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) if(dabs(htot).gt.1.d-10)then -! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then + if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then print*,j,i i_count += 1.D0 print*,htot,hnewtot,dabs(htot-hnewtot) accu += dabs(htot-hnewtot) -! endif + endif endif enddo enddo - print*,'accu = ',accu/i_count + print*,'accu_d = ',accu_d/N_det end From 0011aa2bc2817ad5f519e35c61a2ec667e35223f Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 11:19:38 +0100 Subject: [PATCH 05/14] TC single excitations H matrix elements work with Fock matrix --- src/bi_ort_ints/three_body_ijmk.irp.f | 4 +- src/tc_bi_ortho/slater_tc_opt_single.irp.f | 234 +++++++++++++++++++++ src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 39 ++-- 3 files changed, 262 insertions(+), 15 deletions(-) diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f index 0d5016ce..853972f7 100644 --- a/src/bi_ort_ints/three_body_ijmk.irp.f +++ b/src/bi_ort_ints/three_body_ijmk.irp.f @@ -195,7 +195,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num, ! ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs ! - ! three_e_4_idx_exch13_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! three_e_4_idx_exch13_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign END_DOC @@ -241,7 +241,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num, ! ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs ! - ! three_e_4_idx_exch12_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! three_e_4_idx_exch12_bi_ort(m,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f index df930136..cb9306aa 100644 --- a/src/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -106,12 +106,246 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h htwoe -= buffer_x(i) enddo hthree = 0.d0 + if (three_body_h_tc)then + call three_comp_fock_elem(key_i,h,p,spin,hthree) + endif + + htwoe = htwoe * phase hmono = hmono * phase + hthree = hthree * phase htot = htwoe + hmono + hthree end +subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) + implicit none + integer,intent(in) :: h_fock,p_fock,ispin_fock + integer(bit_kind), intent(in) :: key_i(N_int,2) + double precision, intent(out) :: hthree + integer :: nexc(2),i,ispin,na,nb + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: particle(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_particle(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) + integer(bit_kind) :: det_tmp(N_int,2) + + + nexc(1) = 0 + nexc(2) = 0 + !! Get all the holes and particles of key_i with respect to the ROHF determinant + do i=1,N_int + hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),key_i(i,1)) + particle(i,2) = iand(hole(i,2),key_i(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta + + !! Initialize the matrix element with the reference ROHF Slater determinant Fock element + if(ispin_fock==1)then + hthree = fock_a_tot_3e_bi_orth(p_fock,h_fock) + else + hthree = fock_b_tot_3e_bi_orth(p_fock,h_fock) + endif + det_tmp = ref_bitmask + do ispin=1,2 + na = elec_num_tab(ispin) + nb = elec_num_tab(iand(ispin,1)+1) + do i=1,nexc(ispin) + !DIR$ FORCEINLINE + call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + !DIR$ FORCEINLINE + call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + enddo + enddo +end + +subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes the contribution to the three-electron part of the Fock operator + ! + ! a^dagger_{p_fock} a_{h_fock} of spin ispin_fock + ! + ! on top of a determinant 'key' on which you ADD an electron of spin ispin in orbital iorb + ! + ! in output, the determinant key is changed by the ADDITION of that electron + ! + ! the output hthree is INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint, h_fock,p_fock, ispin_fock + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hthree + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i,jj,j + double precision :: three_e_single_parrallel_spin, direct_int, exchange_int + + + if (iorb < 1) then + print *, irp_here, ': iorb < 1' + print *, iorb, mo_num + stop -1 + endif + if (iorb > mo_num) then + print *, irp_here, ': iorb > mo_num' + print *, iorb, mo_num + stop -1 + endif + + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + ASSERT (tmp(1) == elec_alpha_num) + ASSERT (tmp(2) == elec_beta_num) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + + !! spin of other electrons == ispin + if(ispin == ispin_fock)then + !! in what follows :: jj == other electrons in the determinant + !! :: iorb == electron that has been added of spin ispin + !! :: p_fock, h_fock == hole particle of spin ispin_fock + !! jj = ispin = ispin_fock >> pure parallel spin + do j = 1, na + jj = occ(j,ispin) + hthree += three_e_single_parrallel_spin(jj,iorb,p_fock,h_fock) + enddo + !! spin of jj == other spin than ispin AND ispin_fock + !! exchange between the iorb and (h_fock, p_fock) + do j = 1, nb + jj = occ(j,other_spin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch12_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree += direct_int - exchange_int + enddo + else !! ispin NE to ispin_fock + !! jj = ispin BUT NON EQUAL TO ispin_fock + !! exchange between the jj and iorb + do j = 1, na + jj = occ(j,ispin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch23_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree += direct_int - exchange_int + enddo + !! jj = other_spin than ispin BUT jj == ispin_fock + !! exchange between jj and (h_fock,p_fock) + do j = 1, nb + jj = occ(j,other_spin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch13_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree += direct_int - exchange_int + enddo + endif + + na = na+1 +end + +subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes the contribution to the three-electron part of the Fock operator + ! + ! a^dagger_{p_fock} a_{h_fock} of spin ispin_fock + ! + ! on top of a determinant 'key' on which you REMOVE an electron of spin ispin in orbital iorb + ! + ! in output, the determinant key is changed by the REMOVAL of that electron + ! + ! the output hthree is INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint, h_fock,p_fock, ispin_fock + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hthree + + double precision :: direct_int, exchange_int, three_e_single_parrallel_spin + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i,jj,mm,j,m + integer :: tmp(2) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k>0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + na = na-1 + !! spin of other electrons == ispin + if(ispin == ispin_fock)then + !! in what follows :: jj == other electrons in the determinant + !! :: iorb == electron that has been added of spin ispin + !! :: p_fock, h_fock == hole particle of spin ispin_fock + !! jj = ispin = ispin_fock >> pure parallel spin + do j = 1, na + jj = occ(j,ispin) + hthree -= three_e_single_parrallel_spin(jj,iorb,p_fock,h_fock) + enddo + !! spin of jj == other spin than ispin AND ispin_fock + !! exchange between the iorb and (h_fock, p_fock) + do j = 1, nb + jj = occ(j,other_spin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch12_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree -= direct_int - exchange_int + enddo + else !! ispin NE to ispin_fock + !! jj = ispin BUT NON EQUAL TO ispin_fock + !! exchange between the jj and iorb + do j = 1, na + jj = occ(j,ispin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch23_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree -= direct_int - exchange_int + enddo + !! jj = other_spin than ispin BUT jj == ispin_fock + !! exchange between jj and (h_fock,p_fock) + do j = 1, nb + jj = occ(j,other_spin) + direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + exchange_int = three_e_4_idx_exch13_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR + hthree -= direct_int - exchange_int + enddo + endif + +end + BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] implicit none diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index dda4bd00..f48984ea 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -18,37 +18,50 @@ subroutine test_slater_tc_opt implicit none integer :: i,j double precision :: hmono, htwoe, htot, hthree - double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot + double precision :: hnewmono, hnewtwoe, hnewthree, hnewtot double precision :: accu_d ,i_count, accu accu = 0.d0 accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det -! do i = 14,14 +! do i = 1,1 call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) - call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) -! print*,hthree,hnewthnewree + call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) +! print*,hthree,hnewthree ! print*,htot,hnewtot,dabs(hnewtot-htot) accu_d += dabs(htot-hnewtot) -! if(dabs(htot-hnewtot).gt.1.d-8)then + if(dabs(htot-hnewtot).gt.1.d-8)then print*,i print*,htot,hnewtot,dabs(htot-hnewtot) -! endif - do j = 1, N_det + endif +! do j = 319,319 + do j = 1,N_det if(i==j)cycle - call single_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, htot) - call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) - if(dabs(htot).gt.1.d-10)then - if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i + integer :: degree + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree .ne. 1)cycle + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) +! print*,'j,i',j,i +! print*,htot,hnewtot,dabs(htot-hnewtot) +! print*,hthree,hnewthree,dabs(hthree-hnewthree) + if(dabs(hthree).gt.1.d-15)then +! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then i_count += 1.D0 - print*,htot,hnewtot,dabs(htot-hnewtot) accu += dabs(htot-hnewtot) + if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then + print*,j,i + call debug_det(psi_det(1,1,i),N_int) + call debug_det(psi_det(1,1,j),N_int) +! print*,htot,hnewtot,dabs(htot-hnewtot) + print*,hthree,hnewthree,dabs(hthree-hnewthree) + stop endif endif enddo enddo print*,'accu_d = ',accu_d/N_det + print*,'accu = ',accu/i_count end From 7a144bc1a2ad7cc26dc90310cfd996ea3d464288 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 11:20:26 +0100 Subject: [PATCH 06/14] renamed a file --- src/tc_bi_ortho/{slater_tc_opt.irp.f => slater_tc_opt_diag.irp.f} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/tc_bi_ortho/{slater_tc_opt.irp.f => slater_tc_opt_diag.irp.f} (100%) diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f similarity index 100% rename from src/tc_bi_ortho/slater_tc_opt.irp.f rename to src/tc_bi_ortho/slater_tc_opt_diag.irp.f From 721e0963b9e6cfe6e58e4e7314bd0d8d1230d672 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 11:31:28 +0100 Subject: [PATCH 07/14] working on TC Slater matrix elements --- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index f48984ea..8afb3b25 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -40,22 +40,21 @@ subroutine test_slater_tc_opt if(i==j)cycle integer :: degree call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) - if(degree .ne. 1)cycle +! if(degree .ne. 1)cycle call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) - call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) -! print*,'j,i',j,i -! print*,htot,hnewtot,dabs(htot-hnewtot) -! print*,hthree,hnewthree,dabs(hthree-hnewthree) - if(dabs(hthree).gt.1.d-15)then + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot) +! if(dabs(hthree).gt.1.d-15)then + if(dabs(htot).gt.1.d-15)then ! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then i_count += 1.D0 accu += dabs(htot-hnewtot) - if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then +! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then + if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then print*,j,i call debug_det(psi_det(1,1,i),N_int) call debug_det(psi_det(1,1,j),N_int) -! print*,htot,hnewtot,dabs(htot-hnewtot) - print*,hthree,hnewthree,dabs(hthree-hnewthree) + print*,htot,hnewtot,dabs(htot-hnewtot) +! print*,hthree,hnewthree,dabs(hthree-hnewthree) stop endif endif From f1137bc883fa401d3941dcde119638c78d8bd76a Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 14:57:21 +0100 Subject: [PATCH 08/14] beginning to optimize the double excitations for TC --- src/bi_ort_ints/three_body_ijmkl.irp.f | 6 +-- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 57 ++++++++++++++++++++++ src/tc_bi_ortho/test_normal_order.irp.f | 52 ++++++++++---------- 3 files changed, 87 insertions(+), 28 deletions(-) create mode 100644 src/tc_bi_ortho/slater_tc_opt_double.irp.f diff --git a/src/bi_ort_ints/three_body_ijmkl.irp.f b/src/bi_ort_ints/three_body_ijmkl.irp.f index 6287c5a3..bd5c4977 100644 --- a/src/bi_ort_ints/three_body_ijmkl.irp.f +++ b/src/bi_ort_ints/three_body_ijmkl.irp.f @@ -7,7 +7,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, ! ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs ! - ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign END_DOC @@ -202,7 +202,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num, ! ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs ! - ! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! @@ -251,7 +251,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num, ! ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs ! - ! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f new file mode 100644 index 00000000..9cff8ff3 --- /dev/null +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -0,0 +1,57 @@ +BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC +! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/beta double excitations +! +! from contraction with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_beta a_h2_beta a_h1_alpha + END_DOC + integer :: i,h1,p1,h2,p2 + integer :: hh1,hh2,pp1,pp2,m,mm + integer :: Ne(2) + integer, allocatable :: occ(:,:) + double precision :: contrib + allocate( occ(N_int*bit_kind_size,2) ) + call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + eff_2_e_from_3_e_ab = 0.d0 + do hh1 = 1, n_act_orb !! alpha + h1 = list_act(hh1) + do hh2 = 1, n_act_orb !! beta + h2 = list_act(hh2) + do pp1 = 1, n_act_orb !! alpha + p1 = list_act(pp1) + do pp2 = 1, n_act_orb !! beta + p2 = list_act(pp2) + call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) + eff_2_e_from_3_e_ab(p2,p1,h2,h1) = contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + +subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) + implicit none + integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) + double precision, intent(out) :: contrib + integer :: mm,m + double precision :: direct_int, exchange_int + !! h1,p1 == alpha + !! h2,p2 == beta + contrib = 0.d0 + do mm = 1, Ne(1) !! alpha + m = occ(m,1) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + ! exchange between (h1,p1) and m + exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1) + contrib += direct_int - exchange_int + enddo + + do mm = 1, Ne(2) !! beta + m = occ(m,2) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + ! exchange between (h2,p2) and m + exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1) + contrib += direct_int - exchange_int + enddo +end diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index 8bdc57ee..f3641049 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -10,6 +10,7 @@ program test_normal_order read_wf = .True. touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call provide_all_three_ints_bi_ortho call test end @@ -28,7 +29,7 @@ subroutine test s2 = 2 accu = 0.d0 do h1 = 1, elec_beta_num - do p1 = elec_beta_num+1, mo_num + do p1 = elec_alpha_num+1, mo_num do h2 = 1, elec_beta_num do p2 = elec_beta_num+1, mo_num det_i = ref_bitmask @@ -38,36 +39,37 @@ subroutine test call get_excitation_degree(ref_bitmask,det_i,degree,N_int) call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) hthree *= phase - normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) +! normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) + normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) accu += dabs(hthree-normal) enddo enddo enddo enddo - print*,'accu opposite spin = ',accu +print*,'accu opposite spin = ',accu - s1 = 2 - s2 = 2 - accu = 0.d0 - do h1 = 1, elec_beta_num - do p1 = elec_beta_num+1, mo_num - do h2 = h1+1, elec_beta_num - do p2 = elec_beta_num+1, mo_num - det_i = ref_bitmask - call do_single_excitation(det_i,h1,p1,s1,i_ok) - call do_single_excitation(det_i,h2,p2,s2,i_ok) - if(i_ok.ne.1)cycle - call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) - call get_excitation_degree(ref_bitmask,det_i,degree,N_int) - call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) - hthree *= phase - normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) - accu += dabs(hthree-normal) - enddo - enddo - enddo - enddo - print*,'accu same spin = ',accu +!s1 = 2 +!s2 = 2 +!accu = 0.d0 +!do h1 = 1, elec_beta_num +! do p1 = elec_beta_num+1, mo_num +! do h2 = h1+1, elec_beta_num +! do p2 = elec_beta_num+1, mo_num +! det_i = ref_bitmask +! call do_single_excitation(det_i,h1,p1,s1,i_ok) +! call do_single_excitation(det_i,h2,p2,s2,i_ok) +! if(i_ok.ne.1)cycle +! call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) +! call get_excitation_degree(ref_bitmask,det_i,degree,N_int) +! call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) +! hthree *= phase +! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) +! accu += dabs(hthree-normal) +! enddo +! enddo +! enddo +!enddo +!print*,'accu same spin = ',accu end From ac2ebda9ce4f3efb672dbcd9b989f82f9719dc9a Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 15:49:39 +0100 Subject: [PATCH 09/14] beginning to work on double exc with optimization --- src/tc_bi_ortho/slater_tc_3e.irp.f | 24 ++-- src/tc_bi_ortho/slater_tc_opt.irp.f | 50 ++++++++ src/tc_bi_ortho/slater_tc_opt_double.irp.f | 141 ++++++++++++++++++++- src/tc_bi_ortho/test_normal_order.irp.f | 98 ++++++++++---- 4 files changed, 275 insertions(+), 38 deletions(-) create mode 100644 src/tc_bi_ortho/slater_tc_opt.irp.f diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index 0d5f8542..9740ee2f 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -256,20 +256,16 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) if(Ne(1)+Ne(2).ge.3)then if(s1==s2)then ! same spin excitation ispin = other_spin(s1) -! print*,'htilde ij' - do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1) - mm = occ(m,ispin) -!! direct_int = three_body_ints_bi_ort(mm,p2,p1,mm,h2,h1) -!! exchange_int = three_body_ints_bi_ort(mm,p2,p1,mm,h1,h2) - direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) - exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) -! print*,direct_int,exchange_int - hthree += direct_int - exchange_int - enddo - do m = 1, Ne(s1) ! pure contribution from s1 - mm = occ(m,s1) - hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1) - enddo + do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1) + mm = occ(m,ispin) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) + hthree += direct_int - exchange_int + enddo + do m = 1, Ne(s1) ! pure contribution from s1 + mm = occ(m,s1) + hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1) + enddo else ! different spin excitation do m = 1, Ne(s1) mm = occ(m,s1) ! diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f new file mode 100644 index 00000000..c334b274 --- /dev/null +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -0,0 +1,50 @@ +subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) + + BEGIN_DOC + ! + ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis + !! + ! Returns the detail of the matrix element in terms of single, two and three electron contribution. + !! WARNING !! + ! + ! Non hermitian !! + ! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: degree + + hmono = 0.d0 + htwoe = 0.d0 + htot = 0.d0 + hthree = 0.D0 + + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.gt.2) return + + if(degree == 0)then + call diag_htilde_mu_mat_fock_bi_ortho(Nint, key_i, hmono, htwoe, hthree, htot) + else if (degree == 1)then + call single_htilde_mu_mat_fock_bi_ortho (Nint,key_j, key_i , hmono, htwoe, hthree, htot) + else if(degree == 2)then + call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + if(three_body_h_tc) then + if(.not.double_normal_ord) then + call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + endif + endif + endif + + htot = hmono + htwoe + hthree + if(degree==0) then + htot += nuclear_repulsion + endif + +end + +! --- diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index 9cff8ff3..ef319c47 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -32,6 +32,11 @@ END_PROVIDER subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) implicit none + BEGIN_DOC +! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_beta +! +! on top of a determinant whose occupied orbitals is in (occ, Ne) + END_DOC integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) double precision, intent(out) :: contrib integer :: mm,m @@ -40,7 +45,7 @@ subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) !! h2,p2 == beta contrib = 0.d0 do mm = 1, Ne(1) !! alpha - m = occ(m,1) + m = occ(mm,1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) ! exchange between (h1,p1) and m exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1) @@ -48,10 +53,142 @@ subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) enddo do mm = 1, Ne(2) !! beta - m = occ(m,2) + m = occ(mm,2) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) ! exchange between (h2,p2) and m exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1) contrib += direct_int - exchange_int enddo end + +BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC +! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/alpha double excitations +! +! from contractionelec_alpha_num with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_alpha a_h2_alpha a_h1_alpha +! +! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill +! +! |||| h2>h1, p2>p1 |||| + END_DOC + integer :: i,h1,p1,h2,p2 + integer :: hh1,hh2,pp1,pp2,m,mm + integer :: Ne(2) + integer, allocatable :: occ(:,:) + double precision :: contrib + allocate( occ(N_int*bit_kind_size,2) ) + call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + eff_2_e_from_3_e_aa = 100000000.d0 + do hh1 = 1, n_act_orb !! alpha + h1 = list_act(hh1) + do hh2 = hh1+1, n_act_orb !! alpha + h2 = list_act(hh2) + do pp1 = 1, n_act_orb !! alpha + p1 = list_act(pp1) + do pp2 = pp1+1, n_act_orb !! alpha + p2 = list_act(pp2) + call give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib) + eff_2_e_from_3_e_aa(p2,p1,h2,h1) = contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + +subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib) + implicit none + BEGIN_DOC +! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_alpha +! +! on top of a determinant whose occupied orbitals is in (occ, Ne) + END_DOC + integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) + double precision, intent(out) :: contrib + integer :: mm,m + double precision :: direct_int, exchange_int + double precision :: three_e_double_parrallel_spin + !! h1,p1 == alpha + !! h2,p2 == alpha + contrib = 0.d0 + do mm = 1, Ne(1) !! alpha ==> pure parallele spin contribution + m = occ(mm,1) + contrib += three_e_double_parrallel_spin(m,p2,h2,p1,h1) + enddo + + do mm = 1, Ne(2) !! beta + m = occ(mm,2) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + ! exchange between (h1,p1) and (h2,p2) + exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) + contrib += direct_int - exchange_int + enddo +end + + +BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, mo_num)] + implicit none + BEGIN_DOC +! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for beta/beta double excitations +! +! from contractionelec_beta_num with HF density = a^{dagger}_p1_beta a^{dagger}_p2_beta a_h2_beta a_h1_beta +! +! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill +! +! |||| h2>h1, p2>p1 |||| + END_DOC + integer :: i,h1,p1,h2,p2 + integer :: hh1,hh2,pp1,pp2,m,mm + integer :: Ne(2) + integer, allocatable :: occ(:,:) + double precision :: contrib + allocate( occ(N_int*bit_kind_size,2) ) + call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + eff_2_e_from_3_e_bb = 100000000.d0 + do hh1 = 1, n_act_orb !! beta + h1 = list_act(hh1) + do hh2 = hh1+1, n_act_orb !! beta + h2 = list_act(hh2) + do pp1 = 1, n_act_orb !! beta + p1 = list_act(pp1) + do pp2 = pp1+1, n_act_orb !! beta + p2 = list_act(pp2) + call give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib) + eff_2_e_from_3_e_bb(p2,p1,h2,h1) = contrib + enddo + enddo + enddo + enddo + +END_PROVIDER + +subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib) + implicit none + BEGIN_DOC +! gives the contribution for a double excitation (h1,p1)_beta (h2,p2)_beta +! +! on top of a determinant whose occupied orbitals is in (occ, Ne) + END_DOC + integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) + double precision, intent(out) :: contrib + integer :: mm,m + double precision :: direct_int, exchange_int + double precision :: three_e_double_parrallel_spin + !! h1,p1 == beta + !! h2,p2 == beta + contrib = 0.d0 + do mm = 1, Ne(2) !! beta ==> pure parallele spin contribution + m = occ(mm,1) + contrib += three_e_double_parrallel_spin(m,p2,h2,p1,h1) + enddo + + do mm = 1, Ne(1) !! alpha + m = occ(mm,1) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + ! exchange between (h1,p1) and (h2,p2) + exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) + contrib += direct_int - exchange_int + enddo +end + diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index f3641049..46705f5f 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -39,7 +39,7 @@ subroutine test call get_excitation_degree(ref_bitmask,det_i,degree,N_int) call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) hthree *= phase -! normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) +! !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) accu += dabs(hthree-normal) enddo @@ -48,28 +48,82 @@ subroutine test enddo print*,'accu opposite spin = ',accu -!s1 = 2 -!s2 = 2 -!accu = 0.d0 -!do h1 = 1, elec_beta_num -! do p1 = elec_beta_num+1, mo_num -! do h2 = h1+1, elec_beta_num -! do p2 = elec_beta_num+1, mo_num -! det_i = ref_bitmask -! call do_single_excitation(det_i,h1,p1,s1,i_ok) -! call do_single_excitation(det_i,h2,p2,s2,i_ok) -! if(i_ok.ne.1)cycle -! call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) -! call get_excitation_degree(ref_bitmask,det_i,degree,N_int) -! call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) -! hthree *= phase +! p2=6 +! p1=5 +! h2=2 +! h1=1 + +s1 = 1 +s2 = 1 +accu = 0.d0 +do h1 = 1, elec_alpha_num + do p1 = elec_alpha_num+1, mo_num + do p2 = p1+1, mo_num + do h2 = h1+1, elec_alpha_num + det_i = ref_bitmask + call do_single_excitation(det_i,h1,p1,s1,i_ok) + if(i_ok.ne.1)cycle + call do_single_excitation(det_i,h2,p2,s2,i_ok) + if(i_ok.ne.1)cycle + call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call get_excitation_degree(ref_bitmask,det_i,degree,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + integer :: hh1, pp1, hh2, pp2, ss1, ss2 + call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2) + hthree *= phase ! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) -! accu += dabs(hthree-normal) -! enddo -! enddo -! enddo -!enddo -!print*,'accu same spin = ',accu + normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1) + if(dabs(hthree).lt.1.d-10)cycle + if(dabs(hthree-normal).gt.1.d-10)then + print*,pp2,pp1,hh2,hh1 + print*,p2,p1,h2,h1 + print*,hthree,normal,dabs(hthree-normal) + stop + endif +! print*,hthree,normal,dabs(hthree-normal) + accu += dabs(hthree-normal) + enddo + enddo + enddo +enddo +print*,'accu same spin alpha = ',accu + + +s1 = 2 +s2 = 2 +accu = 0.d0 +do h1 = 1, elec_beta_num + do p1 = elec_beta_num+1, mo_num + do p2 = p1+1, mo_num + do h2 = h1+1, elec_beta_num + det_i = ref_bitmask + call do_single_excitation(det_i,h1,p1,s1,i_ok) + if(i_ok.ne.1)cycle + call do_single_excitation(det_i,h2,p2,s2,i_ok) + if(i_ok.ne.1)cycle + call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call get_excitation_degree(ref_bitmask,det_i,degree,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2) + hthree *= phase +! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) + normal = eff_2_e_from_3_e_bb(p2,p1,h2,h1) + if(dabs(hthree).lt.1.d-10)cycle + if(dabs(hthree-normal).gt.1.d-10)then + print*,pp2,pp1,hh2,hh1 + print*,p2,p1,h2,h1 + print*,hthree,normal,dabs(hthree-normal) + stop + endif +! print*,hthree,normal,dabs(hthree-normal) + accu += dabs(hthree-normal) + enddo + enddo + enddo +enddo +print*,'accu same spin beta = ',accu + + end From 9eba8d692d549e24d30e3236d432e800b7551249 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 16:33:37 +0100 Subject: [PATCH 10/14] optimized all matrix elements with three body terms --- src/tc_bi_ortho/slater_tc_opt.irp.f | 12 +- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 197 +++++++++++++++++++++ src/tc_bi_ortho/test_normal_order.irp.f | 4 +- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 3 +- 4 files changed, 205 insertions(+), 11 deletions(-) diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index c334b274..8ab3388c 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -28,19 +28,13 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree.gt.2) return if(degree == 0)then - call diag_htilde_mu_mat_fock_bi_ortho(Nint, key_i, hmono, htwoe, hthree, htot) + call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1)then - call single_htilde_mu_mat_fock_bi_ortho (Nint,key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint,key_j, key_i , hmono, htwoe, hthree, htot) else if(degree == 2)then - call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) - if(three_body_h_tc) then - if(.not.double_normal_ord) then - call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) - endif - endif + call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) endif - htot = hmono + htwoe + hthree if(degree==0) then htot += nuclear_repulsion endif diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index ef319c47..ca1d0eea 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -1,3 +1,200 @@ + +subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) + + BEGIN_DOC + ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int,phase + + + call get_excitation_degree(key_i, key_j, degree, Nint) + + hmono = 0.d0 + htwoe = 0.d0 + hthree = 0.d0 + htot = 0.d0 + + if(degree.ne.2)then + return + endif + + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + call get_double_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + + if(s1.ne.s2)then + ! opposite spin two-body +! key_j, key_i + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + if(three_body_h_tc)then + if(.not.double_normal_ord)then + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + elseif(double_normal_ord.and.+Ne(1).gt.2)then + htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? + endif + endif + else + ! same spin two-body + ! direct terms + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + ! exchange terms + htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) + if(three_body_h_tc)then + if(.not.double_normal_ord)then + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + elseif(double_normal_ord.and.+Ne(1).gt.2)then + htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? + htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? + endif + endif + endif + hthree *= phase + htwoe *= phase + htot = htwoe + hthree + +end + + + +subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + implicit none + integer(bit_kind), intent(in) :: key_i(N_int,2) + integer, intent(in) :: h1,h2,p1,p2,s1,s2 + double precision, intent(out) :: hthree + integer :: nexc(2),i,ispin,na,nb + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: particle(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_particle(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) + integer(bit_kind) :: det_tmp(N_int,2) + integer :: ipart, ihole + double precision :: direct_int, exchange_int, three_e_double_parrallel_spin + + nexc(1) = 0 + nexc(2) = 0 + !! Get all the holes and particles of key_i with respect to the ROHF determinant + do i=1,N_int + hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),key_i(i,1)) + particle(i,2) = iand(hole(i,2),key_i(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta + if(s1==s2.and.s1==1)then + !!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc + hthree = eff_2_e_from_3_e_aa(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 1 ! i==alpha ==> pure same spin terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ipart=occ_particle(i,ispin) + hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + ihole=occ_hole(i,ispin) + hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + enddo + ispin = 2 ! i==beta ==> alpha/alpha/beta terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and (h2,p2) + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + elseif(s1==s2.and.s1==2)then + !!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc + hthree = eff_2_e_from_3_e_bb(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 2 ! i==beta ==> pure same spin terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ipart=occ_particle(i,ispin) + hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + ihole=occ_hole(i,ispin) + hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + enddo + ispin = 1 ! i==alpha==> beta/beta/alpha terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and (h2,p2) + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + else ! (h1,p1) == alpha/(h2,p2) == beta + hthree = eff_2_e_from_3_e_ab(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 1 ! i==alpha ==> alpha/beta/alpha terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and i + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch13_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch13_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + ispin = 2 ! i==beta ==> alpha/beta/beta terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h2,p2) and i + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch23_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch23_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + endif +end + + BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)] implicit none BEGIN_DOC diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index 46705f5f..118e481a 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -40,13 +40,15 @@ subroutine test call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) hthree *= phase ! !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) - normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) + call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal) +! normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) accu += dabs(hthree-normal) enddo enddo enddo enddo print*,'accu opposite spin = ',accu +stop ! p2=6 ! p1=5 diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 8afb3b25..551eced2 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -50,13 +50,14 @@ subroutine test_slater_tc_opt accu += dabs(htot-hnewtot) ! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i + print*,j,i,degree call debug_det(psi_det(1,1,i),N_int) call debug_det(psi_det(1,1,j),N_int) print*,htot,hnewtot,dabs(htot-hnewtot) ! print*,hthree,hnewthree,dabs(hthree-hnewthree) stop endif + print*,htot,hnewtot,dabs(htot-hnewtot) endif enddo enddo From a5ded6cd59b1c30695836cdca07ae17834972140 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 17:30:08 +0100 Subject: [PATCH 11/14] added providers for the totally symmetrized integrals --- src/tc_bi_ortho/slater_tc_opt_diag.irp.f | 8 +- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 16 +- src/tc_bi_ortho/slater_tc_opt_single.irp.f | 8 +- .../symmetrized_3_e_int_prov.irp.f | 140 ++++++++++++++++++ src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 2 +- 5 files changed, 156 insertions(+), 18 deletions(-) create mode 100644 src/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f diff --git a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f index 4048f481..c0b59969 100644 --- a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -113,7 +113,7 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer :: occ(Nint*bit_kind_size,2) integer :: other_spin integer :: k,l,i,jj,mm,j,m - double precision :: three_e_diag_parrallel_spin, direct_int, exchange_int + double precision :: direct_int, exchange_int if (iorb < 1) then @@ -163,7 +163,7 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) jj = occ(j,ispin) do m = j+1, na mm = occ(m,ispin) - hthree += three_e_diag_parrallel_spin(mm,jj,iorb) + hthree += three_e_diag_parrallel_spin_prov(mm,jj,iorb) enddo enddo !! same-spin/oposite-spin @@ -210,7 +210,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer(bit_kind), intent(inout) :: key(Nint,2) double precision, intent(inout) :: hmono,htwoe,hthree - double precision :: direct_int, exchange_int, three_e_diag_parrallel_spin + double precision :: direct_int, exchange_int integer :: occ(Nint*bit_kind_size,2) integer :: other_spin integer :: k,l,i,jj,mm,j,m @@ -250,7 +250,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) jj = occ(j,ispin) do m = j+1, na mm = occ(m,ispin) - hthree -= three_e_diag_parrallel_spin(mm,jj,iorb) + hthree -= three_e_diag_parrallel_spin_prov(mm,jj,iorb) enddo enddo !! same-spin/oposite-spin diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index ca1d0eea..c16c673d 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -84,7 +84,7 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) integer(bit_kind) :: det_tmp(N_int,2) integer :: ipart, ihole - double precision :: direct_int, exchange_int, three_e_double_parrallel_spin + double precision :: direct_int, exchange_int nexc(1) = 0 nexc(2) = 0 @@ -118,9 +118,9 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) ispin = 1 ! i==alpha ==> pure same spin terms do i = 1, nexc(ispin) ! number of couple of holes/particles ipart=occ_particle(i,ispin) - hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1) ihole=occ_hole(i,ispin) - hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1) enddo ispin = 2 ! i==beta ==> alpha/alpha/beta terms do i = 1, nexc(ispin) ! number of couple of holes/particles @@ -145,9 +145,9 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) ispin = 2 ! i==beta ==> pure same spin terms do i = 1, nexc(ispin) ! number of couple of holes/particles ipart=occ_particle(i,ispin) - hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1) ihole=occ_hole(i,ispin) - hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1) enddo ispin = 1 ! i==alpha==> beta/beta/alpha terms do i = 1, nexc(ispin) ! number of couple of holes/particles @@ -305,13 +305,12 @@ subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib) double precision, intent(out) :: contrib integer :: mm,m double precision :: direct_int, exchange_int - double precision :: three_e_double_parrallel_spin !! h1,p1 == alpha !! h2,p2 == alpha contrib = 0.d0 do mm = 1, Ne(1) !! alpha ==> pure parallele spin contribution m = occ(mm,1) - contrib += three_e_double_parrallel_spin(m,p2,h2,p1,h1) + contrib += three_e_double_parrallel_spin_prov(m,p2,h2,p1,h1) enddo do mm = 1, Ne(2) !! beta @@ -371,13 +370,12 @@ subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib) double precision, intent(out) :: contrib integer :: mm,m double precision :: direct_int, exchange_int - double precision :: three_e_double_parrallel_spin !! h1,p1 == beta !! h2,p2 == beta contrib = 0.d0 do mm = 1, Ne(2) !! beta ==> pure parallele spin contribution m = occ(mm,1) - contrib += three_e_double_parrallel_spin(m,p2,h2,p1,h1) + contrib += three_e_double_parrallel_spin_prov(m,p2,h2,p1,h1) enddo do mm = 1, Ne(1) !! alpha diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f index cb9306aa..ae41591a 100644 --- a/src/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -196,7 +196,7 @@ subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree, integer :: occ(Nint*bit_kind_size,2) integer :: other_spin integer :: k,l,i,jj,j - double precision :: three_e_single_parrallel_spin, direct_int, exchange_int + double precision :: direct_int, exchange_int if (iorb < 1) then @@ -236,7 +236,7 @@ subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree, !! jj = ispin = ispin_fock >> pure parallel spin do j = 1, na jj = occ(j,ispin) - hthree += three_e_single_parrallel_spin(jj,iorb,p_fock,h_fock) + hthree += three_e_single_parrallel_spin_prov(jj,iorb,p_fock,h_fock) enddo !! spin of jj == other spin than ispin AND ispin_fock !! exchange between the iorb and (h_fock, p_fock) @@ -287,7 +287,7 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N integer(bit_kind), intent(inout) :: key(Nint,2) double precision, intent(inout) :: hthree - double precision :: direct_int, exchange_int, three_e_single_parrallel_spin + double precision :: direct_int, exchange_int integer :: occ(Nint*bit_kind_size,2) integer :: other_spin integer :: k,l,i,jj,mm,j,m @@ -315,7 +315,7 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N !! jj = ispin = ispin_fock >> pure parallel spin do j = 1, na jj = occ(j,ispin) - hthree -= three_e_single_parrallel_spin(jj,iorb,p_fock,h_fock) + hthree -= three_e_single_parrallel_spin_prov(jj,iorb,p_fock,h_fock) enddo !! spin of jj == other spin than ispin AND ispin_fock !! exchange between the iorb and (h_fock, p_fock) diff --git a/src/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f b/src/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f new file mode 100644 index 00000000..e8277a74 --- /dev/null +++ b/src/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f @@ -0,0 +1,140 @@ + +BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS + ! + ! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form for same spin matrix elements + ! + ! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin + + three_e_diag_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_diag_parrallel_spin_prov ...' + + integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,m,integral) & + !$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do j = 1, mo_num + do m = j, mo_num + three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do i = 1, mo_num + do j = 1, mo_num + do m = 1, j + three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of for same spin matrix elements + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0, three_e_single_parrallel_spin + + three_e_single_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_single_parrallel_spin_prov ...' + + integral = three_e_single_parrallel_spin(1,1,1,1) + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do m = 1, mo_num + three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + + +! --- + +BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_double_parrallel_spin_prov(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0, three_e_double_parrallel_spin + + three_e_double_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_double_parrallel_spin_prov ...' + call wall_time(wall0) + + integral = three_e_double_parrallel_spin(1,1,1,1,1) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i) + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 551eced2..7d063c61 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -57,7 +57,7 @@ subroutine test_slater_tc_opt ! print*,hthree,hnewthree,dabs(hthree-hnewthree) stop endif - print*,htot,hnewtot,dabs(htot-hnewtot) +! print*,htot,hnewtot,dabs(htot-hnewtot) endif enddo enddo From d0fecfa84577d3f9eee07615bb4399ad33eebe69 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 18:14:29 +0100 Subject: [PATCH 12/14] parallelized the two electron terms for opt doubles tc --- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 24 ++++++++++++ src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 43 +++++++++++++++++++++- 2 files changed, 66 insertions(+), 1 deletion(-) diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index c16c673d..bd2d37a3 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -209,7 +209,13 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, double precision :: contrib allocate( occ(N_int*bit_kind_size,2) ) call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + call give_contrib_for_abab(1,1,1,1,occ,Ne,contrib) eff_2_e_from_3_e_ab = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & + !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_ab) + !$OMP DO SCHEDULE (static) do hh1 = 1, n_act_orb !! alpha h1 = list_act(hh1) do hh2 = 1, n_act_orb !! beta @@ -224,6 +230,8 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, enddo enddo enddo + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER @@ -276,7 +284,13 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, double precision :: contrib allocate( occ(N_int*bit_kind_size,2) ) call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + call give_contrib_for_aaaa(1 ,1 ,1 ,1 ,occ,Ne,contrib) eff_2_e_from_3_e_aa = 100000000.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & + !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_aa) + !$OMP DO SCHEDULE (static) do hh1 = 1, n_act_orb !! alpha h1 = list_act(hh1) do hh2 = hh1+1, n_act_orb !! alpha @@ -291,6 +305,8 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, enddo enddo enddo + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER @@ -341,7 +357,13 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, double precision :: contrib allocate( occ(N_int*bit_kind_size,2) ) call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + call give_contrib_for_bbbb(1,1 ,1 ,1 ,occ,Ne,contrib) eff_2_e_from_3_e_bb = 100000000.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & + !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_bb) + !$OMP DO SCHEDULE (static) do hh1 = 1, n_act_orb !! beta h1 = list_act(hh1) do hh2 = hh1+1, n_act_orb !! beta @@ -356,6 +378,8 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, enddo enddo enddo + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 7d063c61..66ca2e6a 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,7 +11,8 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call test_slater_tc_opt +! call test_slater_tc_opt + call timing_hij end subroutine test_slater_tc_opt @@ -65,3 +66,43 @@ subroutine test_slater_tc_opt print*,'accu = ',accu/i_count end + + +subroutine timing_hij + implicit none + integer :: i,j + double precision :: wall0, wall1 + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot +! allocate(mat_old(N_det,N_det)) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall0) + do i = 1, N_det + do j = 1, N_det + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) +! mat_old(j,i) = htot + enddo + enddo + call wall_time(wall1) + print*,'time for old hij = ',wall1 - wall0 + +! allocate(mat_new(N_det,N_det)) + call wall_time(wall0) + do i = 1, N_det + do j = 1, N_det + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) +! mat_new(j,i) = htot + enddo + enddo + call wall_time(wall1) + print*,'time for new hij = ',wall1 - wall0 + double precision :: accu + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det +! accu += dabs(mat_new(j,i) - mat_old(j,i)) + enddo + enddo + print*,'accu = ',accu + +end From d07bbacd8c160ec1b3d43d0dc6529a26e0e886ed Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 21 Jan 2023 12:21:31 +0100 Subject: [PATCH 13/14] minor modifs in test_tc_bi_ortho.irp.f --- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 35 +++----------------------- 1 file changed, 3 insertions(+), 32 deletions(-) diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 66ca2e6a..1b247da5 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,7 +11,7 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid -! call test_slater_tc_opt + call test_slater_tc_opt call timing_hij end @@ -25,44 +25,15 @@ subroutine test_slater_tc_opt accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det -! do i = 1,1 - call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) - call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) - call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) -! print*,hthree,hnewthree -! print*,htot,hnewtot,dabs(hnewtot-htot) - accu_d += dabs(htot-hnewtot) - if(dabs(htot-hnewtot).gt.1.d-8)then - print*,i - print*,htot,hnewtot,dabs(htot-hnewtot) - endif -! do j = 319,319 do j = 1,N_det - if(i==j)cycle - integer :: degree - call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) -! if(degree .ne. 1)cycle call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot) -! if(dabs(hthree).gt.1.d-15)then if(dabs(htot).gt.1.d-15)then -! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - i_count += 1.D0 - accu += dabs(htot-hnewtot) -! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then - if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i,degree - call debug_det(psi_det(1,1,i),N_int) - call debug_det(psi_det(1,1,j),N_int) - print*,htot,hnewtot,dabs(htot-hnewtot) -! print*,hthree,hnewthree,dabs(hthree-hnewthree) - stop - endif -! print*,htot,hnewtot,dabs(htot-hnewtot) + i_count += 1.D0 + accu += dabs(htot-hnewtot) endif enddo enddo - print*,'accu_d = ',accu_d/N_det print*,'accu = ',accu/i_count end From 20da577c4f7b1302512e2c5739ce415bbe2cc082 Mon Sep 17 00:00:00 2001 From: eginer Date: Sun, 22 Jan 2023 17:00:55 +0100 Subject: [PATCH 14/14] added new timing test --- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 22 ++- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 204 ++++++++++++++++----- 2 files changed, 173 insertions(+), 53 deletions(-) diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index bd2d37a3..9d33523b 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -32,19 +32,23 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, if(degree.ne.2)then return endif - - call bitstring_to_list_ab(key_i, occ, Ne, Nint) + integer :: degree_i,degree_j + call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int) + call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int) call get_double_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) if(s1.ne.s2)then ! opposite spin two-body -! key_j, key_i htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) if(three_body_h_tc)then if(.not.double_normal_ord)then - call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) - elseif(double_normal_ord.and.+Ne(1).gt.2)then + if(degree_i>degree_j)then + call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) + else + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + endif + elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? endif endif @@ -56,8 +60,12 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) if(three_body_h_tc)then if(.not.double_normal_ord)then - call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) - elseif(double_normal_ord.and.+Ne(1).gt.2)then + if(degree_i>degree_j)then + call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) + else + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + endif + elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? endif diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 66ca2e6a..99352162 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -11,13 +11,16 @@ program tc_bi_ortho touch read_wf touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid -! call test_slater_tc_opt - call timing_hij +! call test_slater_tc_opt + call timing_tot +! call timing_diag +! call timing_single +! call timing_double end subroutine test_slater_tc_opt implicit none - integer :: i,j + integer :: i,j,degree double precision :: hmono, htwoe, htot, hthree double precision :: hnewmono, hnewtwoe, hnewthree, hnewtot double precision :: accu_d ,i_count, accu @@ -25,84 +28,193 @@ subroutine test_slater_tc_opt accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det -! do i = 1,1 - call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) - call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) - call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthree, hnewtot) -! print*,hthree,hnewthree -! print*,htot,hnewtot,dabs(hnewtot-htot) - accu_d += dabs(htot-hnewtot) - if(dabs(htot-hnewtot).gt.1.d-8)then - print*,i - print*,htot,hnewtot,dabs(htot-hnewtot) - endif -! do j = 319,319 do j = 1,N_det - if(i==j)cycle - integer :: degree - call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) -! if(degree .ne. 1)cycle call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot) -! if(dabs(hthree).gt.1.d-15)then if(dabs(htot).gt.1.d-15)then -! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then i_count += 1.D0 accu += dabs(htot-hnewtot) -! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then - if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i,degree - call debug_det(psi_det(1,1,i),N_int) - call debug_det(psi_det(1,1,j),N_int) - print*,htot,hnewtot,dabs(htot-hnewtot) -! print*,hthree,hnewthree,dabs(hthree-hnewthree) - stop - endif -! print*,htot,hnewtot,dabs(htot-hnewtot) + if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + print*,j,i,degree + call debug_det(psi_det(1,1,i),N_int) + call debug_det(psi_det(1,1,j),N_int) + print*,htot,hnewtot,dabs(htot-hnewtot) + print*,hthree,hnewthree,dabs(hthree-hnewthree) + stop + endif endif enddo enddo - print*,'accu_d = ',accu_d/N_det print*,'accu = ',accu/i_count end - -subroutine timing_hij +subroutine timing_tot implicit none integer :: i,j double precision :: wall0, wall1 double precision, allocatable :: mat_old(:,:),mat_new(:,:) - double precision :: hmono, htwoe, hthree, htot -! allocate(mat_old(N_det,N_det)) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det +! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + i_count += 1.d0 + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for old hij for total = ',wall1 - wall0 + + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det +! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + i_count += 1.d0 + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij for total = ',wall1 - wall0 + call i_H_j(psi_det(1,1,1), psi_det(1,1,2),N_int,htot) + call wall_time(wall0) + i_count = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_H_j(psi_det(1,1,j), psi_det(1,1,i),N_int,htot) + i_count += 1.d0 + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij STANDARD = ',wall1 - wall0 + +end + +subroutine timing_diag + implicit none + integer :: i,j + double precision :: wall0, wall1 + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) call wall_time(wall0) + i_count = 0.d0 do i = 1, N_det - do j = 1, N_det + do j = i,i + i_count += 1.d0 call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) -! mat_old(j,i) = htot enddo enddo call wall_time(wall1) - print*,'time for old hij = ',wall1 - wall0 + print*,'i_count = ',i_count + print*,'time for old hij for diagonal= ',wall1 - wall0 -! allocate(mat_new(N_det,N_det)) call wall_time(wall0) + i_count = 0.d0 do i = 1, N_det - do j = 1, N_det + do j = i,i + i_count += 1.d0 call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) -! mat_new(j,i) = htot enddo enddo call wall_time(wall1) - print*,'time for new hij = ',wall1 - wall0 - double precision :: accu + print*,'i_count = ',i_count + print*,'time for new hij for diagonal= ',wall1 - wall0 + +end + +subroutine timing_single + implicit none + integer :: i,j + double precision :: wall0, wall1,accu + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + i_count = 0.d0 accu = 0.d0 do i = 1, N_det do j = 1, N_det -! accu += dabs(mat_new(j,i) - mat_old(j,i)) + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.1)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 enddo enddo - print*,'accu = ',accu + print*,'i_count = ',i_count + print*,'time for old hij for singles = ',accu + + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.1)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + print*,'i_count = ',i_count + print*,'time for new hij for singles = ',accu end + +subroutine timing_double + implicit none + integer :: i,j + double precision :: wall0, wall1,accu + double precision, allocatable :: mat_old(:,:),mat_new(:,:) + double precision :: hmono, htwoe, hthree, htot, i_count + integer :: degree + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.2)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + print*,'i_count = ',i_count + print*,'time for old hij for doubles = ',accu + + i_count = 0.d0 + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) + if(degree.ne.2)cycle + i_count += 1.d0 + call wall_time(wall0) + call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + call wall_time(wall1) + accu += wall1 - wall0 + enddo + enddo + call wall_time(wall1) + print*,'i_count = ',i_count + print*,'time for new hij for doubles = ',accu + +end +