routine htilde_mu_mat_opt_bi_ortho works

2024-12-30 15:15:38 +01:00 · 2023-07-10 12:10:53 +02:00 · 2023-07-10 12:10:53 +02:00 · 2d383d09c6
commit 2d383d09c6
parent 8c4a7226cd
5 changed files with 114 additions and 21 deletions
--- a/src/bi_ort_ints/three_body_ijmk.irp.f
+++ b/src/bi_ort_ints/three_body_ijmk.irp.f
@ -86,6 +86,13 @@
        tmp_loc_1 = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i)
        tmp_loc_2 = tmp_aux_2(ipoint,n)

+        tmp1(ipoint,1,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,1,k,i) * tmp_loc_2
+        tmp1(ipoint,2,n) = int2_grad1_u12_bimo_t(ipoint,2,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,2,k,i) * tmp_loc_2
+        tmp1(ipoint,3,n) = int2_grad1_u12_bimo_t(ipoint,3,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,3,k,i) * tmp_loc_2
+        tmp1(ipoint,4,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
+                         + int2_grad1_u12_bimo_t(ipoint,2,n,n) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
+                         + int2_grad1_u12_bimo_t(ipoint,3,n,n) * int2_grad1_u12_bimo_t(ipoint,3,k,i)
+
      enddo
    enddo
    !$OMP END DO
--- a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f
+++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f
@ -90,3 +90,96 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze)

 end 

+subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_chi_array,i_H_phi_array)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Computes $\langle i|H|Phi \rangle  = \sum_J c^R_J \langle i | H | J \rangle$.
+!
+! AND      $\langle Chi|H|  i \rangle  = \sum_J c^L_J \langle J | H | i \rangle$.
+!
+! CONVENTION: i_H_phi_array(0) =   total matrix element,  
+!
+!             i_H_phi_array(1) =   one-electron matrix element,
+!
+!             i_H_phi_array(2) =   two-electron matrix element,
+!
+!             i_H_phi_array(3) = three-electron matrix element,
+!
+! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$
+! is connected.
+!
+! The i_H_psi_minilist is much faster but requires to build the
+! minilists.
+  END_DOC
+  integer, intent(in)            :: Nint, Ndet,Ndet_max,Nstate
+  integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
+  integer(bit_kind), intent(in)  :: key(Nint,2)
+  double precision, intent(in)   :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate)
+  double precision, intent(out)  :: i_H_chi_array(0:3,Nstate),i_H_phi_array(0:3,Nstate)
+
+  integer                        :: i, ii,j
+  double precision               :: phase
+  integer                        :: exc(0:2,2,2)
+  double precision               :: hmono, htwoe, hthree, htot
+  integer, allocatable           :: idx(:)
+
+  ASSERT (Nint > 0)
+  ASSERT (N_int == Nint)
+  ASSERT (Nstate > 0)
+  ASSERT (Ndet > 0)
+  ASSERT (Ndet_max >= Ndet)
+  allocate(idx(0:Ndet))
+
+  i_H_chi_array = 0.d0
+  i_H_phi_array = 0.d0
+
+  call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
+  if (Nstate == 1) then
+
+    do ii=1,idx(0)
+      i = idx(ii)
+      ! computes <Chi|H_tc|i>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
+      i_H_chi_array(0,1) = i_H_chi_array(0,1) + coef_l(i,1)*htot
+      i_H_chi_array(1,1) = i_H_chi_array(1,1) + coef_l(i,1)*hmono
+      i_H_chi_array(2,1) = i_H_chi_array(2,1) + coef_l(i,1)*htwoe
+      i_H_chi_array(3,1) = i_H_chi_array(3,1) + coef_l(i,1)*hthree
+      ! computes <i|H_tc|Phi>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot)
+      i_H_phi_array(0,1) = i_H_phi_array(0,1) + coef_r(i,1)*htot
+      i_H_phi_array(1,1) = i_H_phi_array(1,1) + coef_r(i,1)*hmono
+      i_H_phi_array(2,1) = i_H_phi_array(2,1) + coef_r(i,1)*htwoe
+      i_H_phi_array(3,1) = i_H_phi_array(3,1) + coef_r(i,1)*hthree
+    enddo
+
+  else
+
+    do ii=1,idx(0)
+      i = idx(ii)
+      ! computes <Chi|H_tc|i>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot)
+      do j = 1, Nstate
+        i_H_chi_array(0,j) = i_H_chi_array(0,j) + coef_l(i,j)*htot
+        i_H_chi_array(1,j) = i_H_chi_array(1,j) + coef_l(i,j)*hmono
+        i_H_chi_array(2,j) = i_H_chi_array(2,j) + coef_l(i,j)*htwoe
+        i_H_chi_array(3,j) = i_H_chi_array(3,j) + coef_l(i,j)*hthree
+      enddo
+      ! computes <i|H_tc|Phi>
+      !DIR$ FORCEINLINE
+      call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot)
+      do j = 1, Nstate
+        i_H_phi_array(0,j) = i_H_phi_array(0,j) + coef_r(i,j)*htot
+        i_H_phi_array(1,j) = i_H_phi_array(1,j) + coef_r(i,j)*hmono
+        i_H_phi_array(2,j) = i_H_phi_array(2,j) + coef_r(i,j)*htwoe
+        i_H_phi_array(3,j) = i_H_phi_array(3,j) + coef_r(i,j)*hthree
+      enddo
+    enddo
+
+  endif
+
+end
+
--- a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f
+++ b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f
@ -184,7 +184,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
      ii = occ(i,s1)
      do j = i+1, Ne(s1)
       jj = occ(j,s1)
