Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

plugins/local/bi_ort_ints/total_twoe_pot.irp.f

@@ -259,15 +259,21 @@ BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_transp, (mo_num, mo_num,
  END_DOC
 
  integer :: i,j,k,l
+ print*,'Providing mo_bi_ortho_tc_two_e_transp'
+ double precision :: t0,t1
+ call wall_time(t0)
  do i = 1, mo_num
   do j = 1, mo_num
    do k = 1, mo_num
     do l = 1, mo_num
-     mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e_transp(k,l,i,j)
+     mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e(k,l,i,j)
     enddo
    enddo
   enddo
  enddo
+ call wall_time(t1)
+ 
+ print *, ' WALL TIME for PROVIDING mo_bi_ortho_tc_two_e_transp (min)', (t1-t0)/60.d0
 
 END_PROVIDER 
 ! ---
@@ -326,3 +332,23 @@ END_PROVIDER
 
 ! ---
 
+ BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals_transp , (mo_num,mo_num,mo_num)]
+&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals_transp, (mo_num,mo_num,mo_num)]
+
+  BEGIN_DOC
+  ! tc_2e_3idx_coulomb_integrals_transp (j,k,i) = <jk|ji> 
+  ! tc_2e_3idx_exchange_integrals_transp(j,k,i) = <kj|ji> 
+  END_DOC
+ implicit none
+ integer :: i, j, k
+
+  do i = 1, mo_num
+    do k = 1, mo_num
+      do j = 1, mo_num
+        tc_2e_3idx_coulomb_integrals_transp(j, k,i) = mo_bi_ortho_tc_two_e_transp(j ,k ,j ,i ) 
+        tc_2e_3idx_exchange_integrals_transp(j,k,i) = mo_bi_ortho_tc_two_e_transp(k ,j ,j ,i ) 
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 

plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f

@@ -0,0 +1,108 @@
+subroutine get_d0_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
+  !todo: indices/conjg should be okay for complex
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
+  integer(bit_kind), intent(in) :: phasemask(N_int,2)
+  logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
+  integer(bit_kind) :: det(N_int, 2)
+  double precision, intent(in) :: coefs(N_states,2)
+  double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
+  double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+
+  integer :: i, j, k, s, h1, h2, p1, p2, puti, putj, mm
+  double precision :: phase
+  double precision :: hij,hji
+  double precision, external :: get_phase_bi
+  logical :: ok
+
+  integer, parameter :: bant=1
+  double precision, allocatable :: hij_cache1(:), hij_cache2(:)
+  allocate (hij_cache1(mo_num),hij_cache2(mo_num))
+  double precision, allocatable :: hji_cache1(:), hji_cache2(:)
+  allocate (hji_cache1(mo_num),hji_cache2(mo_num))
+!  print*,'in get_d0_new'
+!  call debug_det(gen,N_int)
+!  print*,'coefs',coefs(1,:)
+
+  if(sp == 3) then ! AB
+    h1 = p(1,1)
+    h2 = p(1,2)
+    do p1=1, mo_num
+      if(bannedOrb(p1, 1)) cycle
+!      call get_mo_two_e_integrals_complex(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
+      do mm = 1, mo_num
+       hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1)
+       hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,p1,h2,h1)
+      enddo
+      !!!!!!!!!! <alpha|H|psi>
+      do p2=1, mo_num
+        if(bannedOrb(p2,2)) cycle
+        if(banned(p1, p2, bant)) cycle ! rentable?
+        if(p1 == h1 .or. p2 == h2) then
+          call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
+          ! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
+!          call i_h_j_complex(det, gen, N_int, hij)
+          call htilde_mu_mat_opt_bi_ortho_no_3e_both(det,gen,N_int, hij,hji)
+        else
+          phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
+          hij = hij_cache1(p2) * phase
+          hji = hji_cache1(p2) * phase
+        end if
+        if (hij == 0.d0.or.hji == 0.d0) cycle
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij  ! HOTSPOT
+          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji  ! HOTSPOT
+        enddo
+      end do
+    end do
+
+  else ! AA BB
+    p1 = p(1,sp)
+    p2 = p(2,sp)
+    do puti=1, mo_num
+      if(bannedOrb(puti, sp)) cycle
+!      call get_mo_two_e_integrals_complex(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map,mo_integrals_map_2)
+!      call get_mo_two_e_integrals_complex(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map,mo_integrals_map_2)
+      do mm = 1, mo_num
+       hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,puti,p2,p1)
+       hij_cache2(mm) = mo_bi_ortho_tc_two_e(mm,puti,p1,p2)
+       hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p2,p1)
+       hji_cache2(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p1,p2)
+      enddo
+      !!!!!!!!!! <alpha|H|psi>
+      do putj=puti+1, mo_num
+        if(bannedOrb(putj, sp)) cycle
+        if(banned(puti, putj, bant)) cycle ! rentable?
+        if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
+          call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
+          !call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
+!          call i_h_j_complex(det, gen, N_int, hij)
+          call htilde_mu_mat_opt_bi_ortho_no_3e_both(det,gen,N_int, hij,hji)
+         if (hij == 0.d0.or.hji == 0.d0) cycle
+        else
+!          hij = (mo_two_e_integral_complex(p1, p2, puti, putj) -  mo_two_e_integral_complex(p2, p1, puti, putj))
+!          hij = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) -  mo_bi_ortho_tc_two_e(p2, p1, puti, putj))
+          hij = (mo_bi_ortho_tc_two_e(puti, putj, p1, p2) -  mo_bi_ortho_tc_two_e(puti, putj, p2, p1))
+          hji = (mo_bi_ortho_tc_two_e_transp(puti, putj, p1, p2) -  mo_bi_ortho_tc_two_e_transp(puti, putj, p2, p1))
+          if (hij == 0.d0.or.hji == 0.d0) cycle
+          phase = get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
+          hij = (hij) * phase
+          hji = (hji) * phase
+        end if
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
+          mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
+        enddo
+      end do
+
+    end do
+  end if
+
+  deallocate(hij_cache1,hij_cache2)
+end
+

plugins/local/cipsi_tc_bi_ortho/get_d1_transp.irp.f

