all the one-body gradients are correct for TC-CASSCF

src/casscf_tc_bi/grad_dm.irp.f

@@ -1,11 +1,29 @@
  BEGIN_PROVIDER [real*8, gradvec_tc_r, (0:3,nMonoEx)]
 &BEGIN_PROVIDER [real*8, gradvec_tc_l,  (0:3,nMonoEx)]
+ BEGIN_DOC
+! gradvec_tc_r(0:3,i) = <Chi_0| H E_q^p(i) |Phi_0> 
+! 
+! gradvec_tc_l(0:3,i) = <Chi_0| E_p^q(i) H |Phi_0> 
+!
+! where the indices "i" corresponds to E_q^p(i) 
+! 
+! i = mat_idx_c_a(q,p) 
+!
+! and gradvec_tc_r/l(0) = full matrix element
+!
+!     gradvec_tc_r/l(1) = one-body part 
+
+!     gradvec_tc_r/l(2) = two-body part 
+
+!     gradvec_tc_r/l(3) = three-body part 
+ END_DOC
  implicit none
  integer                        :: ii,tt,aa,indx
  integer                        :: i,t,a,fff
  double precision :: res_l(0:3), res_r(0:3)
  gradvec_tc_l = 0.d0
  gradvec_tc_r = 0.d0
+ ! computing the core/inactive --> virtual orbitals gradients
   do i=1,n_core_inact_orb
     ii=list_core_inact(i)
     do t=1,n_act_orb
@@ -33,9 +51,15 @@
   end do
   
   do t=1,n_act_orb
+    tt=list_act(t)
     do a=1,n_virt_orb
-      indx = mat_idx_a_v(i,a) 
-!      gradvec_tc_l(indx)=gradvec_ta(t,a)
+      aa=list_virt(a)
+      indx = mat_idx_a_v(t,a) 
+      call gradvec_tc_ta(tt,aa,res_l, res_r)
+      do fff = 0,3
+       gradvec_tc_l(fff,indx)=res_l(fff)
+       gradvec_tc_r(fff,indx)=res_r(fff)
+      enddo
     end do
   end do
 END_PROVIDER 
@@ -65,7 +89,7 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
 !
 ! Corresponds to res_r = <X0|H E_i^t|Phi_0>
 !
-!                res_l = <X0|E_i^t H |Phi_0>
+!                res_l = <X0|E_t^i H |Phi_0>
  END_DOC
  integer, intent(in) :: i,t
  double precision, intent(out) :: res_l(0:3),res_r(0:3)
@@ -83,3 +107,32 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
  enddo
  
 end
+
+subroutine gradvec_tc_ta(t,a,res_l, res_r)
+ implicit none
+ BEGIN_DOC
+! active --> virtual TC gradient 
+!
+! Corresponds to res_r = <X0|H E_t^a|Phi_0>
+!
+!                res_l = <X0|E_a^t H |Phi_0>
+ END_DOC
+ integer, intent(in) :: t,a
+ double precision, intent(out) :: res_l(0:3),res_r(0:3)
+ integer :: rr,r,m
+ double precision :: dm
+ res_r = 0.d0
+ res_l = 0.d0
+! do rr = 1, n_act_orb
+!  r = list_act(rr) 
+!  res_l(1) +=  mo_bi_ortho_tc_one_e(a,r) * tc_transition_matrix_mo(t,r,1,1)
+!  res_r(1) += -mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1)
+! enddo
+ do m = 1, mo_num
+  res_r(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(a,m,1,1) &
+             -mo_bi_ortho_tc_one_e(m,a) * tc_transition_matrix_mo(m,t,1,1)
+  res_l(1) += mo_bi_ortho_tc_one_e(a,m) * tc_transition_matrix_mo(t,m,1,1) &
+             -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,a,1,1)
+ enddo
+ 
+end

src/casscf_tc_bi/grad_old.irp.f

@@ -38,15 +38,19 @@
    enddo
   enddo
  enddo
-!  do indx=1,nMonoEx
-!    ihole=excit(1,indx)
-!    ipart=excit(2,indx)
-!    call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
-!    do ll = 0, 3
-!     gradvec_detail_left_old (ll,indx)=res_l(ll)
-!     gradvec_detail_right_old(ll,indx)=res_r(ll)
-!    enddo
-!  end do
+
+ do tt = 1, n_act_orb
+  ihole = list_act(tt)
+  do aa = 1, n_virt_orb
+   ipart = list_virt(aa)
+   indx = mat_idx_a_v(tt,aa) 
+   call calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
+   do ll = 0, 3
+    gradvec_detail_left_old (ll,indx)=res_l(ll)
+    gradvec_detail_right_old(ll,indx)=res_r(ll)
+   enddo
+  enddo
+ enddo
   
   real*8                         :: norm_grad_left, norm_grad_right
   norm_grad_left=0.d0