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

src/casscf_tc_bi/NEED

@@ -0,0 +1,3 @@
+determinants
+tc_bi_ortho
+fci_tc_bi

src/casscf_tc_bi/det_manip.irp.f

@@ -0,0 +1,125 @@
+use bitmasks
+
+subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart,       &
+      ispin,phase,ierr)
+  BEGIN_DOC
+  ! we create the mono-excitation, and determine, if possible,
+  ! the phase and the number in the list of determinants
+  END_DOC
+  implicit none
+  integer(bit_kind)              :: key1(N_int,2),key2(N_int,2)
+  integer(bit_kind), allocatable :: keytmp(:,:)
+  integer                        :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin
+  real*8                         :: phase
+  logical                        :: found
+  allocate(keytmp(N_int,2))
+  
+  nu=-1
+  phase=1.D0
+  ierr=0
+  call det_copy(key1,key2,N_int)
+  !        write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin
+  !        call print_det(key2,N_int)
+  call do_single_excitation(key2,ihole,ipart,ispin,ierr)
+  !        write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin
+  !        call print_det(key2,N_int)
+  !        write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr
+  if (ierr.eq.1) then
+    ! excitation is possible
+    ! get the phase
+    call get_single_excitation(key1,key2,exc,phase,N_int)
+    ! get the number in the list
+    found=.false.
+    nu=0
+
+    !TODO BOTTLENECK
+    do while (.not.found)
+      nu+=1
+      if (nu.gt.N_det) then
+        ! the determinant is possible, but not in the list
+        found=.true.
+        nu=-1
+      else
+        call det_extract(keytmp,nu,N_int)
+        integer                        :: i,ii
+        found=.true.
+        do ii=1,2
+          do i=1,N_int
+            if (keytmp(i,ii).ne.key2(i,ii)) then
+              found=.false.
+            end if
+          end do
+        end do
+      end if
+    end do
+  end if
+  !
+  ! we found the new string, the phase, and possibly the number in the list
+  !
+end subroutine do_signed_mono_excitation
+
+subroutine det_extract(key,nu,Nint)
+  BEGIN_DOC
+  ! extract a determinant from the list of determinants
+  END_DOC
+  implicit none
+  integer                        :: ispin,i,nu,Nint
+  integer(bit_kind)              :: key(Nint,2)
+  do ispin=1,2
+    do i=1,Nint
+      key(i,ispin)=psi_det(i,ispin,nu)
+    end do
+  end do
+end subroutine det_extract
+
+subroutine det_copy(key1,key2,Nint)
+  use bitmasks ! you need to include the bitmasks_module.f90 features
+  BEGIN_DOC
+  ! copy a determinant from key1 to key2
+  END_DOC
+  implicit none
+  integer                        :: ispin,i,Nint
+  integer(bit_kind)              :: key1(Nint,2),key2(Nint,2)
+  do ispin=1,2
+    do i=1,Nint
+      key2(i,ispin)=key1(i,ispin)
+    end do
+  end do
+end subroutine det_copy
+
+subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2      &
+      ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr)
+  BEGIN_DOC
+  ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q)
+  ! we may create two determinants as result
+  !
+  END_DOC
+  implicit none
+  integer(bit_kind)              :: key_in(N_int,2),key_out1(N_int,2)
+  integer(bit_kind)              :: key_out2(N_int,2)
+  integer                        :: ihole,ipart,ierr,jerr,nu1,nu2
+  integer                        :: ispin
+  real*8                         :: phase1,phase2
+  
+  !        write(6,*) ' applying E_',ipart,ihole,' on determinant '
+  !        call print_det(key_in,N_int)
+  
+  ! spin alpha
+  ispin=1
+  call do_signed_mono_excitation(key_in,key_out1,nu1,ihole           &
+      ,ipart,ispin,phase1,ierr)
+  !        if (ierr.eq.1) then
+  !         write(6,*) ' 1 result is ',nu1,phase1
+  !         call print_det(key_out1,N_int)
+  !        end if
+  ! spin beta
+  ispin=2
+  call do_signed_mono_excitation(key_in,key_out2,nu2,ihole           &
+      ,ipart,ispin,phase2,jerr)
+  !        if (jerr.eq.1) then
+  !         write(6,*) ' 2 result is ',nu2,phase2
+  !         call print_det(key_out2,N_int)
+  !        end if
+  
+end subroutine do_spinfree_mono_excitation
+

