inactive --> active gradient are OK for real !

src/casscf_tc_bi/grad_dm.irp.f

@@ -112,8 +112,7 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
  END_DOC
  integer, intent(in) :: i,t
  double precision, intent(out) :: res_l(0:3),res_r(0:3)
- integer :: rr,r,ss,s,m,mm
- double precision :: dm
+ integer :: rr,r,j,jj,u,uu,v,vv
  res_r = 0.d0
  res_l = 0.d0
  res_r(1)  += -2.d0 * mo_bi_ortho_tc_one_e(i,t) 
@@ -124,6 +123,30 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
   res_r(1) +=  mo_bi_ortho_tc_one_e(i,r) * tc_transition_matrix_mo(t,r,1,1)   
   res_l(1) += -mo_bi_ortho_tc_one_e(r,i) * tc_transition_matrix_mo(r,t,1,1)
  enddo
+
+ do jj = 1, n_core_inact_orb
+  j = list_core_inact(jj)
+  res_r(2) += 2.d0 * (2d0 * mo_bi_ortho_tc_two_e(i,j,t,j) - mo_bi_ortho_tc_two_e(j,i,t,j)) 
+  do rr = 1, n_act_orb
+   r = list_act(rr)
+   res_r(2) += tc_transition_matrix_mo(t,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(i,j,r,j) - mo_bi_ortho_tc_two_e(i,j,j,r)) 
+  enddo
+ enddo
+ do rr = 1, n_act_orb
+  r = list_act(rr)
+  do uu = 1, n_act_orb
+   u = list_act(uu)
+   res_r(2) += -0.5d0 * (  & 
+    tc_transition_matrix_mo(u,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(u,i,r,t) - mo_bi_ortho_tc_two_e(u,i,t,r))  &
+  + tc_transition_matrix_mo(r,u,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(i,r,t,u) - mo_bi_ortho_tc_two_e(i,r,u,t))  &
+    )
+   do vv = 1, n_act_orb
+    v = list_act(vv)
+    res_r(2) +=  0.5d0 * (  & 
+    mo_bi_ortho_tc_two_e(i,r,v,u) * tc_two_rdm(t,r,v,u) + mo_bi_ortho_tc_two_e(r,i,v,u) * tc_two_rdm(r,t,v,u) )
+   enddo
+  enddo
+ enddo
  
 end
 

src/tc_bi_ortho/two_rdm_naive.irp.f

@@ -1,7 +1,7 @@
-BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+BEGIN_PROVIDER [ double precision, tc_two_rdm_chemist, (mo_num, mo_num, mo_num, mo_num)]
  implicit none
  BEGIN_DOC
- ! tc_two_rdm(p,s,q,r) = <Phi| a^dager_p
+ ! tc_two_rdm_chemist(p,s,q,r) = <Phi| a^dagger_p a^dagger_q q_r a_s |Phi> = CHEMIST NOTATION
  END_DOC
  integer :: i,j,istate,m,mm,nn
  integer                        :: exc(0:2,2,2)
@@ -13,7 +13,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
  other_spin(1) = 2
  other_spin(2) = 1
  allocate(occ(N_int*bit_kind_size,2))
- tc_two_rdm = 0.d0
+ tc_two_rdm_chemist = 0.d0
 
  do i = 1, N_det ! psi_left 
   do j = 1, N_det ! psi_right 
@@ -28,7 +28,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
      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,mo_num,contrib)
+     call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,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) 
@@ -39,7 +39,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
       m = occ(mm,s2)
       h2 = m 
       p2 = m 
-      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
      enddo
      ! run over the electrons of same spin than the excitation
      s2 = s1
@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
       m = occ(mm,s2)
       h2 = m 
       p2 = m 
-      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+      call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
      enddo
     endif
    else if(degree == 0)then
@@ -68,7 +68,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        m = occ(mm,s2)
        h2 = m 
        p2 = m 
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
       ! run over the couple of alpha-alpha electrons 
       s2 = s1
@@ -77,7 +77,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        h2 = m 
        p2 = m 
        if(h2.le.h1)cycle
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
     enddo
     s1 = 2
@@ -91,7 +91,7 @@ BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
        h2 = m 
        p2 = m 
        if(h2.le.h1)cycle
-       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm,mo_num,contrib)
+       call update_tc_rdm(h1,p1,h2,p2,s1,s2,tc_two_rdm_chemist,mo_num,contrib)
       enddo
     enddo
    endif
@@ -122,3 +122,22 @@ subroutine update_tc_rdm(h1,p1,h2,p2,s1,s2,array,sze,contrib)
  endif
 
 end
+
+
+BEGIN_PROVIDER [ double precision, tc_two_rdm, (mo_num, mo_num, mo_num, mo_num)]
+ 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
+ do r = 1, mo_num
+  do q = 1, mo_num
+   do s = 1, mo_num
+    do p = 1, mo_num
+    tc_two_rdm(p,q,s,r) = tc_two_rdm_chemist(p,s,q,r) 
+    enddo
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER