From 6b56e213d89cd7ee07a33ec460ed7613fb9a28e6 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 17 Aug 2023 18:17:46 +0200 Subject: [PATCH 1/3] right two-body inactive-virtual gradients implemented --- src/casscf_tc_bi/grad_dm.irp.f | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 6e06f6ce..77336c93 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -79,7 +79,22 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r) res_r = 0.d0 res_l(1) = -2 * mo_bi_ortho_tc_one_e(a,i) res_r(1) = -2 * mo_bi_ortho_tc_one_e(i,a) - + integer :: j,t,r,jj,tt,rr + do jj = 1, n_core_inact_orb + j = list_core_inact(jj) + res_r(2) += -2.d0 * ( 2.d0 * mo_bi_ortho_tc_two_e(j,i,j,a) - mo_bi_ortho_tc_two_e(i,j,j,a)) + res_l(2) += -2.d0 * ( 2.d0 * mo_bi_ortho_tc_two_e(j,a,j,i) - mo_bi_ortho_tc_two_e(j,a,i,j)) + enddo + do tt = 1, n_act_orb + t = list_act(tt) + do rr = 1, n_act_orb + r = list_act(rr) + res_r(2) += -0.5d0 * ( & + tc_transition_matrix_mo(r,t,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,i,t,a) - mo_bi_ortho_tc_two_e(i,r,t,a)) & + +tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(t,i,r,a) - mo_bi_ortho_tc_two_e(i,t,r,a)) & + ) + enddo + enddo end subroutine gradvec_tc_it(i,t,res_l, res_r) From 30c238656f03e5c4b45fc31cb4f67dc3231c2d2d Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 21 Aug 2023 10:06:58 +0200 Subject: [PATCH 2/3] Gradient for inactive-->virtual work --- src/casscf_tc_bi/grad_dm.irp.f | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/casscf_tc_bi/grad_dm.irp.f b/src/casscf_tc_bi/grad_dm.irp.f index 77336c93..d7d7046d 100644 --- a/src/casscf_tc_bi/grad_dm.irp.f +++ b/src/casscf_tc_bi/grad_dm.irp.f @@ -93,6 +93,10 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r) tc_transition_matrix_mo(r,t,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,i,t,a) - mo_bi_ortho_tc_two_e(i,r,t,a)) & +tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(t,i,r,a) - mo_bi_ortho_tc_two_e(i,t,r,a)) & ) + res_l(2) += -0.5d0 * ( & + tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(t,a,r,i) - mo_bi_ortho_tc_two_e(t,a,i,r)) & + +tc_transition_matrix_mo(t,r,1,1) *(2.d0 * mo_bi_ortho_tc_two_e(r,a,t,i) - mo_bi_ortho_tc_two_e(r,a,i,t)) & + ) enddo enddo end From b2fa6b0b9c0da217d87ff7fc348fef87a76a0e3b Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 21 Aug 2023 11:17:58 +0200 Subject: [PATCH 3/3] 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