src/casscf_tc_bi/grad_dm.irp.f

@@ -193,7 +193,7 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
    r = list_act(rr)
    res_r_inact_test += -tc_transition_matrix_mo(r,t,1,1) * & 
    (2.d0 * mo_bi_ortho_tc_two_e(r,j,a,j) - mo_bi_ortho_tc_two_e(r,j,j,a))
-   res_l_inact_test += -tc_transition_matrix_mo(t,r,1,1) * & 
+   res_l_inact_test -= -tc_transition_matrix_mo(t,r,1,1) * & 
    (2.d0 * mo_bi_ortho_tc_two_e(a,j,r,j) - mo_bi_ortho_tc_two_e(j,a,r,j))
   enddo
  enddo
@@ -207,7 +207,7 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
     u = list_act(uu)
     res_r_act_test += - (mo_bi_ortho_tc_two_e(v,r,u,a) * tc_two_rdm(r,v,t,u) & 
                         +mo_bi_ortho_tc_two_e(v,r,a,u) * tc_two_rdm(r,v,u,t))
-    res_l_act_test += - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) & 
+    res_l_act_test -= - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) & 
                         +mo_bi_ortho_tc_two_e(a,u,v,r) * tc_two_rdm(u,t,r,v))
    enddo
   enddo

src/casscf_tc_bi/test_tc_casscf.irp.f