-!       ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
+!      !ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
       hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR
      enddo
     enddo
--- a/src/tc_bi_ortho/test_tc_fock.irp.f
+++ b/src/tc_bi_ortho/test_tc_fock.irp.f
@ -152,9 +152,7 @@ subroutine routine_tot()
 !   do i = 1, elec_num_tab(s1)
 !     do a = elec_num_tab(s1)+1, mo_num ! virtual 
   do i = 1, elec_beta_num
-     do a = elec_beta_num+1, elec_alpha_num! virtual 
-!   do i = elec_beta_num+1, elec_alpha_num
-!     do a = elec_alpha_num+1, mo_num! virtual 
+     do a = elec_beta_num+1, mo_num! virtual 
       print*,i,a 
 
       det_i = ref_bitmask
@ -167,7 +165,7 @@ subroutine routine_tot()
 
       call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
       print*,htilde_ij
-       if(dabs(htilde_ij).lt.1.d-10)cycle
+!       if(dabs(htilde_ij).lt.1.d-10)cycle
       print*, ' excited det'
       call debug_det(det_i, N_int)

@ -184,9 +182,12 @@ subroutine routine_tot()
 !       endif
       err_ai = dabs(dabs(ref) - dabs(new))
       if(err_ai .gt. 1d-7) then
+         print*,'---------'
         print*,'s1 = ',s1
         print*, ' warning on', i, a
         print*, ref,new,err_ai
+         print*,hmono, htwoe, hthree
+         print*,'---------'
       endif
       print*, ref,new,err_ai
       err_tot += err_ai
--- a/src/tc_scf/fock_tc.irp.f
+++ b/src/tc_scf/fock_tc.irp.f
@ -208,10 +208,10 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]

    if(three_body_h_tc) then
      !call wall_time(tt0)
-      !PROVIDE fock_a_tot_3e_bi_orth
-      !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
-      PROVIDE fock_3e_uhf_mo_a
-      Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
+      PROVIDE fock_a_tot_3e_bi_orth
+      Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
+!      PROVIDE fock_3e_uhf_mo_a
+!      Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a
      !call wall_time(tt1)
      !print*, ' 3-e term:', tt1-tt0
    endif
@ -241,21 +241,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]

  if(bi_ortho) then

-    !allocate(tmp(ao_num,ao_num))
-    !tmp = Fock_matrix_tc_ao_beta
-    !if(three_body_h_tc) then
-    !  tmp += fock_3e_uhf_ao_b
-    !endif
-    !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1))
-    !deallocate(tmp)
-
    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
-      !PROVIDE fock_b_tot_3e_bi_orth
-      !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
-      PROVIDE fock_3e_uhf_mo_b
-      Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
+       PROVIDE fock_b_tot_3e_bi_orth
+       Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
+!      PROVIDE fock_3e_uhf_mo_b
+!      Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b
    endif

  else