@@ -0,0 +1,350 @@
+subroutine get_d1_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
+  !todo: indices should be okay for complex?
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in)  :: mask(N_int, 2), gen(N_int, 2)
+  integer(bit_kind), intent(in)  :: phasemask(N_int,2)
+  logical, intent(in)            :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
+  integer(bit_kind)              :: det(N_int, 2)
+  double precision, intent(in)   :: coefs(N_states,2)
+  double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
+  double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
+  integer, intent(in)            :: h(0:2,2), p(0:4,2), sp
+  double precision, external     :: get_phase_bi
+  double precision, external     :: mo_two_e_integral_complex
+  logical                        :: ok
+
+  logical, allocatable           :: lbanned(:,:)
+  integer                        :: puti, putj, ma, mi, s1, s2, i, i1, i2, j
+  integer                        :: hfix, pfix, h1, h2, p1, p2, ib, k, l, mm
+
+  integer, parameter             :: turn2(2) = (/2,1/)
+  integer, parameter             :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+
+  integer                        :: bant
+  double precision, allocatable :: hij_cache(:,:)
+  double precision               :: hij, tmp_rowij(N_states, mo_num), tmp_rowij2(N_states, mo_num),phase
+  double precision, allocatable :: hji_cache(:,:)
+  double precision               :: hji, tmp_rowji(N_states, mo_num), tmp_rowji2(N_states, mo_num)
+!  PROVIDE mo_integrals_map N_int
+!  print*,'in get_d1_new'
+!  call debug_det(gen,N_int)
+!  print*,'coefs',coefs(1,:)
+
+  allocate (lbanned(mo_num, 2))
+  allocate (hij_cache(mo_num,2))
+  allocate (hji_cache(mo_num,2))
+  lbanned = bannedOrb
+
+  do i=1, p(0,1)
+    lbanned(p(i,1), 1) = .true.
+  end do
+  do i=1, p(0,2)
+    lbanned(p(i,2), 2) = .true.
+  end do
+
+  ma = 1
+  if(p(0,2) >= 2) ma = 2
+  mi = turn2(ma)
+
+  bant = 1
+
+  if(sp == 3) then
+    !move MA
+    if(ma == 2) bant = 2
+    puti = p(1,mi)
+    hfix = h(1,ma)
+    p1 = p(1,ma)
+    p2 = p(2,ma)
+    if(.not. bannedOrb(puti, mi)) then
+!      call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
+!      call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
+      do mm = 1, mo_num
+       hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
+       hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
+       hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,p2)
+       hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,p1)
+      enddo
+      !! <alpha|H|psi>
+      tmp_rowij = 0.d0
+      tmp_rowji = 0.d0
+      do putj=1, hfix-1
+        if(lbanned(putj, ma)) cycle
+        if(banned(putj, puti,bant)) cycle
+        hij = hij_cache(putj,1) - hij_cache(putj,2)
+        hji = hji_cache(putj,1) - hji_cache(putj,2)
+        if (hij /= 0.d0.and.hji/=0.d0) then
+          phase = get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
+          hij = hij * phase
+          hji = hji * phase
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
+            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
+          enddo
+        endif
+      end do
+      do putj=hfix+1, mo_num
+        if(lbanned(putj, ma)) cycle
+        if(banned(putj, puti,bant)) cycle
+        hij = hij_cache(putj,2) - hij_cache(putj,1)
+        hji = hji_cache(putj,2) - hji_cache(putj,1)
+        if (hij /= 0.d0.and.hji/=0.d0) then
+          phase = get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
+          hij = hij * phase
+          hji = hji * phase
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
+            tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
+          enddo
+        endif
+      end do
+
+      if(ma == 1) then
+        mat_r(1:N_states,1:mo_num,puti) = mat_r(1:N_states,1:mo_num,puti) + tmp_rowij(1:N_states,1:mo_num)
+        mat_l(1:N_states,1:mo_num,puti) = mat_l(1:N_states,1:mo_num,puti) + tmp_rowji(1:N_states,1:mo_num)
+      else
+        do l=1,mo_num
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k,puti,l) = mat_r(k,puti,l) + tmp_rowij(k,l)
+            mat_l(k,puti,l) = mat_l(k,puti,l) + tmp_rowji(k,l)
+          enddo
+        enddo
+      end if
+
+    end if
+
+    !MOVE MI
+    pfix = p(1,mi)
+    tmp_rowij  = 0.d0
+    tmp_rowij2 = 0.d0
+    tmp_rowji  = 0.d0
+    tmp_rowji2 = 0.d0
+!    call get_mo_two_e_integrals_complex(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
+!    call get_mo_two_e_integrals_complex(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
+    do mm = 1, mo_num
+     hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1)
+     hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2)
+     hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,pfix,p1)
+     hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,pfix,p2)
+    enddo
+    putj = p1
+    !! <alpha|H|psi>
+    do puti=1,mo_num !HOT
+      if(lbanned(puti,mi)) cycle
+      !p1 fixed
+      putj = p1
+      if(.not. banned(putj,puti,bant)) then
+        hij = hij_cache(puti,2)
+        hji = hji_cache(puti,2)
+        if (hij /= 0.d0.and.hji/=0.d0) then
+          phase = get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
+          hij = hij * phase
+          hji = hji * phase
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
+            tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
+          enddo
+        endif
+      end if
+!      
+      putj = p2
+      if(.not. banned(putj,puti,bant)) then
+        hij = hij_cache(puti,1)
+        hji = hji_cache(puti,1)
+        if (hij /= 0.d0.and.hji/=0.d0) then
+          phase = get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
+          hij = hij * phase
+          hji = hji * phase
+          do k=1,N_states
+            tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
+            tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
+          enddo
+        endif
+      end if
+    end do
+
+    if(mi == 1) then
+      mat_r(:,:,p1) = mat_r(:,:,p1) + tmp_rowij(:,:)
+      mat_r(:,:,p2) = mat_r(:,:,p2) + tmp_rowij2(:,:)
+      mat_l(:,:,p1) = mat_l(:,:,p1) + tmp_rowji(:,:)
+      mat_l(:,:,p2) = mat_l(:,:,p2) + tmp_rowji2(:,:)
+    else
+      do l=1,mo_num
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij(k,l)
+          mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij2(k,l)
+          mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji(k,l)
+          mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji2(k,l)
+        enddo
+      enddo
+    end if
+
+  else  ! sp /= 3
+
+    if(p(0,ma) == 3) then
+      do i=1,3
+        hfix = h(1,ma)
+        puti = p(i, ma)
+        p1 = p(turn3(1,i), ma)
+        p2 = p(turn3(2,i), ma)
+!        call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
+!        call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
+        do mm = 1, mo_num
+         hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
+         hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
+         hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,p2)
+         hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,p1)
+        enddo
+    !! <alpha|H|psi>
+        tmp_rowij = 0.d0
+        tmp_rowji = 0.d0
+        do putj=1,hfix-1
+          if(banned(putj,puti,1)) cycle
+          if(lbanned(putj,ma)) cycle
+          hij = hij_cache(putj,1) - hij_cache(putj,2)
+          hji = hji_cache(putj,1) - hji_cache(putj,2)
+          if (hij /= 0.d0.and.hji/=0.d0) then
+            phase = get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
+            hij = hij * phase
+            hji = hji * phase
+            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
+            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
+          endif
+        end do
+        do putj=hfix+1,mo_num
+          if(banned(putj,puti,1)) cycle
+          if(lbanned(putj,ma)) cycle
+          hij = hij_cache(putj,2) - hij_cache(putj,1)
+          hji = hji_cache(putj,2) - hji_cache(putj,1)
+          if (hij /= 0.d0.and.hji/=0.d0) then
+            phase = get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
+            hij = hij * phase
+            hji = hji * phase
+            tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
+            tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
+          endif
+        end do
+
+        mat_r(:, :puti-1, puti) = mat_r(:, :puti-1, puti) + tmp_rowij(:,:puti-1)
+        mat_l(:, :puti-1, puti) = mat_l(:, :puti-1, puti) + tmp_rowji(:,:puti-1)
+        do l=puti,mo_num
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, puti, l) = mat_r(k, puti,l) + tmp_rowij(k,l)
+            mat_l(k, puti, l) = mat_l(k, puti,l) + tmp_rowji(k,l)
+          enddo
+        enddo
+      end do
+    else
+      hfix = h(1,mi)
+      pfix = p(1,mi)
+      p1 = p(1,ma)
+      p2 = p(2,ma)
+      tmp_rowij =  0.d0
+      tmp_rowij2 = 0.d0
+      tmp_rowji =  0.d0
+      tmp_rowji2 = 0.d0
+!      call get_mo_two_e_integrals_complex(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
+!      call get_mo_two_e_integrals_complex(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
+      do mm = 1, mo_num
+       hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,pfix)
+       hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,pfix)
+       hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,pfix)
+       hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,pfix)
+      enddo
+      putj = p2
+    !! <alpha|H|psi>
+      do puti=1,mo_num
+        if(lbanned(puti,ma)) cycle
+        putj = p2
+        if(.not. banned(puti,putj,1)) then
+          hij = hij_cache(puti,1)
+          hji = hji_cache(puti,1)
+          if (hij /= 0.d0.and.hji/=0.d0) then
+            phase = get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
+            hij = hij * phase
+            hji = hji * phase
+            !DIR$ LOOP COUNT AVG(4)
+            do k=1,N_states
+              tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
+              tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
+            enddo
+          endif
+        end if
+
+        putj = p1
+        if(.not. banned(puti,putj,1)) then
+          hij = hij_cache(puti,2)
+          hji = hji_cache(puti,2)
+          if (hij /= 0.d0.and.hji/=0.d0) then
+            phase = get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
+            hij = hij * phase
+            hji = hji * phase
+            do k=1,N_states
+              tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
+              tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
+            enddo
+          endif
+        end if
+      end do
+      mat_r(:,:p2-1,p2) = mat_r(:,:p2-1,p2) + tmp_rowij(:,:p2-1)
+      mat_l(:,:p2-1,p2) = mat_l(:,:p2-1,p2) + tmp_rowji(:,:p2-1)
+      do l=p2,mo_num
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij(k,l)
+          mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji(k,l)
+        enddo
+      enddo
+      mat_r(:,:p1-1,p1) = mat_r(:,:p1-1,p1) + tmp_rowij2(:,:p1-1)
+      mat_l(:,:p1-1,p1) = mat_l(:,:p1-1,p1) + tmp_rowji2(:,:p1-1)
+      do l=p1,mo_num
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij2(k,l)
+          mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji2(k,l)
+        enddo
+      enddo
+    end if
+  end if
+  deallocate(lbanned,hij_cache, hji_cache)
+
+ !! MONO
+    if(sp == 3) then
+      s1 = 1
+      s2 = 2
+    else
+      s1 = sp
+      s2 = sp
+    end if
+
+    do i1=1,p(0,s1)
+      ib = 1
+      if(s1 == s2) ib = i1+1
+      do i2=ib,p(0,s2)
+        p1 = p(i1,s1)
+        p2 = p(i2,s2)
+        if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
+        call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
+        ! gen is a selector; mask is ionized generator; det is alpha
+        ! hij is contribution to <psi|H|alpha>
+!        call i_h_j_complex(gen, det, N_int, hij)
+        call htilde_mu_mat_opt_bi_ortho_no_3e_both(det, gen, N_int, hij,hji)
+!        call htilde_mu_mat_opt_bi_ortho_no_3e(gen, det, N_int, hji)
+        !DIR$ LOOP COUNT AVG(4)
+        do k=1,N_states
+          ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * dconjg(hij)
+          mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij
+          mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji
+        enddo
+      end do
+    end do
+end
+

plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f

@@ -25,9 +25,6 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
 
   integer :: bant
   bant = 1
-!  print*, 'in get_d2_new'
-!  call debug_det(gen,N_int)
-!  print*,'coefs',coefs(1,:)
 
   tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles
 

plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f

