From b3b080929b38c152ac1bd868bad4dd311108e7e9 Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Tue, 11 Jul 2023 01:56:28 +0200
Subject: [PATCH] fixed stupid bug in TC 1-RDM and one-e gradient: o-v, o-a are
 ok

---
 src/casscf_tc_bi/grad_dm.irp.f  | 28 ++++++++++++++++------------
 src/casscf_tc_bi/grad_old.irp.f | 21 +++++++++++++++++----
 src/tc_bi_ortho/tc_prop.irp.f   |  8 +++++---
 3 files changed, 38 insertions(+), 19 deletions(-)

diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f
index 0fc2e4eb..7f6155ab 100644
--- a/src/casscf_tc_bi/grad_dm.irp.f
+++ b/src/casscf_tc_bi/grad_dm.irp.f
@@ -11,8 +11,7 @@
     do t=1,n_act_orb
       tt=list_act(t)
       indx = mat_idx_c_a(i,t) 
-      call gradvec_tc_it(ii,tt,res_l)
-      call gradvec_tc_it(tt,ii,res_r)
+      call gradvec_tc_it(ii,tt,res_l,res_r)
       do fff = 0,3
        gradvec_tc_l(fff,indx)=res_l(fff)
        gradvec_tc_r(fff,indx)=res_r(fff)
@@ -56,23 +55,28 @@ subroutine gradvec_tc_ia(i,a,res)
  
 end
 
-subroutine gradvec_tc_it(i,t,res)
+subroutine gradvec_tc_it(i,t,res_l, res_r)
  implicit none
  BEGIN_DOC
 ! doubly occupied --> active TC gradient 
 !
-! Corresponds to <X0|H E_i^t|Phi_0>
+! Corresponds to res_r = <X0|H E_i^t|Phi_0>
+!
+!                res_l = <X0|E_i^t H |Phi_0>
  END_DOC
  integer, intent(in) :: i,t
- double precision, intent(out) :: res(0:3)
- integer :: rr,r,ss,s
+ double precision, intent(out) :: res_l(0:3),res_r(0:3)
+ integer :: rr,r,ss,s,m
  double precision :: dm
- res = 0.d0
- res(1) = -2 * mo_bi_ortho_tc_one_e(i,t)
- do rr = 1, n_act_orb
-  r = list_act(rr)
-  dm = tc_transition_matrix_mo(t,r,1,1)
-  res(1) += mo_bi_ortho_tc_one_e(i,r) * dm
+ res_r = 0.d0
+ do m = 1, mo_num
+  res_r(1) += mo_bi_ortho_tc_one_e(i,m) * tc_transition_matrix_mo(t,m,1,1) & 
+             -mo_bi_ortho_tc_one_e(m,t) * tc_transition_matrix_mo(m,i,1,1)
+ enddo
+ res_l = 0.d0
+ do m = 1, mo_num
+  res_l(1) += mo_bi_ortho_tc_one_e(t,m) * tc_transition_matrix_mo(i,m,1,1) & 
+             -mo_bi_ortho_tc_one_e(m,i) * tc_transition_matrix_mo(m,t,1,1)
  enddo
  
 end
diff --git a/src/casscf_tc_bi/grad_old.irp.f b/src/casscf_tc_bi/grad_old.irp.f
index 6610dee3..ea6747b1 100644
--- a/src/casscf_tc_bi/grad_old.irp.f
+++ b/src/casscf_tc_bi/grad_old.irp.f
@@ -25,6 +25,19 @@
    enddo
   enddo
  enddo
+
+ do ii = 1, n_core_inact_orb
+  ihole = list_core_inact(ii)
+  do tt = 1, n_act_orb
+   ipart = list_act(tt)
+   indx = mat_idx_c_a(ii,tt) 
+   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
 !  do indx=1,nMonoEx
 !    ihole=excit(1,indx)
 !    ipart=excit(2,indx)
@@ -57,7 +70,7 @@ END_PROVIDER
 subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
   BEGIN_DOC
   ! eq 18 of Siegbahn et al, Physica Scripta 1980
-  ! we calculate res_l  = <Phi| H^tc E_pq | Psi>, and res_r = <Phi| E_qp H^tc | Psi>
+  ! we calculate res_r  = <Phi| H^tc E_pq | Psi>, and res_r = <Phi| E_qp H^tc | Psi>
   ! q=hole, p=particle
   ! res_l(0) =   total matrix element
   ! res_l(1) =   one-electron part 
@@ -89,12 +102,12 @@ subroutine calc_grad_elem_h_tc(ihole,ipart,res_l, res_r)
       if (ierr.eq.1) then
 
         call i_H_tc_psi_phi(det_mu_ex,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)
+            ,N_det,psi_det_size,N_states,i_H_chi_array,i_H_phi_array)
 !        print*,i_H_chi_array(1,1),i_H_phi_array(1,1)
         do istate=1,N_states
          do ll = 0,3
-          res_l(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
-          res_r(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase
+          res_l(ll)+=i_H_phi_array(ll,istate)*psi_l_coef_bi_ortho(mu,istate)*phase
+          res_r(ll)+=i_H_chi_array(ll,istate)*psi_r_coef_bi_ortho(mu,istate)*phase
          enddo
         end do
       end if
diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f
index 5bb0e2c0..a13dc9a2 100644
--- a/src/tc_bi_ortho/tc_prop.irp.f
+++ b/src/tc_bi_ortho/tc_prop.irp.f
@@ -29,7 +29,7 @@
          tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         enddo
         do p = 1, n_occ_ab(2) ! browsing the beta electrons
-         m = occ(p,1)
+         m = occ(p,2)
          tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         enddo
        else
@@ -38,12 +38,14 @@
           ! Single alpha
           h = exc(1,1,1) ! hole in psi_det(1,1,j) 
           p = exc(1,2,1) ! particle in psi_det(1,1,j) 
-          tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+          tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= &
+          phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         else
           ! Single beta
           h = exc(1,1,2) ! hole in psi_det(1,1,j) 
           p = exc(1,2,2) ! particle in psi_det(1,1,j) 
-          tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
+          tc_transition_matrix_mo_beta(p,h,istate,jstate)+=  &
+          phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
         endif
        endif
      enddo