@@ -0,0 +1,252 @@
+program tc_bi_ortho
+
+  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
+
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  read_wf = .True.
+  touch read_wf
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  print*, ' nb of states = ', N_states
+  print*, ' nb of det    = ', N_det
+
+!  call routine_i_h_psi
+!  call routine_grad
+ call routine_grad_num_dm_one_body
+end
+
+subroutine routine_i_h_psi
+ implicit none
+ integer :: i,j
+ double precision :: i_H_chi_array(0:3,N_states),i_H_phi_array(0:3,N_states)
+ double precision               :: hmono, htwoe, hthree, htot
+ double precision               :: accu_l_hmono, accu_l_htwoe, accu_l_hthree, accu_l_htot
+ double precision               :: accu_r_hmono, accu_r_htwoe, accu_r_hthree, accu_r_htot
+ double precision               :: test_l_hmono, test_l_htwoe, test_l_hthree, test_l_htot
+ double precision               :: test_r_hmono, test_r_htwoe, test_r_hthree, test_r_htot
+
+ test_l_hmono = 0.d0
+ test_l_htwoe = 0.d0
+ test_l_hthree= 0.d0
+ test_l_htot  = 0.d0
+ test_r_hmono = 0.d0
+ test_r_htwoe = 0.d0
+ test_r_hthree= 0.d0
+ test_r_htot  = 0.d0
+
+ do i = 1, N_det
+   call i_H_tc_psi_phi(psi_det(1,1,i),psi_det,psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,&
+                       N_int,N_det,N_det,N_states,i_H_chi_array,i_H_phi_array)
+   accu_l_hmono = 0.d0
+   accu_l_htwoe = 0.d0
+   accu_l_hthree= 0.d0
+   accu_l_htot  = 0.d0
+   accu_r_hmono = 0.d0
+   accu_r_htwoe = 0.d0
+   accu_r_hthree= 0.d0
+   accu_r_htot  = 0.d0
+   do j = 1, N_det
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
+    accu_l_hmono  += psi_l_coef_bi_ortho(j,1) * hmono
+    accu_l_htwoe  += psi_l_coef_bi_ortho(j,1) * htwoe
+    accu_l_hthree += psi_l_coef_bi_ortho(j,1) * hthree
+    accu_l_htot   += psi_l_coef_bi_ortho(j,1) * htot
+    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+    accu_r_hmono  += psi_r_coef_bi_ortho(j,1) * hmono
+    accu_r_htwoe  += psi_r_coef_bi_ortho(j,1) * htwoe
+    accu_r_hthree += psi_r_coef_bi_ortho(j,1) * hthree
+    accu_r_htot   += psi_r_coef_bi_ortho(j,1) * htot
+   enddo
+   test_l_htot   += dabs(i_H_chi_array(0,1)-accu_l_htot)
+   test_l_hmono  += dabs(i_H_chi_array(1,1)-accu_l_hmono)
+   test_l_htwoe  += dabs(i_H_chi_array(2,1)-accu_l_htwoe)
+   test_l_hthree += dabs(i_H_chi_array(3,1)-accu_l_hthree)
+
+   test_r_htot   += dabs(i_H_phi_array(0,1)-accu_r_htot)
+   test_r_hmono  += dabs(i_H_phi_array(1,1)-accu_r_hmono)
+   test_r_htwoe  += dabs(i_H_phi_array(2,1)-accu_r_htwoe)
+   test_r_hthree += dabs(i_H_phi_array(3,1)-accu_r_hthree)
+
+ enddo
+ 
+ test_l_htot   *= 1.D0/dble(N_det)
+ test_l_hmono  *= 1.D0/dble(N_det)
+ test_l_htwoe  *= 1.D0/dble(N_det)
+ test_l_hthree *= 1.D0/dble(N_det)
+
+ test_r_htot   *= 1.D0/dble(N_det)
+ test_r_hmono  *= 1.D0/dble(N_det)
+ test_r_htwoe  *= 1.D0/dble(N_det)
+ test_r_hthree *= 1.D0/dble(N_det)
+ 
+ print*,'**************************'
+ print*,'test_l_htot    = ',test_l_htot
+ print*,'test_l_hmono   = ',test_l_hmono
+ print*,'test_l_htwoe   = ',test_l_htwoe
+ print*,'test_l_hthree  = ',test_l_hthree
+ print*,'**************************'
+ print*,'test_r_htot    = ',test_r_htot
+ print*,'test_r_hmono   = ',test_r_hmono
+ print*,'test_r_htwoe   = ',test_r_htwoe
+ print*,'test_r_hthree  = ',test_r_hthree
+
+end
+
+subroutine routine_grad_num
+ implicit none
+ integer :: indx,ihole,ipart
+ integer :: p,q
+ double precision :: accu_l, accu_r
+ double precision :: contrib_l, contrib_r
+
+ accu_l = 0.d0
+ accu_r = 0.d0
+  do indx=1,nMonoEx
+   q = excit(1,indx)
+   p = excit(2,indx)
+   contrib_l  =  dabs(dabs(gradvec_detail_left_old(0,indx)) - 2.D0 * dabs( Fock_matrix_tc_mo_tot(q,p)))
+   contrib_r  =  dabs(dabs(gradvec_detail_right_old(0,indx)) -2.D0 * dabs( Fock_matrix_tc_mo_tot(p,q)))
+   if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+    print*,indx,q,p
+    print*,gradvec_detail_left_old(0,indx),gradvec_detail_right_old(0,indx) 
+    print*,2.D0* Fock_matrix_tc_mo_tot(q,p), 2.d0* Fock_matrix_tc_mo_tot(p,q)
+   endif
+   accu_l += contrib_l
+   accu_r += contrib_r
+  enddo
+ print*,'accu_l,accu_r'
+ print*,accu_l,accu_r
+ 
+! do i = 1, nMonoEx
+!  print*,i,gradvec_old(i)
+! enddo
+
+end
+
+subroutine routine_grad_num_dm_one_body
+ implicit none
+ integer :: indx,ii,i,a,aa,tt,t,ibody
+ double precision :: accu_l, accu_r,ref_r, new_r, ref_l, new_l
+ double precision :: contrib_l, contrib_r
+ double precision :: res_l(0:3),res_r(0:3)
+
+ ibody = 2 ! check only the two-body term
+ provide gradvec_detail_left_old gradvec_tc_l 
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing inactive-->virtual'
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do ii = 1, n_core_inact_orb
+   do aa = 1, n_virt_orb
+    indx = mat_idx_c_v(ii,aa) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    i = list_core_inact(ii)
+    a = list_virt(aa)
+!    if(i==1.and.a==9)then
+!     print*,i,a,ref_r, new_r
+!     stop
+!    endif
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     print*,indx,i,a,ii,aa
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+     print*,gradvec_detail_left_old(0,indx),gradvec_tc_l(0,indx)
+     print*,gradvec_detail_right_old(0,indx),gradvec_tc_r(0,indx)
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+  print*,'**************************'
+  print*,'**************************'
+ endif
+ 
+ ibody = 2 ! check only the two-body term
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing inactive-->active'
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do ii = 1, n_core_inact_orb
+   do tt = 1, n_act_orb
+    indx = mat_idx_c_a(ii,tt) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     i = list_core_inact(ii)
+     t = list_act(tt)
+     print*,indx,i,t
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+ endif
+
+ if(.True.)then
+  print*,'**************************'
+  print*,'**************************'
+  print*,'testing active-->virtual '
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do tt = 1, n_act_orb
+   do aa = 1, n_virt_orb
+    indx = mat_idx_a_v(tt,aa) 
+    ref_l = gradvec_detail_left_old(ibody,indx)
+    new_l = gradvec_tc_l(ibody,indx) 
+    contrib_l  =  dabs(dabs(ref_l) - dabs(new_l))
+    ref_r = gradvec_detail_right_old(ibody,indx)
+    new_r = gradvec_tc_r(ibody,indx) 
+    contrib_r  =  dabs(dabs(ref_r) - dabs(new_r))
+    if(contrib_l.gt.1.d-10.or.contrib_r.gt.1.d-10)then
+     print*,'---------'
+     print*,'warning !'
+     a = list_virt(aa)
+     t = list_act(tt)
+     print*,indx,t,a
+     print*,ref_l, new_l, contrib_l
+     print*,ref_r, new_r, contrib_r
+!     print*,gradvec_detail_right_old(0,indx),gradvec_tc_r(0,indx)
+     print*,'---------'
+    endif
+    accu_l += contrib_l
+    accu_r += contrib_r
+   enddo
+  enddo
+  print*,'accu_l,accu_r'
+  print*,accu_l,accu_r
+ endif
+ 
+
+end