@@ -0,0 +1,235 @@
+
+subroutine get_d2_new_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
+  !todo: indices/conjg should be correct for complex
+  use bitmasks
+  implicit none
+
+  integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
+  integer(bit_kind), intent(in) :: phasemask(N_int,2)
+  logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
+  double precision, intent(in) :: coefs(N_states,2)
+  double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
+  double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
+  integer, intent(in) :: h(0:2,2), p(0:4,2), sp
+
+  double precision, external :: get_phase_bi
+
+  integer :: i, j, k, tip, ma, mi, puti, putj
+  integer :: h1, h2, p1, p2, i1, i2
+  double precision :: phase
+  double precision :: hij,hji
+
+  integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0,  3,4, 0,0, 0,0,  2,4, 1,4, 0,0,  2,3, 1,3, 1,2 /), (/2,3,4/))
+  integer, parameter :: turn2(2) = (/2, 1/)
+  integer, parameter :: turn3(2,3) = reshape((/2,3,  1,3, 1,2/), (/2,3/))
+
+  integer :: bant
+  bant = 1
+
+  tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles
+
+  ma = sp !1:(alpha,alpha); 2:(b,b); 3:(a,b)
+  if(p(0,1) > p(0,2)) ma = 1 ! more alpha particles than beta particles
+  if(p(0,1) < p(0,2)) ma = 2 ! fewer alpha particles than beta particles
+  mi = mod(ma, 2) + 1
+
+  if(sp == 3) then ! if one alpha and one beta xhole 
+    !(where xholes refer to the ionizations from the generator, not the holes occupied in the ionized generator)
+    if(ma == 2) bant = 2 ! if more beta particles than alpha particles
+
+    if(tip == 3) then ! if 3 of one particle spin and 1 of the other particle spin
+      puti = p(1, mi)
+      if(bannedOrb(puti, mi)) return
+      h1 = h(1, ma)
+      h2 = h(2, ma)
+
+      !! <alpha|H|psi>
+      do i = 1, 3    ! loop over all 3 combinations of 2 particles with spin ma
+        putj = p(i, ma)
+        if(banned(putj,puti,bant)) cycle
+        i1 = turn3(1,i)
+        i2 = turn3(2,i)
+        p1 = p(i1, ma)
+        p2 = p(i2, ma)
+        
+     ! |G> = |psi_{gen,i}>
+     ! |G'> = a_{x1} a_{x2} |G>
+     ! |alpha> = a_{puti}^{\dagger} a_{putj}^{\dagger} |G'>
+     ! |alpha> = t_{x1,x2}^{puti,putj} |G>
+     ! hij = <psi_{selectors,i}|H|alpha>
+     ! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}>
+        !todo: <i|H|j>  =  (<h1,h2|p1,p2> - <h1,h2|p2,p1>) * phase
+        !    <psi|H|j> +=  dconjg(c_i) * <i|H|j>
+        !      <j|H|i>  =  (<p1,p2|h1,h2> - <p2,p1|h1,h2>) * phase
+        !    <j|H|psi> +=  <j|H|i> * c_i
+
+!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! 
+        ! take the transpose of what's written above because later use the complex conjugate 
+
+!        hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e( h1, h2, p2, p1)
+!        hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e_transp( h1, h2, p2, p1)
+        hij = mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp( p1, p2, h2, h1)
+        hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)       - mo_bi_ortho_tc_two_e( p1, p2, h2, h1)
+        if (hij == 0.d0.or.hji==0.d0) cycle
+
+        ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!        hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
+        phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
+        hij = hij * phase
+        hji = hji * phase
+
+        if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
+            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
+          enddo
+        else            ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
+          enddo
+        end if
+      end do
+    else ! if 2 alpha and 2 beta particles
+      h1 = h(1,1)
+      h2 = h(1,2)
+      !! <alpha|H|psi>
+      do j = 1,2 ! loop over all 4 combinations of one alpha and one beta particle
+        putj = p(j, 2)
+        if(bannedOrb(putj, 2)) cycle
+        p2 = p(turn2(j), 2)
+        do i = 1,2
+          puti = p(i, 1)
+          if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle
+          p1 = p(turn2(i), 1)
+    ! hij = <psi_{selectors,i}|H|alpha> 
+!          hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
+!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! 
+        ! take the transpose of what's written above because later use the complex conjugate 
+!          hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2 )
+!          hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2 )
+          hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2 )
+          hji = mo_bi_ortho_tc_two_e( p1, p2, h1, h2)
+          if (hij /= 0.d0.or.hji==0.d0) then
+            ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!            hij = dconjg(hij) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
+            phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
+            hij = hij * phase
+            hji = hji * phase
+            !DIR$ LOOP COUNT AVG(4)
+            do k=1,N_states
+              mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
+              mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
+            enddo
+          endif
+        end do
+      end do
+    end if
+
+  else ! if holes are (a,a) or (b,b)
+    if(tip == 0) then ! if particles are (a,a,a,a) or (b,b,b,b)
+      h1 = h(1, ma)
+      h2 = h(2, ma)
+      !! <alpha|H|psi>
+      do i=1,3
+        puti = p(i, ma)
+        if(bannedOrb(puti,ma)) cycle
+        do j=i+1,4
+          putj = p(j, ma)
+          if(bannedOrb(putj,ma)) cycle
+          if(banned(puti,putj,1)) cycle
+
+          i1 = turn2d(1, i, j)
+          i2 = turn2d(2, i, j)
+          p1 = p(i1, ma)
+          p2 = p(i2, ma)
+!          hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)
+!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! 
+        ! take the transpose of what's written above because later use the complex conjugate 
+          hij = mo_bi_ortho_tc_two_e_transp(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e_transp(p1, p2, h2,h1 )
+          hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p1, p2, h2,h1 )
+          if (hij == 0.d0.or.hji == 0.d0) cycle
+
+          ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!          hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
+          phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
+          hij = hij * phase
+          hji = hji * phase
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hij
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hji
+          enddo
+        end do
+      end do
+    else if(tip == 3) then ! if particles are (a,a,a,b) (ma=1,mi=2) or (a,b,b,b) (ma=2,mi=1)
+      h1 = h(1, mi)
+      h2 = h(1, ma)
+      p1 = p(1, mi)
+      !! <alpha|H|psi>
+      do i=1,3
+        puti = p(turn3(1,i), ma)
+        if(bannedOrb(puti,ma)) cycle
+        putj = p(turn3(2,i), ma)
+        if(bannedOrb(putj,ma)) cycle
+        if(banned(puti,putj,1)) cycle
+        p2 = p(i, ma)
+
+!        hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
+!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! 
+        ! take the transpose of what's written above because later use the complex conjugate 
+        hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2)
+        hji = mo_bi_ortho_tc_two_e(p1, p2,h1, h2 )
+        if (hij == 0.d0) cycle
+
+        ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!        hij = dconjg(hij) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
+        phase = get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
+        hij = hij * phase
+        hji = hji * phase
+        if (puti < putj) then
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
+          enddo
+        else
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
+            mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
+          enddo
+        endif
+      end do
+    else ! tip == 4  (a,a,b,b)
+      puti = p(1, sp)
+      putj = p(2, sp)
+      if(.not. banned(puti,putj,1)) then
+        p1 = p(1, mi)
+        p2 = p(2, mi)
+        h1 = h(1, mi)
+        h2 = h(2, mi)
+      !! <alpha|H|psi>
+!        hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2))
+!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! 
+        ! take the transpose of what's written above because later use the complex conjugate 
+        hij = (mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp(p2,p1,h1, h2))
+        hji = (mo_bi_ortho_tc_two_e(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e(p2,p1,h1, h2))
+        if (hij /= 0.d0.or.hji==0.d0) then
+          ! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
+!          hij = dconjg(hij) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
+          phase = get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
+          hij = hij * phase
+          hji = hji* phase
+          !DIR$ LOOP COUNT AVG(4)
+          do k=1,N_states
+            mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
+            mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
+          enddo
+        end if
+      end if
+    end if
+  end if
+end

plugins/local/cipsi_tc_bi_ortho/pt2.irp.f

@@ -67,6 +67,7 @@ subroutine tc_pt2
   call pt2_alloc(pt2_data_err, N_states)
   call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection
   call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+  call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2)
 
 end
 

plugins/local/cipsi_tc_bi_ortho/selection.irp.f

@@ -636,10 +636,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
     negMask(i,2) = not(mask(i,2))
   end do
 
-!  print*,'in selection '
   do i = 1, N_sel
-!    call debug_det(det(1,1,i),N_int)
-!    print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i))
     if(interesting(i) < 0) then
       stop 'prefetch interesting(i) and det(i)'
     endif
@@ -691,11 +688,23 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
 
         call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
         if(nt == 4) then
-          call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          if(transpose_two_e_int)then
+           call get_d2_new_transp(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          else
+           call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          endif
         elseif(nt == 3) then
-          call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          if(transpose_two_e_int)then
+           call get_d1_transp(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          else
+           call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+          endif
         else
+         if(transpose_two_e_int)then
+          call get_d0_transp (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+         else
           call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
+         endif
         endif
     elseif(nt == 4) then
         call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
@@ -887,79 +896,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
         do istate = 1,N_states
           delta_E = E0(istate) - Hii + E_shift
-          double precision               :: alpha_h_psi_tmp, psi_h_alpha_tmp, error
-          if(debug_tc_pt2 == 1)then !! Using the old version
-            psi_h_alpha = 0.d0
-            alpha_h_psi = 0.d0
-            do iii = 1, N_det_selectors
-              call htilde_mu_mat_bi_ortho_tot_slow(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
-              call htilde_mu_mat_bi_ortho_tot_slow(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
-              call get_excitation_degree(psi_selectors(1,1,iii), det,degree,N_int)
-              if(degree == 0)then
-               print*,'problem !!!'
-               print*,'a determinant is already in the wave function !!'
-               print*,'it corresponds to the selector number ',iii
-               call debug_det(det,N_int)
-               stop
-              endif
-!              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
-!              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
-              psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
-              alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
-            enddo
-          else if(debug_tc_pt2 == 2)then !! debugging the new version
-!            psi_h_alpha_tmp = 0.d0
-!            alpha_h_psi_tmp = 0.d0
-!            do iii = 1, N_det_selectors ! old version
-!              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
-!              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
-!              psi_h_alpha_tmp += i_h_alpha * psi_selectors_coef_tc(iii,1,1) ! left function
-!              alpha_h_psi_tmp += alpha_h_i * psi_selectors_coef_tc(iii,2,1) ! right function
-!            enddo
-            psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version
-            alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version
-            psi_h_alpha = 0.d0
-            alpha_h_psi = 0.d0
-            do iii = 1, N_det ! old version
-              call htilde_mu_mat_opt_bi_ortho_no_3e(psi_det(1,1,iii), det, N_int, i_h_alpha)
-              call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_det(1,1,iii), N_int, alpha_h_i)
-              psi_h_alpha += i_h_alpha * psi_l_coef_bi_ortho(iii,1) ! left function
-              alpha_h_psi += alpha_h_i * psi_r_coef_bi_ortho(iii,1) ! right function
-            enddo
-            if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then
-              error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi)
-              if(error.gt.1.d-2)then
-                call debug_det(det, N_int)
-                print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E
-                print*,psi_h_alpha , alpha_h_psi
-                print*,psi_h_alpha_tmp , alpha_h_psi_tmp
-                print*,'selectors '
-                do iii = 1, N_det_selectors ! old version
-                 print*,'iii',iii,psi_selectors_coef_tc(iii,1,1),psi_selectors_coef_tc(iii,2,1)
-                 call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
-                 call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
-                 print*,i_h_alpha,alpha_h_i
-                 call debug_det(psi_selectors(1,1,iii),N_int)
-                enddo
-!                print*,'psi_det '
-!                do iii = 1, N_det! old version
-!                 print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1)
-!                 call debug_det(psi_det(1,1,iii),N_int)
-!                enddo
-                stop
-              endif
-            endif
-          else
-           psi_h_alpha = mat_l(istate, p1, p2)
-           alpha_h_psi = mat_r(istate, p1, p2)
-          endif
+          psi_h_alpha = mat_l(istate, p1, p2)
+          alpha_h_psi = mat_r(istate, p1, p2)
           val = 4.d0 * psi_h_alpha * alpha_h_psi
           tmp = dsqrt(delta_E * delta_E + val)
-!          if (delta_E < 0.d0) then
-!              tmp = -tmp
-!          endif
           e_pert(istate) = 0.25 * val / delta_E
-!          e_pert(istate) = 0.5d0 * (tmp - delta_E)
           if(dsqrt(tmp).gt.1.d-4.and.dabs(psi_h_alpha).gt.1.d-4)then
            coef(istate)   = e_pert(istate) / psi_h_alpha
           else
@@ -976,15 +917,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
            if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0
           endif
 
-!         if(selection_tc     ==  1 )then
-!          if(e_pert(istate).lt.0.d0)then
-!           e_pert(istate) = 0.d0
-!          endif
-!         else if(selection_tc == -1)then
-!          if(e_pert(istate).gt.0.d0)then
-!           e_pert(istate) = 0.d0
-!          endif
-!         endif
         enddo
 
 

plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -88,6 +88,9 @@ subroutine run_stochastic_cipsi
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
     call pt2_alloc(pt2_data_err, N_states)
+    if(transpose_two_e_int)then
+     provide mo_bi_ortho_tc_two_e_transp
+    endif
     call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
 !    stop
 

plugins/local/fci_tc_bi/pt2_tc.irp.f

@@ -13,6 +13,8 @@ program tc_pt2_prog
 
   pruning = -1.d0
   touch pruning
+  read_wf = .True.
+  touch read_wf 
 
 !  pt2_relative_error = 0.01d0
 !  touch pt2_relative_error

plugins/local/old_delta_tc_qmc/NEED

@@ -0,0 +1 @@
+

plugins/local/old_delta_tc_qmc/README.rst

@@ -0,0 +1,4 @@
+================
+old_delta_tc_qmc
+================
+

plugins/local/old_delta_tc_qmc/dressing_vectors_lr.irp.f

@@ -27,7 +27,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
 
   i = 1
   j = 1
-  call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
+  call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
   call hmat_bi_ortho         (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
 
   delta = 0.d0
@@ -39,7 +39,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
     do j = 1, ndet
 
       ! < I |Htilde | J >
-      call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
+      call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
       ! < I |H | J >
       call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
 
@@ -78,7 +78,7 @@ subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta)
 
   i = 1
   j = 1
-  call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
+  call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
 
   delta = 0.d0
  !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8)   &
@@ -88,7 +88,7 @@ subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta)
     do j = 1, ndet
 
       ! < I |Htilde | J >
