From 67c08c4089fce386ab1f6526db0b34094164d8ef Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 3 Nov 2022 18:14:08 +0100 Subject: [PATCH] naive version of three e fock matrix works for closed shell --- src/tc_bi_ortho/test_tc_fock.irp.f | 71 +++++---- src/tc_scf/fock_three_bi_ortho_new.irp.f | 176 +++++++++++++++++++++++ 2 files changed, 216 insertions(+), 31 deletions(-) create mode 100644 src/tc_scf/fock_three_bi_ortho_new.irp.f diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index badc40b6..7813b08f 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -120,44 +120,53 @@ subroutine routine_3() implicit none integer :: i, a, i_ok, s1 double precision :: hmono, htwoe, hthree, htilde_ij - double precision :: err_ai, err_tot + double precision :: err_ai, err_tot, ref, new integer(bit_kind), allocatable :: det_i(:,:) allocate(det_i(N_int,2)) err_tot = 0.d0 - s1 = 1 + do s1 = 1, 2 - det_i = ref_bitmask - call debug_det(det_i, N_int) - print*, ' HF det' - call debug_det(det_i, N_int) - - do i = 1, elec_alpha_num ! occupied - do a = elec_alpha_num+1, mo_num ! virtual - - - det_i = ref_bitmask - call do_single_excitation(det_i, i, a, s1, i_ok) - if(i_ok == -1) then - print*, 'PB !!' - print*, i, a - stop - endif - !print*, ' excited det' - !call debug_det(det_i, N_int) - - call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) - err_ai = dabs(htilde_ij) - if(err_ai .gt. 1d-7) then - print*, ' warning on', i, a - print*, hmono, htwoe, htilde_ij - endif - err_tot += err_ai - - write(22, *) htilde_ij - enddo + det_i = ref_bitmask + call debug_det(det_i, N_int) + print*, ' HF det' + call debug_det(det_i, N_int) + + do i = 1, elec_alpha_num ! occupied + do a = elec_alpha_num+1, mo_num ! virtual + + + det_i = ref_bitmask + call do_single_excitation(det_i, i, a, s1, i_ok) + if(i_ok == -1) then + print*, 'PB !!' + print*, i, a + stop + endif + !print*, ' excited det' + !call debug_det(det_i, N_int) + + call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) + if(dabs(hthree).lt.1.d-10)cycle + ref = hthree + if(s1 == 1)then + new = fock_a_tot_3e_bi_orth(a,i) + else if(s1 == 2)then + new = fock_b_tot_3e_bi_orth(a,i) + endif + err_ai = dabs(dabs(ref) - dabs(new)) + if(err_ai .gt. 1d-7) then + print*,'s1 = ',s1 + print*, ' warning on', i, a + print*, ref,new,err_ai + endif + err_tot += err_ai + + write(22, *) htilde_ij + enddo + enddo enddo print *, ' err_tot = ', err_tot diff --git a/src/tc_scf/fock_three_bi_ortho_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new.irp.f new file mode 100644 index 00000000..b705a05c --- /dev/null +++ b/src/tc_scf/fock_three_bi_ortho_new.irp.f @@ -0,0 +1,176 @@ +BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_abb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution + END_DOC + fock_a_abb_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_23_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + fock_a_abb_3e_bi_orth(a,i) += direct_int - exch_23_int + enddo + enddo + + enddo + enddo + fock_a_abb_3e_bi_orth = - fock_a_abb_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_aba_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution + END_DOC + fock_a_aba_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = 1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + fock_a_aba_3e_bi_orth(a,i) += direct_int - exch_13_int + enddo + enddo + + enddo + enddo + fock_a_aba_3e_bi_orth = - fock_a_aba_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_aaa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution + END_DOC + fock_a_aaa_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + ! positive terms :: cycle contrib + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + fock_a_aaa_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + fock_a_aaa_3e_bi_orth(a,i) += - exch_13_int - exch_23_int - exch_12_int + enddo + enddo + + enddo + enddo + fock_a_aaa_3e_bi_orth = - fock_a_aaa_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC + ! fock_a_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions + END_DOC + fock_a_tot_3e_bi_orth = fock_a_aaa_3e_bi_orth + fock_a_abb_3e_bi_orth + fock_a_aba_3e_bi_orth + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_baa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution + END_DOC + fock_b_baa_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_23_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + fock_b_baa_3e_bi_orth(a,i) += direct_int - exch_23_int + enddo + enddo + + enddo + enddo + fock_b_baa_3e_bi_orth = - fock_b_baa_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_bab_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution + END_DOC + fock_b_bab_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num + do k = 1, elec_beta_num + ! b a b b a b + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + fock_b_bab_3e_bi_orth(a,i) += direct_int - exch_13_int + enddo + enddo + + enddo + enddo + fock_b_bab_3e_bi_orth = - fock_b_bab_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_bbb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution + END_DOC + fock_b_bbb_3e_bi_orth = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + ! positive terms :: cycle contrib + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + fock_b_bbb_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + fock_b_bbb_3e_bi_orth(a,i) += - exch_13_int - exch_23_int - exch_12_int + enddo + enddo + + enddo + enddo + fock_b_bbb_3e_bi_orth = - fock_b_bbb_3e_bi_orth +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)] + implicit none + BEGIN_DOC + ! fock_b_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions + END_DOC + fock_b_tot_3e_bi_orth = fock_b_bbb_3e_bi_orth + fock_b_bab_3e_bi_orth + fock_b_baa_3e_bi_orth + +END_PROVIDER