src/cisd/lccsd_prov.irp.f

@@ -3,7 +3,7 @@
  implicit none
  double precision, allocatable :: Dress_jj(:), H_jj(:), u_in(:,:)
  double precision :: ebefore, eafter, ecorr, thresh
- integer :: i,it
+ integer :: i,it,degree
  logical :: converged
  external H_u_0_nstates_openmp
  allocate(Dress_jj(N_det),H_jj(N_det),u_in(N_det,N_states_diag))
@@ -31,7 +31,11 @@
   print*,'ecorr = ',ecorr
   Dress_jj(1) = 0.d0
   do i = 2, N_det
-    if(ecorr + H_jj(i) .gt. H_jj(1))then
+    if(ecorr + H_jj(i) .lt. H_jj(1))then
+     print*,'Warning, some dets are not dressed: ' 
+     call get_excitation_degree(ref_bitmask,psi_det(1,1,i),degree,N_int)
+     print*,'degree, Delta E, coef', degree, H_jj(i)-H_jj(1), u_in(i,1)/u_in(1,1)
+    else
      Dress_jj(i) = ecorr
     endif
   enddo

src/tc_bi_ortho/h_mat_triple.irp.f

@@ -0,0 +1,391 @@
+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 
+  !
+  ! Assumes that the determinants are in psi_det
+  !
+  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+  END_DOC
+
+  use bitmasks
+  implicit none
+
+  integer,          intent(in)  :: N_st,sze
+  double precision, intent(in)  :: u_0(sze,N_st)
+  double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
+  call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
+  integer :: i,j,degree,ist
+  double precision :: hmono, htwoe, hthree, htot
+  do i = 1, N_det
+   do j = 1, N_det
+    call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+    if(degree .ne. 3)cycle
+    call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
+    do ist = 1, N_st
+     v_0(i,ist) += htot * u_0(j,ist)
+    enddo
+   enddo
+  enddo
+end
+
+subroutine H_tc_s2_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
+  BEGIN_DOC
+  ! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS 
+  !
+  ! Assumes that the determinants are in psi_det
+  !
+  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+  END_DOC
+
+  use bitmasks
+  implicit none
+
+  integer,          intent(in)  :: N_st,sze
+  double precision, intent(in)  :: u_0(sze,N_st)
+  double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
+  call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
+  integer :: i,j,degree,ist
+  double precision :: hmono, htwoe, hthree, htot
+  !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
+  !$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0)       &
+  !$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
+  do i = 1, N_det
+   do j = 1, N_det
+    call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+    if(degree .ne. 3)cycle
+    call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
+    do ist = 1, N_st
+     v_0(i,ist) += htot * u_0(j,ist)
+    enddo
+   enddo
+  enddo
+ !$OMP END PARALLEL DO
+end
+
+! ---
+
+subroutine H_tc_s2_dagger_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
+  BEGIN_DOC
+  ! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS 
+  !
+  ! Assumes that the determinants are in psi_det
+  !
+  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+  END_DOC
+
+  use bitmasks
+  implicit none
+
+  integer,          intent(in)  :: N_st,sze
+  double precision, intent(in)  :: u_0(sze,N_st)
+  double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
+  call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
+  integer :: i,j,degree,ist
+  double precision :: hmono, htwoe, hthree, htot
+  do i = 1, N_det
+   do j = 1, N_det
+    call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+    if(degree .ne. 3)cycle
+    call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
+    do ist = 1, N_st
+     v_0(i,ist) += htot * u_0(j,ist)
+    enddo
+   enddo
+  enddo
+end
+
+subroutine H_tc_s2_dagger_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
+  BEGIN_DOC
+  ! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS 
+  !
+  ! Assumes that the determinants are in psi_det
+  !
+  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+  END_DOC
+
+  use bitmasks
+  implicit none
+
+  integer,          intent(in)  :: N_st,sze
+  double precision, intent(in)  :: u_0(sze,N_st)
+  double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
+  call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
+  integer :: i,j,degree,ist
+  double precision :: hmono, htwoe, hthree, htot
+  !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
+  !$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0)       &
+  !$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
+  do i = 1, N_det
+   do j = 1, N_det
+    call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
+    if(degree .ne. 3)cycle
+    call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
+    do ist = 1, N_st
+     v_0(i,ist) += htot * u_0(j,ist)
+    enddo
+   enddo
+  enddo
+ !$OMP END PARALLEL DO
+end
+
+! ---
+subroutine triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+  use bitmasks
+  BEGIN_DOC
+! <key_j | H_tilde | key_i> for triple excitation  
+!!
+!! WARNING !!
+! 
+! Genuine triple excitations of the same spin are not yet implemented
+  END_DOC
+  implicit none
+  integer(bit_kind), intent(in)  :: key_j(N_int,2),key_i(N_int,2)
+  integer, intent(in)            :: Nint
+  double precision, intent(out)  :: hmono, htwoe, hthree, htot
+  integer                        :: degree
+  integer                        :: h1, p1, h2, p2, s1, s2, h3, p3, s3
+  integer                        :: holes_array(100,2),particles_array(100,2),degree_array(2)
+  double precision               :: phase,sym_3_e_int_from_6_idx_tensor
+
+  hmono  = 0.d0
+  htwoe  = 0.d0
+  hthree = 0.d0
+  htot   = 0.d0
+  call get_excitation_general(key_j, key_i, Nint,degree_array,holes_array, particles_array,phase)
+  degree = degree_array(1) + degree_array(2)
+  if(degree .ne. 3)return
+  if(degree_array(1)==3.or.degree_array(2)==3)then
+   if(degree_array(1) == 3)then
+    h1 = holes_array(1,1) 
+    h2 = holes_array(2,1) 
+    h3 = holes_array(3,1) 
+    p1 = particles_array(1,1) 
+    p2 = particles_array(2,1) 
+    p3 = particles_array(3,1) 
+   else
+    h1 = holes_array(1,2) 
+    h2 = holes_array(2,2) 
+    h3 = holes_array(3,2) 
+    p1 = particles_array(1,2) 
+    p2 = particles_array(2,2) 
+    p3 = particles_array(3,2) 
+   endif
+   hthree = sym_3_e_int_from_6_idx_tensor(p3, p2, p1, h3, h2, h1)
+  else 
+   if(degree_array(1) == 2.and.degree_array(2) == 1)then ! double alpha + single beta
+    h1 = holes_array(1,1) 
+    h2 = holes_array(2,1) 
+    h3 = holes_array(1,2) 
+    p1 = particles_array(1,1) 
+    p2 = particles_array(2,1) 
+    p3 = particles_array(1,2) 
+   else if(degree_array(2) == 2 .and. degree_array(1) == 1)then ! double beta + single alpha 
+    h1 = holes_array(1,2) 
+    h2 = holes_array(2,2) 
+    h3 = holes_array(1,1) 
+    p1 = particles_array(1,2) 
+    p2 = particles_array(2,2) 
+    p3 = particles_array(1,1) 
+   else 
+    print*,'PB !!'
+    stop
+   endif
+   hthree = three_body_ints_bi_ort(p3,p2,p1,h3,h2,h1) - three_body_ints_bi_ort(p3,p2,p1,h3,h1,h2)
+  endif
+  hthree  *= phase
+  htot = hthree
+ end
+

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -19,6 +19,9 @@ subroutine provide_all_three_ints_bi_ortho()
     if(three_e_4_idx_term) then
       PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort 
     endif
+    if(pure_three_body_h_tc)then
+     provide three_body_ints_bi_ort
+    endif
 
     if(.not. double_normal_ord .and. three_e_5_idx_term) then
       PROVIDE three_e_5_idx_direct_bi_ort 
@@ -85,14 +88,26 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
   hthree = 0.D0
 
   call get_excitation_degree(key_i, key_j, degree, Nint)
-  if(degree.gt.2) return
+  if(.not.pure_three_body_h_tc)then
-
+   if(degree.gt.2) return
-  if(degree == 0) then
+   if(degree == 0) then
-    call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
+     call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
-  else if (degree == 1) then
+   else if (degree == 1) then
-    call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
+     call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
-  else if(degree == 2) then
+   else if(degree == 2) then
-    call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+     call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+   endif
+  else 
+   if(degree.gt.3) return
+   if(degree == 0) then
+     call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
+   else if (degree == 1) then
+     call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
+   else if(degree == 2) then
+     call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+   else
+     call triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+   endif
   endif
 
   if(degree==0) then

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -225,6 +225,8 @@ end
     external                         H_tc_dagger_u_0_opt
     external                         H_tc_s2_dagger_u_0_opt
     external                         H_tc_s2_u_0_opt
+    external                         H_tc_s2_dagger_u_0_with_pure_three_omp
+    external                         H_tc_s2_u_0_with_pure_three_omp
 
     allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
 
@@ -250,7 +252,11 @@ end
     converged = .False.
     i_it = 0
     do while (.not.converged)
-      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+      if(.not.pure_three_body_h_tc)then
+       call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
+      else 
+       call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_with_pure_three_omp)
+      endif
       i_it += 1
       if(i_it .gt. 5) exit
     enddo
@@ -275,7 +281,11 @@ end
     converged = .False.
     i_it = 0
     do while (.not. converged)
-      call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
+      if(.not.pure_three_body_h_tc)then
+       call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
+      else
+       call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_with_pure_three_omp)
+      endif
       i_it += 1
       if(i_it .gt. 5) exit
     enddo

src/tc_bi_ortho/tc_natorb.irp.f

@@ -32,19 +32,17 @@
   thr_d   = 1.d-6
   thr_nd  = 1.d-6
   thr_deg = 1.d-3
-  if(n_core_orb.ne.0)then
+!  if(n_core_orb.ne.0)then
-!   print*,'core orbitals'
+!   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
-!   pause
+!                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
-   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
+!  else
-                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+!   call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
-  else
+!                               , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
-   call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
+!  endif
-                               , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
-  endif
-!   call non_hrmt_bieig( mo_num, dm_tmp&
-!                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
-!                      , mo_num, natorb_tc_eigval )
 
+  call non_hrmt_bieig(mo_num, dm_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag & 
+                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo                  & 
+                      , mo_num, natorb_tc_eigval )
   accu = 0.d0
   do i = 1, mo_num
     print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)

src/tc_bi_ortho/test_s2_tc.irp.f

@@ -14,12 +14,14 @@ program test_tc
   read_wf = .True.
   touch read_wf
 
-  call routine_test_s2
+  call provide_all_three_ints_bi_ortho()
-  call routine_test_s2_davidson
+  call routine_h_triple_left
+  call routine_h_triple_right
+!  call routine_test_s2_davidson
 
 end
 
-subroutine routine_test_s2
+subroutine routine_h_triple_right
  implicit none
  logical           :: do_right 
  integer :: sze ,i, N_st, j
@@ -29,67 +31,65 @@ subroutine routine_test_s2
  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 Left '
+ print*,'Checking first the Right '
- do_right = .False.
- do i = 1, sze
-  u_0(i,1) = psi_l_coef_bi_ortho(i,1)
- enddo
- call H_tc_u_0_nstates_openmp(v_0_ref,u_0,N_st,sze, do_right)
- s_0_ref = 0.d0
- do i = 1, sze
-  do j = 1, sze
-    call get_s2(psi_det(1,1,i),psi_det(1,1,j),N_int,sij)
-    s_0_ref(i,1) += u_0(j,1) * sij
-  enddo
- enddo
- call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,u_0,N_st,sze, do_right)
- accu_e = 0.d0
- accu_s = 0.d0
- accu_e_0 = 0.d0
- accu_s_0 = 0.d0
- do i = 1, sze
-  accu_e_0 += v_0_ref(i,1) * psi_r_coef_bi_ortho(i,1)
-  accu_s_0 += s_0_ref(i,1) * psi_r_coef_bi_ortho(i,1)
-  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
- print*,'accu_e_0 = ',accu_e_0
- print*,'accu_s_0 = ',accu_s_0
-
- print*,'Checking then the right '
- do_right = .True.
  do i = 1, sze
   u_0(i,1) = psi_r_coef_bi_ortho(i,1)
  enddo
- call H_tc_u_0_nstates_openmp(v_0_ref,u_0,N_st,sze, do_right)
+ double precision :: wall0,wall1
- s_0_ref = 0.d0
+ call wall_time(wall0)
- do i = 1, sze
+ call H_tc_s2_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
-  do j = 1, sze
+ call wall_time(wall1)
-    call get_s2(psi_det(1,1,i),psi_det(1,1,j),N_int,sij)
+ print*,'time for omp',wall1 - wall0
-    s_0_ref(i,1) += u_0(j,1) * sij
+ call wall_time(wall0)
-  enddo
+ call H_tc_s2_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
- enddo
+ call wall_time(wall1)
- call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,u_0,N_st,sze, do_right)
+ print*,'time serial ',wall1 - wall0
  accu_e = 0.d0
  accu_s = 0.d0
- accu_e_0 = 0.d0
- accu_s_0 = 0.d0
  do i = 1, sze
-  accu_e_0 += v_0_ref(i,1) * psi_l_coef_bi_ortho(i,1)
-  accu_s_0 += s_0_ref(i,1) * psi_l_coef_bi_ortho(i,1)
   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
- print*,'accu_e_0 = ',accu_e_0
- print*,'accu_s_0 = ',accu_s_0
-
 
 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(:)

src/tc_bi_ortho/test_tc_two_rdm.irp.f

@@ -25,49 +25,43 @@ end
 
 subroutine  test
  implicit none
- integer :: h1,p1,h2,p2,i,j,istate
+ integer :: h1,p1,h2,p2,i,j,istate,s1,s2
- double precision :: rdm, integral, accu,ref
+ double precision :: rdm, integral, accu,ref, accu_new ,rdm_new
  double precision :: hmono, htwoe, hthree, htot
  accu = 0.d0
+ accu_new = 0.d0
  do h1 = 1, mo_num
   do p1 = 1, mo_num
    do h2 = 1, mo_num
     do p2 = 1, mo_num
      integral = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
-     rdm = tc_two_rdm(p1,h1,p2,h2)
+     rdm = tc_two_rdm(p2,p1,h2,h1)
-!     if(dabs(rdm).gt.1.d-10)then
-!      print*,h1,p1,h2,p2
-!      print*,rdm,integral,rdm*integral
-!     endif
      accu += integral * rdm
+     rdm_new = 0.d0
+     do s2 = 1, 2
+      do s1 = 1, 2
+       rdm_new += tc_two_rdm_s1s2(p2,p1,h2,h1,s1,s2)
+      enddo
+     enddo
+     accu_new += integral * rdm_new
     enddo
    enddo
   enddo
  enddo
  accu *= 0.5d0