-      call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
+      call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
 
       delta(i) = delta(i) + psicoef(j) * htc_tot
     enddo

plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f

@@ -1,4 +1,4 @@
-program slater_tc
+program old_delta_tc_qmc
   implicit none
   BEGIN_DOC
 ! TODO : Put the documentation of the program here

plugins/local/slater_tc/h_mat_triple.irp.f

@@ -1,196 +1,3 @@
-subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
- use bitmasks
- BEGIN_DOC
-! returns the array, for each spin, of holes/particles between key_i and key_j 
-!
-! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
- END_DOC
-  include 'utils/constants.include.F'
- implicit none
- integer, intent(in)            :: Nint
- integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
- integer, intent(out)           :: holes_array(100,2),particles_array(100,2),degree_array(2)
- double precision, intent(out)  :: phase
- integer :: ispin,k,i,pos
- integer(bit_kind) :: key_hole, key_particle
- integer(bit_kind) :: xorvec(N_int_max,2)
- holes_array = -1
- particles_array = -1
- degree_array = 0
-  do i = 1, N_int
-   xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
-   xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
-   degree_array(1) += popcnt(xorvec(i,1))
-   degree_array(2) += popcnt(xorvec(i,2))
-  enddo
-  degree_array(1) = shiftr(degree_array(1),1)
-  degree_array(2) = shiftr(degree_array(2),1)
-  
- do ispin = 1, 2
-  k = 1
-  !!! GETTING THE HOLES 
-  do i = 1, N_int
-   key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
-   do while(key_hole .ne.0_bit_kind)
-    pos = trailz(key_hole)
-    holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
-    key_hole = ibclr(key_hole,pos)
-    k += 1
-    if(k .gt.100)then
-     print*,'WARNING in get_excitation_general'
-     print*,'More than a 100-th excitation for spin ',ispin
-     print*,'stoping ...'
-     stop
-    endif
-   enddo 
-  enddo
- enddo
- do ispin = 1, 2
-  k = 1
-  !!! GETTING THE PARTICLES
-  do i = 1, N_int
-   key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
-   do while(key_particle .ne.0_bit_kind)
-    pos = trailz(key_particle)
-    particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
-    key_particle = ibclr(key_particle,pos)
-    k += 1
-    if(k .gt.100)then
-     print*,'WARNING in get_excitation_general '
-     print*,'More than a 100-th excitation for spin ',ispin
-     print*,'stoping ...'
-     stop
-    endif
-   enddo
-  enddo 
- enddo
- integer :: h,p, i_ok
- integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
- integer                        :: exc(0:2,2,2)
- double precision :: phase_tmp
- allocate(det_i(Nint,2),det_ip(N_int,2))
- det_i = key_i
- phase = 1.d0
- do ispin = 1, 2
-  do i = 1, degree_array(ispin)
-   h = holes_array(i,ispin)
-   p = particles_array(i,ispin)
-   det_ip = det_i
-   call do_single_excitation(det_ip,h,p,ispin,i_ok)
-   if(i_ok == -1)then
-     print*,'excitation was not possible '
-     stop
-   endif
-   call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
-   phase *= phase_tmp
-   det_i = det_ip
-  enddo
- enddo
-
-end
-
-subroutine get_holes_general(key_i, key_j,Nint, holes_array)
- use bitmasks
- BEGIN_DOC
-! returns the array, per spin, of holes between key_i and key_j 
-!
-! with the following convention: a_{hole}|key_i> --> |key_j>
- END_DOC
- implicit none
- integer, intent(in)            :: Nint
- integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
- integer, intent(out)           :: holes_array(100,2)
- integer(bit_kind) :: key_hole
- integer :: ispin,k,i,pos
- holes_array = -1
- do ispin = 1, 2
-  k = 1
-  do i = 1, N_int
-   key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
-   do while(key_hole .ne.0_bit_kind)
-    pos = trailz(key_hole)
-    holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
-    key_hole = ibclr(key_hole,pos)
-    k += 1
-    if(k .gt.100)then
-     print*,'WARNING in get_holes_general'
-     print*,'More than a 100-th excitation for spin ',ispin
-     print*,'stoping ...'
-     stop
-    endif
-   enddo 
-  enddo
- enddo
-end
-
-subroutine get_particles_general(key_i, key_j,Nint,particles_array)
- use bitmasks
- BEGIN_DOC
-! returns the array, per spin, of particles between key_i and key_j 
-!
-! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
- END_DOC
- implicit none
- integer, intent(in)            :: Nint
- integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
- integer, intent(out)           :: particles_array(100,2)
- integer(bit_kind) :: key_particle
- integer :: ispin,k,i,pos
- particles_array = -1
- do ispin = 1, 2
-  k = 1
-  do i = 1, N_int
-   key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
-   do while(key_particle .ne.0_bit_kind)
-    pos = trailz(key_particle)
-    particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
-    key_particle = ibclr(key_particle,pos)
-    k += 1
-    if(k .gt.100)then
-     print*,'WARNING in get_holes_general'
-     print*,'More than a 100-th excitation for spin ',ispin
-     print*,'Those are the two determinants'
-     call debug_det(key_i, N_int)
-     call debug_det(key_j, N_int)
-     print*,'stoping ...'
-     stop
-    endif
-   enddo 
-  enddo
- enddo
-end
-
-subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
- implicit none
- integer, intent(in)            :: degree(2), Nint
- integer(bit_kind), intent(in)  :: key_i(Nint,2)
- integer, intent(in)            :: holes_array(100,2),particles_array(100,2)
- double precision, intent(out)  :: phase
- integer :: i,ispin,h,p, i_ok
- integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
- integer                        :: exc(0:2,2,2)
- double precision :: phase_tmp
- allocate(det_i(Nint,2),det_ip(N_int,2))
- det_i = key_i
- phase = 1.d0
- do ispin = 1, 2
-  do i = 1, degree(ispin)
-   h = holes_array(i,ispin)
-   p = particles_array(i,ispin)
-   det_ip = det_i
-   call do_single_excitation(det_ip,h,p,ispin,i_ok)
-   if(i_ok == -1)then
-     print*,'excitation was not possible '
-     stop
-   endif
-   call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
-   phase *= phase_tmp
-   det_i = det_ip
-  enddo
- enddo
-
-end
-
 subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
   BEGIN_DOC
   ! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS 

plugins/local/slater_tc/slater_tc_opt.irp.f

@@ -181,3 +181,48 @@ end
 
 ! ---
 
+subroutine htilde_mu_mat_opt_bi_ortho_no_3e_both(key_j, key_i, Nint, hji,hij)
+
+  BEGIN_DOC
+  !
+  ! <key_j |H_tilde | key_i> 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 WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS
+  !! 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) :: hji,hij
+  integer                       :: degree
+
+  hji = 0.d0
+  hij = 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_no_3e(Nint, key_i,hji)
+    hij = hji
+  else if (degree == 1) then
+    call single_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint,key_j, key_i , hji,hij)
+  else if(degree == 2) then
+    call double_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
+  endif
+
+  if(degree==0) then
+    hji += nuclear_repulsion
+    hij += nuclear_repulsion
+  endif
+
+end
+
+! ---
+

plugins/local/slater_tc/slater_tc_opt_diag.irp.f

@@ -19,13 +19,13 @@
   PROVIDE HF_bitmask
   PROVIDE mo_l_coef mo_r_coef
 
-  call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot)
+  call diag_htc_bi_orth_2e_brute(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_slow(N_int, HF_bitmask, hthree)
+    call diag_htc_bi_orth_3e_brute(N_int, HF_bitmask, hthree)
     ref_tc_energy_3e = hthree
   else
     ref_tc_energy_3e = 0.d0
@@ -524,3 +524,310 @@ end
 
 ! ---
 
+subroutine diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
+
+  BEGIN_DOC
+  !
+  ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS 
+  !
+  END_DOC
+
+  use bitmasks
+
+  implicit none
+  integer,           intent(in)  :: Nint
+  integer(bit_kind), intent(in)  :: 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
+  double precision               :: get_mo_two_e_integral_tc_int
+  integer(bit_kind)              :: key_i_core(Nint,2)
+
+  PROVIDE mo_bi_ortho_tc_two_e
+
+  hmono = 0.d0
+  htwoe = 0.d0
+  htot  = 0.d0
+
+  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+
+  do ispin = 1, 2
+    do i = 1, Ne(ispin)
+      ii = occ(i,ispin)
+      hmono += mo_bi_ortho_tc_one_e(ii,ii)
+    enddo
+  enddo
+
+  ! alpha/beta two-body
+  ispin = 1
+  jspin = 2
+  do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
+    ii = occ(i,ispin) 
+    do j = 1, Ne(jspin) ! electron 2 
+      jj = occ(j,jspin) 
+      htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii) 
+    enddo
+  enddo
+ 
+  ! alpha/alpha two-body
+  do i = 1, Ne(ispin)
+    ii = occ(i,ispin) 
+    do j = i+1, Ne(ispin)
+      jj = occ(j,ispin) 
+      htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
+    enddo
+  enddo
+ 
+  ! beta/beta two-body
+  do i = 1, Ne(jspin)
+    ii = occ(i,jspin) 
+    do j = i+1, Ne(jspin)
+      jj = occ(j,jspin) 
+      htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
+    enddo
+  enddo
+
+  htot = hmono + htwoe 
+
+end
+
+! ---                                                                                           
+
+subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
+
+  BEGIN_DOC
+  !  diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
+  END_DOC
+
+  use bitmasks
+
+  implicit none
+  integer,           intent(in) :: Nint
+  integer(bit_kind), intent(in) :: key_i(Nint,2)
+  double precision, intent(out) :: hthree
+  integer                       :: occ(Nint*bit_kind_size,2)
+  integer                       :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
+  integer(bit_kind)             :: key_i_core(Nint,2)
+  double precision              :: direct_int, exchange_int, ref
+  double precision, external    :: sym_3_e_int_from_6_idx_tensor
+  double precision, external    :: three_e_diag_parrallel_spin
+
+  PROVIDE mo_l_coef mo_r_coef
+
+  if(core_tc_op) then
+    do i = 1, Nint
+      key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1))
+      key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2))
+    enddo
+    call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
+  else
+    call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+  endif
+
+  hthree = 0.d0
+
+  if((Ne(1)+Ne(2)) .ge. 3) then
+
+    ! alpha/alpha/beta three-body
+    do i = 1, Ne(1)
+      ii = occ(i,1)
+      do j = i+1, Ne(1)
+        jj = occ(j,1)
+        do m = 1, Ne(2)
+          mm = occ(m,2)
+          !direct_int   = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor
+          !exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor
+          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
+        enddo
+      enddo
+    enddo
+
+    ! beta/beta/alpha three-body
+    do i = 1, Ne(2)
+      ii = occ(i,2)
+      do j = i+1, Ne(2)
+        jj = occ(j,2)
+        do m = 1, Ne(1)
+          mm = occ(m,1)
+          !direct_int   = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor
+          !exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor
+          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
+        enddo
+      enddo
+    enddo
+
+    ! alpha/alpha/alpha three-body
+    do i = 1, Ne(1)
+      ii = occ(i,1) ! 1
+      do j = i+1, Ne(1)
+        jj = occ(j,1) ! 2
+        do m = j+1, Ne(1)
+          mm = occ(m,1) ! 3
+          !hthree += 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
+        enddo
+      enddo
+    enddo
+
+    ! beta/beta/beta three-body
+    do i = 1, Ne(2)
+      ii = occ(i,2) ! 1
+      do j = i+1, Ne(2)
+        jj = occ(j,2) ! 2
+        do m = j+1, Ne(2)
+          mm = occ(m,2) ! 3
+          !hthree += 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
+        enddo
+      enddo
+    enddo
+
+  endif
+
+end
+
+
+
+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 <mji|-L|mji> 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 <mjk|-L|mji> 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) = <mlk|-L|mji> ::: 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 
+

