From 10b602c50340372b662582b5b4009a1f19899ab0 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 3 Nov 2022 18:53:44 +0100 Subject: [PATCH] three body terms in TCSCF work for energy gradients but not for the print of energy --- src/tc_bi_ortho/slater_tc_3e.irp.f | 6 ------ src/tc_bi_ortho/test_tc_fock.irp.f | 5 +++-- src/tc_scf/fock_tc.irp.f | 20 +++++++++++++------- src/tc_scf/fock_tc_mo_tot.irp.f | 3 +++ src/tc_scf/fock_three.irp.f | 1 + src/tc_scf/fock_three_bi_ortho_new.irp.f | 11 ++++++----- 6 files changed, 26 insertions(+), 20 deletions(-) diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index f4899a80..0d5f8542 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -169,8 +169,6 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) ! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is ! < h1 j i | p1 j i > - < h1 j i | p1 i j > ! -! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor -! exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1) hthree += direct_int - exchange_int @@ -182,14 +180,10 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) ii = occ(i,ispin) ! other spin do j = 1, Ne(s1) ! same spin jj = occ(j,s1) ! same spin -! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor -! exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1) ! < h1 j i | p1 j i > - < h1 j i | j p1 i > hthree += direct_int - exchange_int -! print*,'h1,p1,ii,jj = ',h1,p1,ii,jj -! print*,direct_int, exchange_int enddo enddo ! diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index 7813b08f..d585ce6c 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -134,8 +134,8 @@ subroutine routine_3() 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 + do i = 1, elec_num_tab(s1) + do a = elec_num_tab(s1)+1, mo_num ! virtual det_i = ref_bitmask @@ -162,6 +162,7 @@ subroutine routine_3() print*, ' warning on', i, a print*, ref,new,err_ai endif + print*, ref,new,err_ai err_tot += err_ai write(22, *) htilde_ij diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 7eba9231..ca500fbb 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -76,13 +76,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)] END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis - END_DOC - Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta) -END_PROVIDER +!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] +! implicit none +! BEGIN_DOC +! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis +! END_DOC +! Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta) +!END_PROVIDER ! --- @@ -94,6 +94,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] if(bi_ortho)then call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) + if(three_body_h_tc)then + Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + endif else call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) @@ -110,6 +113,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] if(bi_ortho)then call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) + if(three_body_h_tc)then + Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + endif else call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f index 5eaa4726..a99c7698 100644 --- a/src/tc_scf/fock_tc_mo_tot.irp.f +++ b/src/tc_scf/fock_tc_mo_tot.irp.f @@ -113,6 +113,9 @@ enddo enddo endif + if(.not.bi_ortho .and. three_body_h_tc)then + Fock_matrix_tc_mo_tot += fock_3_mat + endif END_PROVIDER diff --git a/src/tc_scf/fock_three.irp.f b/src/tc_scf/fock_three.irp.f index e4348892..e9ad4ce5 100644 --- a/src/tc_scf/fock_three.irp.f +++ b/src/tc_scf/fock_three.irp.f @@ -16,6 +16,7 @@ BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)] fock_3_mat(j,i) = -contrib enddo enddo + else if(bi_ortho.and.three_body_h_tc)then !! !$OMP END DO !! !$OMP END PARALLEL !! do i = 1, mo_num diff --git a/src/tc_scf/fock_three_bi_ortho_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new.irp.f index b705a05c..004a2aa4 100644 --- a/src/tc_scf/fock_three_bi_ortho_new.irp.f +++ b/src/tc_scf/fock_three_bi_ortho_new.irp.f @@ -33,8 +33,9 @@ BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth, (mo_num, mo_num)] do i = 1, mo_num do a = 1, mo_num - do j = 1, elec_beta_num - do k = 1, elec_alpha_num + do j = 1, elec_alpha_num ! a + do k = 1, elec_beta_num ! b + ! a b a a b a 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 @@ -82,7 +83,7 @@ BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] 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 + fock_a_tot_3e_bi_orth = fock_a_abb_3e_bi_orth + fock_a_aba_3e_bi_orth + fock_a_aaa_3e_bi_orth END_PROVIDER @@ -121,8 +122,8 @@ BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth, (mo_num, mo_num)] do i = 1, mo_num do a = 1, mo_num - do j = 1, elec_alpha_num - do k = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_alpha_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