- print*,'accu = ',accu
+ accu_new *= 0.5d0
-! print*,mo_bi_ortho_tc_two_e(2,15,2,1)
+ print*,'accu     = ',accu
-! print*,mo_bi_ortho_tc_two_e(15,2,2,1)
+ print*,'accu_new = ',accu_new
-! print*,mo_bi_ortho_tc_two_e(2,1,2,15)
-! print*,mo_bi_ortho_tc_two_e(2,1,15,2)
  ref = 0.d0
  do i = 1, N_det
   do j = 1, N_det
-!   if(i.eq.j)cycle
    call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
    do istate = 1,N_states
-!  print*,'i,j',i,j
-!  print*,psi_l_coef_bi_ortho(i,istate) , psi_r_coef_bi_ortho(j,istate) , htwoe 
-!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
-!   if(i.ne.j)then
-!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) ,  htwoe 
-!  print*,psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * htwoe 
-!   endif
     ref += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * state_average_weight(istate) * htwoe 
    enddo
   enddo
  enddo
- print*,' ref = ',ref 
+ print*,' ref     = ',ref 
  print*,'delta= ',ref-accu
 
 end

src/tc_bi_ortho/two_rdm_naive.irp.f

@@ -1,4 +1,5 @@
-BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
+ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist_s1s2, (mo_num, mo_num, mo_num, mo_num, 2,2)]
  implicit none
  BEGIN_DOC
  ! tc_two_rdm_chemist(p,s,q,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = CHEMIST NOTATION
@@ -14,6 +15,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
  other_spin(2) = 1
  allocate(occ(N_int*bit_kind_size,2))
  tc_two_rdm_chemist = 0.d0
+ tc_two_rdm_chemist_s1s2 = 0.d0
 
  do i = 1, N_det ! psi_left 
   do j = 1, N_det ! psi_right 
@@ -21,14 +23,16 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
    if(degree.gt.2)cycle
    if(degree.gt.0)then
     ! get excitation operators: from psi_det(j) --> psi_det(i)
-    call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
+    ! T_{j-->i} = a^p1_s1 a_h1_s1 
-    call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+     call get_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,degree,phase,N_int)
-    contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
+     call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-    do istate = 2, N_states
+     contrib = psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * phase * state_average_weight(1)
-     contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
+     do istate = 2, N_states
-    enddo
+      contrib += psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,istate) * phase * state_average_weight(istate)
+     enddo
     if(degree == 2)then
      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
     else if(degree==1)then
      ! occupation of the determinant psi_det(j)
      call bitstring_to_list_ab(psi_det(1,1,j), occ, n_occ_ab, N_int) 
@@ -40,6 +44,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
       h2 = m 
       p2 = m 
       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
      enddo
      ! run over the electrons of same spin than the excitation
      s2 = s1
@@ -48,6 +53,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
       h2 = m 
       p2 = m 
       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
      enddo
     endif
    else if(degree == 0)then
@@ -69,6 +75,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
        h2 = m 
        p2 = m 
        call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
       enddo
       ! run over the couple of alpha-alpha electrons 
       s2 = s1
@@ -78,6 +85,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
        p2 = m 
        if(h2.le.h1)cycle
        call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
       enddo
     enddo
     s1 = 2
@@ -92,6 +100,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num,
        p2 = m 
        if(h2.le.h1)cycle
        call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist_s1s2(1,1,1,1,s1,s2) ,mo_num,contrib)
       enddo
     enddo
    endif
@@ -124,12 +133,13 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
 end
 
 
-BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+&BEGIN_PROVIDER [ double precision, tc_two_rdm_s1s2, (mo_num, mo_num, mo_num, mo_num,2,2)]
  implicit none
  BEGIN_DOC
  ! tc_two_rdm(p,q,s,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = PHYSICIST NOTATION
  END_DOC
- integer :: p,q,r,s
+ integer :: p,q,r,s,s1,s2
  do r = 1, mo_num
   do q = 1, mo_num
    do s = 1, mo_num
@@ -139,5 +149,18 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
    enddo
   enddo
  enddo
+ do s2 = 1, 2
+  do s1 = 1, 2
+    do r = 1, mo_num
+     do q = 1, mo_num
+      do s = 1, mo_num
+       do p = 1, mo_num
+        tc_two_rdm_s1s2(p,q,s,r,s1,s2) = tc_two_rdm_chemist_s1s2(p,s,q,r,s1,s2) 
+       enddo
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
 
 END_PROVIDER