plugins/local/slater_tc/slater_tc_opt_double.irp.f

@@ -505,3 +505,63 @@ subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
 
 end
 
+subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
+
+  BEGIN_DOC
+  ! <key_j |H_tilde | key_i> and <key_i |H_tilde | key_j> 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) :: hji,hij
+  double precision :: hmono, htwoe_ji, htwoe_ij
+  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_ji  = 0.d0
+  htwoe_ij  = 0.d0
+  hji   = 0.d0
+  hij   = 0.d0
+
+  if(degree.ne.2)then
+   return
+  endif
+  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
+    htwoe_ji  = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+    htwoe_ij  = mo_bi_ortho_tc_two_e_transp(p2,p1,h2,h1)
+  else
+   ! same spin two-body
+   ! direct terms
+   htwoe_ji  = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
+   htwoe_ij  = mo_bi_ortho_tc_two_e_transp(p2,p1,h2,h1)
+   ! exchange terms
+   htwoe_ji -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
+   htwoe_ij -= mo_bi_ortho_tc_two_e_transp(p1,p2,h2,h1)
+  endif
+  htwoe_ji  *= phase
+  hji    =  htwoe_ji
+  htwoe_ij  *= phase
+  hij    =  htwoe_ij
+
+end

plugins/local/slater_tc/slater_tc_opt_single.irp.f

@@ -618,3 +618,145 @@ subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, sp
 
 end
 
+
+subroutine  single_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
+
+  BEGIN_DOC
+  ! <key_j |H_tilde | key_i> and <key_i |H_tilde | key_j> 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) :: hji,hij
+
+  double precision              :: hmono, htwoe
+  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
+  hji   = 0.d0
+  hij   = 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_no_3e_both(Nint, key_i, key_j, h1, p1, s1, phase, hji,hij)
+
+end
+
+! ---
+
+subroutine get_single_excitation_from_fock_tc_no_3e_both(Nint, key_i, key_j, h, p, spin, phase, hji,hij)
+
+  use bitmasks
+
+  implicit none
+  integer,           intent(in) :: Nint
+  integer,           intent(in) :: h, p, spin
+  double precision,  intent(in) :: phase
+  integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out) :: hji,hij
+  double precision :: hmono_ji,htwoe_ji
+  double precision :: hmono_ij,htwoe_ij
+
+  integer(bit_kind)             :: differences(Nint,2)
+  integer(bit_kind)             :: hole(Nint,2)
+  integer(bit_kind)             :: partcl(Nint,2)
+  integer                       :: occ_hole(Nint*bit_kind_size,2)
+  integer                       :: occ_partcl(Nint*bit_kind_size,2)
+  integer                       :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
+  integer                       :: i0,i
+  double precision              :: buffer_c_ji(mo_num), buffer_x_ji(mo_num)
+  double precision              :: buffer_c_ij(mo_num), buffer_x_ij(mo_num)
+
+  do i = 1, mo_num
+    buffer_c_ji(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
+    buffer_x_ji(i) = tc_2e_3idx_exchange_integrals(i,p,h)
+    buffer_c_ij(i) = tc_2e_3idx_coulomb_integrals_transp(i,p,h)
+    buffer_x_ij(i) = tc_2e_3idx_exchange_integrals_transp(i,p,h)
+  enddo
+
+  do i = 1, Nint
+   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, Nint)
+  call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
+  hmono_ji = mo_bi_ortho_tc_one_e(p,h)
+  htwoe_ji = fock_op_2_e_tc_closed_shell(p,h)
+  hmono_ij = mo_bi_ortho_tc_one_e(h,p)
+  htwoe_ij = fock_op_2_e_tc_closed_shell(h,p)
+
+  ! holes :: direct terms
+  do i0 = 1, n_occ_ab_hole(1)
+    i = occ_hole(i0,1)
+    htwoe_ji -= buffer_c_ji(i)
+    htwoe_ij -= buffer_c_ij(i)
+  enddo
+  do i0 = 1, n_occ_ab_hole(2)
+    i = occ_hole(i0,2)
+    htwoe_ji -= buffer_c_ji(i)
+    htwoe_ij -= buffer_c_ij(i)
+  enddo
+ 
+  ! holes :: exchange terms
+  do i0 = 1, n_occ_ab_hole(spin)
+    i = occ_hole(i0,spin)
+    htwoe_ji += buffer_x_ji(i)
+    htwoe_ij += buffer_x_ij(i)
+  enddo
+ 
+  ! particles :: direct terms
+  do i0 = 1, n_occ_ab_partcl(1)
+    i = occ_partcl(i0,1)
+    htwoe_ji += buffer_c_ji(i)
+    htwoe_ij += buffer_c_ij(i)
+  enddo
+  do i0 = 1, n_occ_ab_partcl(2)
+    i = occ_partcl(i0,2)
+    htwoe_ji += buffer_c_ji(i)
+    htwoe_ij += buffer_c_ij(i)
+  enddo
+ 
+  ! particles :: exchange terms
+  do i0 = 1, n_occ_ab_partcl(spin)
+    i = occ_partcl(i0,spin)
+    htwoe_ji -= buffer_x_ji(i)
+    htwoe_ij -= buffer_x_ij(i)
+  enddo
+  htwoe_ji = htwoe_ji * phase
+  hmono_ji = hmono_ji * phase
+  hji  = htwoe_ji + hmono_ji 
+
+  htwoe_ij = htwoe_ij * phase
+  hmono_ij = hmono_ij * phase
+  hij  = htwoe_ij + hmono_ij 
+
+end
+

plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f

@@ -1,140 +0,0 @@
-
-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 <mji|-L|mji> 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 <mjk|-L|mji> 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) = <mlk|-L|mji> ::: 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 
-

plugins/local/slater_tc_no_opt/.gitignore

@@ -0,0 +1,59 @@
+IRPF90_temp/
+IRPF90_man/
+build.ninja
+irpf90.make
+ezfio_interface.irp.f
+irpf90_entities
+tags
+Makefile
+ao_basis
+ao_one_e_ints
+ao_two_e_erf_ints
+ao_two_e_ints
+aux_quantities
+becke_numerical_grid
+bitmask
+cis
+cisd
+cipsi
+davidson
+davidson_dressed
+davidson_undressed
+density_for_dft
+determinants
+dft_keywords
+dft_utils_in_r
+dft_utils_one_e
+dft_utils_two_body
+dressing
+dummy
+electrons
+ezfio_files
+fci
+generators_cas
+generators_full
+hartree_fock
+iterations
+kohn_sham
+kohn_sham_rs
+mo_basis
+mo_guess
+mo_one_e_ints
+mo_two_e_erf_ints
+mo_two_e_ints
+mpi
+mrpt_utils
+nuclei
+perturbation
+pseudo
+psiref_cas
+psiref_utils
+scf_utils
+selectors_cassd
+selectors_full
+selectors_utils
+single_ref_method
+slave
+tools
+utils
+zmq

plugins/local/slater_tc_no_opt/NEED

@@ -0,0 +1,8 @@
+determinants
+normal_order_old
+bi_ort_ints
+bi_ortho_mos
+tc_keywords
+non_hermit_dav
+dav_general_mat
+tc_scf

plugins/local/slater_tc_no_opt/README.rst

@@ -0,0 +1,4 @@
+================
+slater_tc_no_opt
+================
+

plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f

@@ -1,7 +1,7 @@
 
 ! ---
 
-subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
+subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
 
   BEGIN_DOC
   !  diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS

plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f

@@ -0,0 +1,7 @@
+program slater_tc_no_opt
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  print *, 'Hello world'
+end

plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f

@@ -61,7 +61,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
   if(degree.gt.2) return
 
   if(degree == 0) then
-    call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
+    call diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
   else if (degree == 1) then
     call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
   else if(degree == 2) then
@@ -76,7 +76,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
     else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then
       call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
     else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then
-      call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
+      call diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
     endif
   endif
 
@@ -95,75 +95,6 @@ end
 
 ! ---
 
-subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
-
-  BEGIN_DOC
-  !
-  ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS 
-  !
-  END_DOC
-
-  use bitmasks
-
-  implicit none
-  integer,           intent(in)  :: Nint
-  integer(bit_kind), intent(in)  :: 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
-  double precision               :: get_mo_two_e_integral_tc_int
-  integer(bit_kind)              :: key_i_core(Nint,2)
-
-  PROVIDE mo_bi_ortho_tc_two_e
-
-  hmono = 0.d0
-  htwoe = 0.d0
-  htot  = 0.d0
-
-  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
-
-  do ispin = 1, 2 
-    do i = 1, Ne(ispin)
-      ii = occ(i,ispin) 
-      hmono += mo_bi_ortho_tc_one_e(ii,ii)
-    enddo
-  enddo
-
-  ! alpha/beta two-body
-  ispin = 1
-  jspin = 2 
-  do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
-    ii = occ(i,ispin) 
-    do j = 1, Ne(jspin) ! electron 2 
-      jj = occ(j,jspin) 
-      htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii) 
-    enddo
-  enddo
- 
-  ! alpha/alpha two-body
-  do i = 1, Ne(ispin)
-    ii = occ(i,ispin) 
-    do j = i+1, Ne(ispin)
-      jj = occ(j,ispin) 
-      htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
-    enddo
-  enddo
- 
-  ! beta/beta two-body
-  do i = 1, Ne(jspin)
-    ii = occ(i,jspin) 
-    do j = i+1, Ne(jspin)
-      jj = occ(j,jspin) 
-      htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
-    enddo
-  enddo
-
-  htot = hmono + htwoe 
-
-end
-
-! ---
-
 subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
 
   BEGIN_DOC

