From b2fa6b0b9c0da217d87ff7fc348fef87a76a0e3b Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 21 Aug 2023 11:17:58 +0200 Subject: [PATCH] inactive --> active gradient are OK for real ! --- src/casscf_tc_bi/grad_dm.irp.f | 27 +++++++++++++++++++-- src/tc_bi_ortho/two_rdm_naive.irp.f | 37 ++++++++++++++++++++++------- 2 files changed, 53 insertions(+), 11 deletions(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index d7d7046d..f62acbdd 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/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 diff --git a/src/tc_bi_ortho/two_rdm_naive.irp.f b/src/tc_bi_ortho/two_rdm_naive.irp.f index 8fd34975..3963d09e 100644 --- a/src/tc_bi_ortho/two_rdm_naive.irp.f +++ b/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) = = 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) = = 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