From a5ded6cd59b1c30695836cdca07ae17834972140 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 17:30:08 +0100 Subject: [PATCH] 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