plugins/local/slater_tc_no_opt/test_tc_bi_ortho.irp.f

@@ -88,7 +88,7 @@ subroutine test_slater_tc_opt
  i_count = 0.d0
  do i = 1, N_det
   do j = 1,N_det
-   call htilde_mu_mat_bi_ortho_slow(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, 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(htot).gt.1.d-15)then
      i_count += 1.D0
@@ -124,7 +124,7 @@ subroutine timing_tot
   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_slow(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, hmono, htwoe, hthree, htot)
   enddo
  enddo
  call wall_time(wall1)
@@ -171,7 +171,7 @@ subroutine timing_diag
  do i = 1, N_det
   do j = i,i 
    i_count += 1.d0
-   call htilde_mu_mat_bi_ortho_slow(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, hmono, htwoe, hthree, htot)
   enddo
  enddo
  call wall_time(wall1)
@@ -208,7 +208,7 @@ subroutine timing_single
    if(degree.ne.1)cycle
    i_count += 1.d0
    call wall_time(wall0)
-   call htilde_mu_mat_bi_ortho_slow(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, hmono, htwoe, hthree, htot)
    call wall_time(wall1)
    accu += wall1 - wall0
   enddo
@@ -250,7 +250,7 @@ subroutine timing_double
    if(degree.ne.2)cycle
    i_count += 1.d0
    call wall_time(wall0)
-   call htilde_mu_mat_bi_ortho_slow(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, hmono, htwoe, hthree, htot)
    call wall_time(wall1)
    accu += wall1 - wall0
   enddo

plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f

@@ -2,7 +2,7 @@
  BEGIN_PROVIDER [ double precision, e_tilde_00]
  implicit none
  double precision :: hmono,htwoe,hthree,htot
- call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
+ call htilde_mu_mat_opt_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
  e_tilde_00 = htot
  END_PROVIDER 
 
@@ -18,16 +18,15 @@
  do i = 1, N_det
   call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
   if(degree == 1 .or. degree == 2)then
-   call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
-   call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
    delta_e = e_tilde_00 - e_i0
    coef_pt1 = htilde_ij / delta_e
-   call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+   call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
    e_pt2_tc_bi_orth += coef_pt1 * htilde_ij
    if(degree == 1)then
     e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij
    else 
-!    print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij
     e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij
    endif
   endif
@@ -37,7 +36,7 @@
  BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00]
  implicit none
  double precision :: hmono,htwoe,hthree,htilde_ij
- call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
+ call htilde_mu_mat_opt_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
  e_tilde_bi_orth_00 += nuclear_repulsion
  END_PROVIDER 
 
@@ -57,7 +56,7 @@
  e_corr_double_bi_orth = 0.d0
  do i = 1, N_det
   call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
-  call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+  call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
   if(degree == 1)then
    e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
    e_corr_single_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1))
@@ -80,7 +79,7 @@
  do i = 1, N_det
   accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1)
   do j = 1, N_det
-   call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
    e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1)
   enddo
  enddo
@@ -99,8 +98,8 @@ BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)]
   if(degree==0)then
    coef_pt1_bi_ortho(i) = 1.d0
   else
-   call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
-   call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+   call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
    delta_e = e_tilde_00 - e_i0
    coef_pt1 = htilde_ij / delta_e
    coef_pt1_bi_ortho(i)= coef_pt1

plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f

@@ -1,129 +0,0 @@
-program pt2_tc_cisd
-
-  BEGIN_DOC
-  !
-  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together 
-  !        with the energy. Saves the left-right wave functions at the end. 
-  !
-  END_DOC
-
-  implicit none
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  print*, ' nb of states = ', N_states
-  print*, ' nb of det    = ', N_det
-  call routine_diag()
-
-  call routine
-end
-
-subroutine routine
- implicit none
- integer :: i,h1,p1,h2,p2,s1,s2,degree
- double precision :: h0i,hi0,e00,ei,delta_e
- double precision :: norm,e_corr,coef,e_corr_pos,e_corr_neg,e_corr_abs
-
- integer                        :: exc(0:2,2,2)
- double precision               :: phase
- double precision :: eh1,ep1,eh2,ep2
-
- norm = 0.d0 
- e_corr = 0.d0
- e_corr_abs = 0.d0
- e_corr_pos = 0.d0
- e_corr_neg = 0.d0
- call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,1), N_int, e00) 
- do i = 2, N_det
-  call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0) 
-  call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i) 
-  call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, ei) 
-  call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i),degree,N_int)
-  call get_excitation(psi_det(1,1,1), psi_det(1,1,i),exc,degree,phase,N_int)
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-  eh1 = Fock_matrix_tc_diag_mo_tot(h1)
-  ep1 = Fock_matrix_tc_diag_mo_tot(p1)
-  delta_e = eh1 - ep1
-  if (degree==2)then
-   eh2 = Fock_matrix_tc_diag_mo_tot(h2)
-   ep2 = Fock_matrix_tc_diag_mo_tot(p2)
-   delta_e +=  eh2 - ep2
-  endif
-!  delta_e = e00 - ei
-  coef = hi0/delta_e
-  norm += coef*coef
-  e_corr = coef* h0i
-  if(e_corr.lt.0.d0)then
-   e_corr_neg += e_corr
-  elseif(e_corr.gt.0.d0)then
-   e_corr_pos += e_corr
-  endif
-  e_corr_abs += dabs(e_corr)
- enddo
- print*,'e_corr_abs = ',e_corr_abs
- print*,'e_corr_pos = ',e_corr_pos
- print*,'e_corr_neg = ',e_corr_neg
- print*,'norm       = ',dsqrt(norm)
-
-end
-
-subroutine routine_diag()
-
-  implicit none
-  integer          :: i, j, k
-  double precision :: dE
-
-  ! provide eigval_right_tc_bi_orth
-  ! provide overlap_bi_ortho
-  ! provide htilde_matrix_elmt_bi_ortho
-
-  if(N_states .eq. 1) then
-
-    print*,'eigval_right_tc_bi_orth   = ',eigval_right_tc_bi_orth(1)
-    print*,'e_tc_left_right           = ',e_tc_left_right
-    print*,'e_tilde_bi_orth_00        = ',e_tilde_bi_orth_00
-    print*,'e_pt2_tc_bi_orth          = ',e_pt2_tc_bi_orth
-    print*,'e_pt2_tc_bi_orth_single   = ',e_pt2_tc_bi_orth_single
-    print*,'e_pt2_tc_bi_orth_double   = ',e_pt2_tc_bi_orth_double
-    print*,'***'                      
-    print*,'e_corr_bi_orth            = ',e_corr_bi_orth
-    print*,'e_corr_bi_orth_proj       = ',e_corr_bi_orth_proj
-    print*,'e_corr_bi_orth_proj_abs   = ',e_corr_bi_orth_proj_abs
-    print*,'e_corr_single_bi_orth     = ',e_corr_single_bi_orth
-    print*,'e_corr_double_bi_orth     = ',e_corr_double_bi_orth
-    print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs
-    print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs
-    print*,'Left/right eigenvectors'
-    do i = 1,N_det
-      write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
-    enddo
-
-  else
-
-    print*,'eigval_right_tc_bi_orth : '
-    do i = 1, N_states
-      print*, i, eigval_right_tc_bi_orth(i)
-    enddo
-
-    print*,''
-    print*,'******************************************************'
-    print*,'TC Excitation energies (au)                     (eV)'
-    do i = 2, N_states
-      dE = eigval_right_tc_bi_orth(i) - eigval_right_tc_bi_orth(1)
-      print*, i, dE, dE/0.0367502d0
-    enddo
-    print*,''
-
-  endif
-
-end
-
-
-

plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f

@@ -1,36 +0,0 @@
-
-! ---
-
-program tc_cisd_sc2
-
-  BEGIN_DOC
-  ! TODO : Put the documentation of the program here
-  END_DOC
-
-  implicit none
-
-  print *, 'Hello world'
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  call test
-
-end
-
-! ---
-
-subroutine test()
- implicit none
-! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
-! allocate(dressing_dets(N_det),e_corr_dets(N_det))
-! e_corr_dets = 0.d0
-! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
-  provide eigval_tc_cisd_sc2_bi_ortho
-end

plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f

@@ -1,145 +0,0 @@
- BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
-&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
-&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)]
- implicit none
- integer :: it,n_real,degree,i,istate
- double precision :: e_before, e_current,thr, hmono,htwoe,hthree,accu
- double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:)
- double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
- allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
- allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det))
- allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),eigval_tmp(N_states))
- dressing_dets = 0.d0
- do i = 1, N_det
-  call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
-  call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
-  if(degree == 1 .or. degree == 2)then
-   call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
-  endif
- enddo
- reigvec_tc_bi_orth_tmp = 0.d0
- do i = 1, N_det 
-  reigvec_tc_bi_orth_tmp(i,1) = psi_r_coef_bi_ortho(i,1) 
- enddo
- vec_tmp = 0.d0
- do istate = 1, N_states
-  vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
- enddo
- do istate = N_states+1, n_states_diag
-  vec_tmp(istate,istate) = 1.d0
- enddo
- print*,'Diagonalizing the TC CISD '
- call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
- do i = 1, N_det 
-  e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
- enddo
- E_before = eigval_tmp(1)
- print*,'Starting from ',E_before
-
- e_current = 10.d0
- thr = 1.d-5
- it = 0
- dressing_dets = 0.d0
-  double precision, allocatable :: H_jj(:),vec_tmp(:,:),eigval_tmp(:)
-  external                         htc_bi_ortho_calc_tdav_slow
-  external                         htcdag_bi_ortho_calc_tdav_slow
-  logical                       :: converged
- do while (dabs(E_before-E_current).gt.thr)
-  it += 1
-  E_before = E_current
-!  h_sc2 = htilde_matrix_elmt_bi_ortho
-  call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
-  do i = 1, N_det
-!   print*,'dressing_dets(i) = ',dressing_dets(i)
-   h_sc2(i,i) += dressing_dets(i)
-  enddo
-  print*,'********************'
-  print*,'iteration       ',it
-!  call non_hrmt_real_diag(N_det,h_sc2,& 
-!       leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& 
-!       n_real,eigval_right_tmp)
-!  print*,'eigval_right_tmp(1)',eigval_right_tmp(1)
-  vec_tmp = 0.d0
-  do istate = 1, N_states
-   vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
-  enddo
-  do istate = N_states+1, n_states_diag
-   vec_tmp(istate,istate) = 1.d0
-  enddo
-  call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
-  print*,'outside Davidson'
-  print*,'eigval_tmp(1) = ',eigval_tmp(1)
-  do i = 1, N_det 
-   reigvec_tc_bi_orth_tmp(i,1) = vec_tmp(i,1)
-   e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
-  enddo
-!  E_current = eigval_right_tmp(1)
-  E_current = eigval_tmp(1)
-  print*,'it, E(SC)^2 = ',it,E_current
- enddo
- eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states)
- reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
- leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
- 
-END_PROVIDER 
-
-subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets)
- implicit none
-  use bitmasks
- integer, intent(in) :: ndet
- integer(bit_kind), intent(in)  :: dets(N_int,2,ndet)
- double precision, intent(in)   :: e_corr_dets(ndet)
- double precision, intent(out) :: dressing_dets(ndet)
- integer, allocatable  :: degree(:),hole(:,:),part(:,:),spin(:,:)
- integer(bit_kind), allocatable :: hole_part(:,:,:)
- integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2
- integer(bit_kind) :: xorvec(2,N_int)
-
- double precision :: phase
- dressing_dets = 0.d0
- allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet))
- do i = 2, ndet
-  call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int)
-  do j = 1, N_int
-   hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i))
-   hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i))
-  enddo
-  if(degree(i) == 1)then
-   call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
-  else if(degree(i) == 2)then
-   call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
-  endif
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-  hole(1,i) = h1
-  hole(2,i) = h2
-  part(1,i) = p1
-  part(2,i) = p2
-  spin(1,i) = s1
-  spin(2,i) = s2
- enddo
- 
- integer :: same
- if(elec_alpha_num+elec_beta_num<3)return
- do i = 2, ndet
-  do j = i+1, ndet
-   same = 0
-   if(degree(i) == degree(j) .and. degree(i)==1)cycle
-   do k = 1, N_int
-    xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j))
-    xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j))
-    same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2)) 
-   enddo
-!   print*,'i,j',i,j
-!   call debug_det(dets(1,1,i),N_int) 
-!   call debug_det(hole_part(1,1,i),N_int) 
-!   call debug_det(dets(1,1,j),N_int) 
-!   call debug_det(hole_part(1,1,j),N_int) 
-!   print*,'same = ',same
-   if(same.eq.0)then
-    dressing_dets(i) += e_corr_dets(j)  
-    dressing_dets(j) += e_corr_dets(i)  
-   endif
-  enddo
- enddo
- 
-end

plugins/local/tc_bi_ortho/tc_utils.irp.f

@@ -25,7 +25,7 @@ subroutine write_tc_energy()
     E_2e_tmp(i) = 0.d0
     E_3e_tmp(i) = 0.d0
     do j = 1, N_det
-      call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+      call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
       E_TC_tmp(i) = E_TC_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htot
       E_1e_tmp(i) = E_1e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * hmono
       E_2e_tmp(i) = E_2e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htwoe
@@ -70,7 +70,7 @@ subroutine write_tc_energy()
     E_3e = 0.d0
     do i = 1, N_det
       do j = 1, N_det
-        call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+        call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
         E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
         E_1e = E_1e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hmono
         E_2e = E_2e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htwoe
@@ -109,8 +109,8 @@ subroutine write_tc_var()
 
     SIGMA_TC = 0.d0
     do j = 2, N_det
-      call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,1), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot_1j)
-      call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot_j1)
+      call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot_1j)
+      call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot_j1)
       SIGMA_TC = SIGMA_TC + htot_1j * htot_j1
     enddo
 
@@ -132,7 +132,7 @@ subroutine write_tc_gs_var_HF()
 
   SIGMA_TC = 0.d0
   do j = 2, N_det
-    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
     SIGMA_TC = SIGMA_TC + htot * htot
   enddo
 

plugins/local/tc_bi_ortho/test_natorb.irp.f

@@ -1,64 +0,0 @@
-
-! ---
-
-program test_natorb
-
-  BEGIN_DOC
-  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
-  END_DOC
-
-  implicit none
-
-  print *, 'Hello world'
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  call routine()
-  ! call test()
-
-end
-
-! ---
-
-subroutine routine()
-
- implicit none
- double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:)
- allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num))
- double precision, allocatable :: eigval_ref(:),leigvec_ref(:,:),reigvec_ref(:,:)
- allocate(eigval_ref(mo_num),leigvec_ref(mo_num,mo_num),reigvec_ref(mo_num,mo_num))
- 
- double precision :: thr_deg
- integer :: i,n_real,j
- print*,'fock_matrix'
- do i = 1, mo_num
-  fock_diag(i) = Fock_matrix_mo(i,i)
-  print*,i,fock_diag(i)
- enddo
- thr_deg = 1.d-6
- mat_ref = -one_e_dm_mo
- print*,'diagonalization by block'
- call diag_mat_per_fock_degen(fock_diag,mat_ref,mo_num,thr_deg,leigvec,reigvec,eigval)
- call non_hrmt_bieig( mo_num, mat_ref&
-                     , leigvec_ref, reigvec_ref& 
-                     , n_real, eigval_ref)
- print*,'TEST ***********************************'
- double precision :: accu_l, accu_r
- do i = 1, mo_num
-  accu_l = 0.d0
-  accu_r = 0.d0
-  do j = 1, mo_num
-   accu_r += reigvec_ref(j,i) * reigvec(j,i)
-   accu_l += leigvec_ref(j,i) * leigvec(j,i)
-  enddo
-  print*,i
-  write(*,'(I3,X,100(F16.10,X))')i,eigval(i),eigval_ref(i),accu_l,accu_r
- enddo
-end

plugins/local/tc_bi_ortho/test_normal_order.irp.f

@@ -1,173 +0,0 @@
-
-! ---
-
-program test_normal_order
-
-  BEGIN_DOC
-  ! TODO : Put the documentation of the program here
-  END_DOC
-
-  implicit none
-
-  print *, 'Hello world'
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  call provide_all_three_ints_bi_ortho()
-  call test()
-
-end
-
-! ---
-
-subroutine test
- implicit none
-  use bitmasks ! you need to include the bitmasks_module.f90 features
- integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2)
- integer                        :: exc(0:2,2,2)
- integer(bit_kind), allocatable :: det_i(:,:)
- double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal,hthree_tmp
- integer,           allocatable :: occ(:,:)
- allocate( occ(N_int*bit_kind_size,2) )
- call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
- allocate(det_i(N_int,2))
- s1 = 1 
- s2 = 2
- accu = 0.d0
- do h1 = 1, elec_beta_num
-  do p1 = elec_alpha_num+1, mo_num
-   do h2 = 1, elec_beta_num
-    do p2 = elec_beta_num+1, mo_num
-     hthree = 0.d0
-
-     det_i = ref_bitmask
- s1 = 1 
- s2 = 2
-     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_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,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_tmp *= phase
-     hthree += 0.5d0 * hthree_tmp
-     det_i = ref_bitmask
- s1 = 2 
- s2 = 1
-     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_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,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_tmp *= phase
-     hthree += 0.5d0 * hthree_tmp
-
-
-!     normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
-     call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, normal)
-     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
-!     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
-!    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_slow(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)
-!     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_slow(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
-
-

plugins/local/tc_bi_ortho/test_s2_tc.irp.f

@@ -1,170 +0,0 @@
-
-! ---
-
-program test_tc 
-
-  implicit none
- 
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  call provide_all_three_ints_bi_ortho()
-  call routine_h_triple_left
-  call routine_h_triple_right
-!  call routine_test_s2_davidson
-
-end
-
-subroutine routine_h_triple_right
- implicit none
- logical           :: do_right 
- integer :: sze ,i, N_st, j
- double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
- double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
- double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
- sze = N_det
- N_st = 1
- allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
- print*,'Checking first the Right '
- do i = 1, sze
-  u_0(i,1) = psi_r_coef_bi_ortho(i,1)
- enddo
- double precision :: wall0,wall1
- call wall_time(wall0)
- call H_tc_s2_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
- call wall_time(wall1)
- print*,'time for omp',wall1 - wall0
- call wall_time(wall0)
- call H_tc_s2_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
- call wall_time(wall1)
- print*,'time serial ',wall1 - wall0
- accu_e = 0.d0
- accu_s = 0.d0
- do i = 1, sze
-  accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
-  accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
- enddo
- print*,'accu_e   = ',accu_e
- print*,'accu_s   = ',accu_s
-
-end
-
-subroutine routine_h_triple_left
- implicit none
- logical           :: do_right 
- integer :: sze ,i, N_st, j
- double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
- double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
- double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
- sze = N_det
- N_st = 1
- allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
- print*,'Checking the Left '
- do i = 1, sze
-  u_0(i,1) = psi_l_coef_bi_ortho(i,1)
- enddo
- double precision :: wall0,wall1
- call wall_time(wall0)
- call H_tc_s2_dagger_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
- call wall_time(wall1)
- print*,'time for omp',wall1 - wall0
- call wall_time(wall0)
- call H_tc_s2_dagger_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
- call wall_time(wall1)
- print*,'time serial ',wall1 - wall0
- accu_e = 0.d0
- accu_s = 0.d0
- do i = 1, sze
-  accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
-  accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
- enddo
- print*,'accu_e   = ',accu_e
- print*,'accu_s   = ',accu_s
-
-end
-
-
-subroutine routine_test_s2_davidson
- implicit none
- double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:)
- integer :: i,istate
- logical :: converged 
- external H_tc_s2_dagger_u_0_opt
- external H_tc_s2_u_0_opt
- allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag))
- do i = 1, N_det
-   call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
- enddo
- ! Preparing the left-eigenvector
- print*,'Computing the left-eigenvector '
- vec_tmp = 0.d0
- do istate = 1, N_states
-  vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
- enddo
- do istate = N_states+1, n_states_diag
-  vec_tmp(istate,istate) = 1.d0
- enddo
- do istate = 1, N_states
-  leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
- enddo
- integer :: n_it_max
- n_it_max = 1
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
- double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
- integer :: sze,N_st
- logical           :: do_right 
- sze = N_det
- N_st = 1
- do_right = .False.
- allocate(s_0_new(N_det,1),v_0_new(N_det,1))
- call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
- double precision :: accu_e_0, accu_s_0
- accu_e_0 = 0.d0
- accu_s_0 = 0.d0
- do i = 1, sze
-  accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
-  accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
- enddo
- print*,'energies = ',energies
- print*,'s2       = ',s2
- print*,'accu_e_0',accu_e_0
- print*,'accu_s_0',accu_s_0
-
- ! Preparing the right-eigenvector
- print*,'Computing the right-eigenvector '
- vec_tmp = 0.d0
- do istate = 1, N_states
-  vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
- enddo
- do istate = N_states+1, n_states_diag
-  vec_tmp(istate,istate) = 1.d0
- enddo
- do istate = 1, N_states
-  leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
- enddo
- n_it_max = 1
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
- sze = N_det
- N_st = 1
- do_right = .True.
- v_0_new = 0.d0
- s_0_new = 0.d0
- call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
- accu_e_0 = 0.d0
- accu_s_0 = 0.d0
- do i = 1, sze
-  accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
-  accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
- enddo
- print*,'energies = ',energies
- print*,'s2       = ',s2
- print*,'accu_e_0',accu_e_0
- print*,'accu_s_0',accu_s_0
-
-end

plugins/local/tc_bi_ortho/test_tc_fock.irp.f

@@ -1,171 +0,0 @@
-
-! ---
-
-program test_tc_fock
-
-  BEGIN_DOC
-  ! TODO : Put the documentation of the program here
-  END_DOC
-
-  implicit none
-
-  print *, 'Hello world'
-
-  my_grid_becke = .True.
-  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = tc_grid1_a
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
-  read_wf = .True.
-  touch read_wf
-
-  !call routine_1
-  !call routine_2
-!  call routine_3()
-
-  call routine_tot
-
-end
-
-! ---
-
-subroutine routine_3()
-
-  use bitmasks ! you need to include the bitmasks_module.f90 features
-
-  implicit none
-  integer                        :: i, a, i_ok, s1
-  double precision               :: hmono, htwoe, hthree, htilde_ij
-  double precision               :: err_ai, err_tot, ref, new
-  integer(bit_kind), allocatable :: det_i(:,:)
-
-  allocate(det_i(N_int,2))
-
-  err_tot = 0.d0
- 
-  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_num_tab(s1)
-     do a = elec_num_tab(s1)+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_slow(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
-       print*, ref,new,err_ai
-       err_tot += err_ai
- 
-       write(22, *) htilde_ij
-     enddo
-   enddo
-  enddo
-
-  print *, ' err_tot = ', err_tot
-
-  deallocate(det_i)
-
-end subroutine routine_3
-
-! ---
-subroutine routine_tot()
-
-  use bitmasks ! you need to include the bitmasks_module.f90 features
-
-  implicit none
-  integer                        :: i, a, i_ok, s1,other_spin(2)
-  double precision               :: hmono, htwoe, hthree, htilde_ij
-  double precision               :: err_ai, err_tot, ref, new
-  integer(bit_kind), allocatable :: det_i(:,:)
-
-  allocate(det_i(N_int,2))
-  other_spin(1) = 2
-  other_spin(2) = 1
-
-  err_tot = 0.d0
- 
-!  do s1 = 1, 2
-   s1 = 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_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, mo_num! virtual 
-       print*,i,a 
- 
-       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
- 
-       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
-       print*, ' excited det'
-       call debug_det(det_i, N_int)
-
-       if(s1 == 1)then
-        new = Fock_matrix_tc_mo_alpha(a,i)
-       else
-        new = Fock_matrix_tc_mo_beta(a,i)
-       endif
-       ref = htilde_ij
-!       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*,'---------'
-         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
- 
-       write(22, *) htilde_ij
-     enddo
-   enddo
-!  enddo
-
-  print *, ' err_tot = ', err_tot
-
-  deallocate(det_i)
-
-end subroutine routine_3

plugins/local/tc_keywords/EZFIO.cfg

@@ -184,12 +184,6 @@ doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
-[debug_tc_pt2]
-type: integer
-doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
-interface: ezfio,provider,ocaml
-default: -1
-
 [only_spin_tc_right]
 type: logical
 doc: If |true|, only the right part of WF is used to compute spin dens
@@ -268,3 +262,8 @@ doc: Thresholds on the Imag part of TC energy
 interface: ezfio,provider,ocaml
 default: 1.e-7
 
+[transpose_two_e_int]
+type: logical
+doc: If |true|, you duplicate the two-electron TC integrals with the transpose matrix. Acceleates the PT2.  
+interface: ezfio,provider,ocaml
+default: False

plugins/local/tc_progs/NEED

@@ -0,0 +1 @@
+tc_bi_ortho

plugins/local/tc_progs/print_tc_wf.irp.f

@@ -61,12 +61,12 @@ subroutine routine
  do i = 1, N_det
   call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
    if(degree == 1 .or. degree == 2)then
-    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
-    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
     delta_e = e_tilde_00 - e_i0
     coef_pt1 = htilde_ij / delta_e
  
-    call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
+    call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
     contrib_pt = coef_pt1 * htilde_ij
     e_pt2 += contrib_pt
  

plugins/local/tc_progs/tc_som.irp.f

@@ -49,8 +49,8 @@ subroutine main()
   U_SOM = 0.d0 
   do i = 1, N_det
     if(i == i_HF) cycle
-    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i_HF), psi_det(1,1,i), N_int, hmono_1, htwoe_1, hthree_1, htot_1)
-    call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,i_HF), N_int, hmono_2, htwoe_2, hthree_2, htot_2)
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i_HF), psi_det(1,1,i), N_int, hmono_1, htwoe_1, hthree_1, htot_1)
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,i_HF), N_int, hmono_2, htwoe_2, hthree_2, htot_2)
     U_SOM += htot_1 * htot_2
   enddo
   U_SOM = 0.5d0 * U_SOM

src/determinants/slater_rules_general.irp.f

@@ -0,0 +1,192 @@
+subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
+ use bitmasks
+ BEGIN_DOC
+! returns the array, for each spin, of holes/particles between key_i and key_j 
+!
+! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
+ END_DOC
+  include 'utils/constants.include.F'
+ implicit none
+ integer, intent(in)            :: Nint
+ integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
+ integer, intent(out)           :: holes_array(100,2),particles_array(100,2),degree_array(2)
+ double precision, intent(out)  :: phase
+ integer :: ispin,k,i,pos
+ integer(bit_kind) :: key_hole, key_particle
+ integer(bit_kind) :: xorvec(N_int_max,2)
+ holes_array = -1
+ particles_array = -1
+ degree_array = 0
+  do i = 1, N_int
+   xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
+   xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
+   degree_array(1) += popcnt(xorvec(i,1))
+   degree_array(2) += popcnt(xorvec(i,2))
+  enddo
+  degree_array(1) = shiftr(degree_array(1),1)
+  degree_array(2) = shiftr(degree_array(2),1)
+  
+ do ispin = 1, 2
+  k = 1
+  !!! GETTING THE HOLES 
+  do i = 1, N_int
+   key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
+   do while(key_hole .ne.0_bit_kind)
+    pos = trailz(key_hole)
+    holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
+    key_hole = ibclr(key_hole,pos)
+    k += 1
+    if(k .gt.100)then
+     print*,'WARNING in get_excitation_general'
+     print*,'More than a 100-th excitation for spin ',ispin
+     print*,'stoping ...'
+     stop
+    endif
+   enddo 
+  enddo
+ enddo
+ do ispin = 1, 2
+  k = 1
+  !!! GETTING THE PARTICLES
+  do i = 1, N_int
+   key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
+   do while(key_particle .ne.0_bit_kind)
+    pos = trailz(key_particle)
+    particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
+    key_particle = ibclr(key_particle,pos)
+    k += 1
+    if(k .gt.100)then
+     print*,'WARNING in get_excitation_general '
+     print*,'More than a 100-th excitation for spin ',ispin
+     print*,'stoping ...'
+     stop
+    endif
+   enddo
+  enddo 
+ enddo
+ integer :: h,p, i_ok
+ integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
+ integer                        :: exc(0:2,2,2)
+ double precision :: phase_tmp
+ allocate(det_i(Nint,2),det_ip(N_int,2))
+ det_i = key_i
+ phase = 1.d0
+ do ispin = 1, 2
+  do i = 1, degree_array(ispin)
+   h = holes_array(i,ispin)
+   p = particles_array(i,ispin)
+   det_ip = det_i
+   call do_single_excitation(det_ip,h,p,ispin,i_ok)
+   if(i_ok == -1)then
+     print*,'excitation was not possible '
+     stop
+   endif
+   call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
+   phase *= phase_tmp
+   det_i = det_ip
+  enddo
+ enddo
+
+end
+
+subroutine get_holes_general(key_i, key_j,Nint, holes_array)
+ use bitmasks
+ BEGIN_DOC
+! returns the array, per spin, of holes between key_i and key_j 
+!
+! with the following convention: a_{hole}|key_i> --> |key_j>
+ END_DOC
+ implicit none
+ integer, intent(in)            :: Nint
+ integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
+ integer, intent(out)           :: holes_array(100,2)
+ integer(bit_kind) :: key_hole
+ integer :: ispin,k,i,pos
+ holes_array = -1
+ do ispin = 1, 2
+  k = 1
+  do i = 1, N_int
+   key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
+   do while(key_hole .ne.0_bit_kind)
+    pos = trailz(key_hole)
+    holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
+    key_hole = ibclr(key_hole,pos)
+    k += 1
+    if(k .gt.100)then
+     print*,'WARNING in get_holes_general'
+     print*,'More than a 100-th excitation for spin ',ispin
+     print*,'stoping ...'
+     stop
+    endif
+   enddo 
+  enddo
+ enddo
+end
+
+subroutine get_particles_general(key_i, key_j,Nint,particles_array)
+ use bitmasks
+ BEGIN_DOC
+! returns the array, per spin, of particles between key_i and key_j 
+!
+! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
+ END_DOC
+ implicit none
+ integer, intent(in)            :: Nint
+ integer(bit_kind), intent(in)  :: key_j(Nint,2),key_i(Nint,2)
+ integer, intent(out)           :: particles_array(100,2)
+ integer(bit_kind) :: key_particle
+ integer :: ispin,k,i,pos
+ particles_array = -1
+ do ispin = 1, 2
+  k = 1
+  do i = 1, N_int
+   key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
+   do while(key_particle .ne.0_bit_kind)
+    pos = trailz(key_particle)
+    particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
+    key_particle = ibclr(key_particle,pos)
+    k += 1
+    if(k .gt.100)then
+     print*,'WARNING in get_holes_general'
+     print*,'More than a 100-th excitation for spin ',ispin
+     print*,'Those are the two determinants'
+     call debug_det(key_i, N_int)
+     call debug_det(key_j, N_int)
+     print*,'stoping ...'
+     stop
+    endif
+   enddo 
+  enddo
+ enddo
+end
+
+subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
+ implicit none
+ integer, intent(in)            :: degree(2), Nint
+ integer(bit_kind), intent(in)  :: key_i(Nint,2)
+ integer, intent(in)            :: holes_array(100,2),particles_array(100,2)
+ double precision, intent(out)  :: phase
+ integer :: i,ispin,h,p, i_ok
+ integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
+ integer                        :: exc(0:2,2,2)
+ double precision :: phase_tmp
+ allocate(det_i(Nint,2),det_ip(N_int,2))
+ det_i = key_i
+ phase = 1.d0
+ do ispin = 1, 2
+  do i = 1, degree(ispin)
+   h = holes_array(i,ispin)
+   p = particles_array(i,ispin)
+   det_ip = det_i
+   call do_single_excitation(det_ip,h,p,ispin,i_ok)
+   if(i_ok == -1)then
+     print*,'excitation was not possible '
+     stop
+   endif
+   call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
+   phase *= phase_tmp
+   det_i = det_ip
+  enddo
+ enddo
+
+end