LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

all gradients are OK !

Browse Source
This commit is contained in:
eginer 2023-08-22 13:48:01 +02:00
parent 99d856570a
commit 50f2fb0bfa
2 changed files with 102 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
src/casscf_tc_bi/grad_dm.irp.f
View File

@ -50,11 +50,13 @@
end do
end do
! print*,'DM grad'
do t=1,n_act_orb
tt=list_act(t)
do a=1,n_virt_orb
aa=list_virt(a)
indx = mat_idx_a_v(t,a)
! print*,indx,t,a
call gradvec_tc_ta(tt,aa,res_l, res_r)
do fff = 0,3
gradvec_tc_l(fff,indx)=res_l(fff)
@ -83,7 +85,7 @@ subroutine gradvec_tc_ia(i,a,res_l, res_r)
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))
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)
@ -93,12 +95,14 @@ 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 * ( &
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)) &
+tc_transition_matrix_mo(r,t,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
res_r(0) = res_r(1) + res_r(2) + res_r(3)
res_l(0) = res_l(1) + res_l(2) + res_l(3)
end
subroutine gradvec_tc_it(i,t,res_l, res_r)
@ -116,20 +120,22 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
res_r = 0.d0
res_l = 0.d0
res_r(1) += -2.d0 * mo_bi_ortho_tc_one_e(i,t)
res_l(1) += 2.D0 * mo_bi_ortho_tc_one_e(t,i)
res_l(1) -= -2.D0 * mo_bi_ortho_tc_one_e(t,i)
do rr = 1, n_act_orb
r = list_act(rr)
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)
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))
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))
res_l(2) -= -2.d0 * (2d0 * mo_bi_ortho_tc_two_e(t,j,i,j) - mo_bi_ortho_tc_two_e(t,j,j,i))
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))
res_l(2) -= tc_transition_matrix_mo(r,t,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(r,j,i,j) - mo_bi_ortho_tc_two_e(j,r,j,i))
enddo
enddo
do rr = 1, n_act_orb
@ -140,14 +146,21 @@ subroutine gradvec_tc_it(i,t,res_l, res_r)
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)) &
)
res_l(2) -= -0.5d0 * ( &
tc_transition_matrix_mo(r,u,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(r,t,u,i) - mo_bi_ortho_tc_two_e(t,r,u,i)) &
+ tc_transition_matrix_mo(u,r,1,1) * (2.d0 * mo_bi_ortho_tc_two_e(t,u,i,r) - mo_bi_ortho_tc_two_e(u,t,i,r)) &
)
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) )
res_l(2) -= 0.5d0 * ( &
mo_bi_ortho_tc_two_e(v,u,i,r) * tc_two_rdm(v,u,t,r) + mo_bi_ortho_tc_two_e(v,u,r,i) * tc_two_rdm(v,u,r,t) )
enddo
enddo
enddo
res_r(0) = res_r(1) + res_r(2) + res_r(3)
res_l(0) = res_l(1) + res_l(2) + res_l(3)
end
subroutine gradvec_tc_ta(t,a,res_l, res_r)
@ -161,14 +174,90 @@ subroutine gradvec_tc_ta(t,a,res_l, res_r)
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
integer :: rr,r,j,jj,u,uu,v,vv
double precision :: res_r_inact_test, res_r_act_test
double precision :: res_l_inact_test, res_l_act_test
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)
res_r(1) -= mo_bi_ortho_tc_one_e(r,a) * tc_transition_matrix_mo(r,t,1,1)
enddo
res_r_inact_test = 0.d0
res_l_inact_test = 0.d0
do jj = 1, n_core_inact_orb
j = list_core_inact(jj)
do rr = 1, n_act_orb
r = list_act(rr)
res_r_inact_test += -tc_transition_matrix_mo(r,t,1,1) * &
(2.d0 * mo_bi_ortho_tc_two_e(r,j,a,j) - mo_bi_ortho_tc_two_e(r,j,j,a))
res_l_inact_test += -tc_transition_matrix_mo(t,r,1,1) * &
(2.d0 * mo_bi_ortho_tc_two_e(a,j,r,j) - mo_bi_ortho_tc_two_e(j,a,r,j))
enddo
enddo
res_r_act_test = 0.d0
res_l_act_test = 0.d0
do rr = 1, n_act_orb
r = list_act(rr)
do vv = 1, n_act_orb
v = list_act(vv)
do uu = 1, n_act_orb
u = list_act(uu)
res_r_act_test += - (mo_bi_ortho_tc_two_e(v,r,u,a) * tc_two_rdm(r,v,t,u) &
+mo_bi_ortho_tc_two_e(v,r,a,u) * tc_two_rdm(r,v,u,t))
res_l_act_test += - (mo_bi_ortho_tc_two_e(u,a,v,r) * tc_two_rdm(t,u,r,v) &
+mo_bi_ortho_tc_two_e(a,u,v,r) * tc_two_rdm(u,t,r,v))
enddo
enddo
enddo
res_r_act_test *= 0.5d0
res_l_act_test *= 0.5d0
res_r(2) = res_r_inact_test + res_r_act_test
res_l(2) = res_l_inact_test + res_l_act_test
integer :: m,x,y
double precision :: res_r_inact, res_r_act
if(.False.)then
! test quantities
res_r_inact = 0.d0
res_r_act = 0.d0
do m = 1, mo_num
do x = 1, mo_num
do jj = 1, n_core_inact_orb
j = list_core_inact(jj)
res_r_inact += 0.5d0 * mo_bi_ortho_tc_two_e(t,j,m,x) * tc_two_rdm(a,j,m,x) &
-0.5d0 * mo_bi_ortho_tc_two_e(m,j,a,x) * tc_two_rdm(m,j,t,x) &
+0.5d0 * mo_bi_ortho_tc_two_e(j,t,m,x) * tc_two_rdm(j,a,m,x) &
-0.5d0 * mo_bi_ortho_tc_two_e(x,j,m,a) * tc_two_rdm(x,j,m,t)
enddo
do rr = 1, n_act_orb
r = list_act(rr)
res_r_act += 0.5d0 * mo_bi_ortho_tc_two_e(t,r,m,x) * tc_two_rdm(a,r,m,x) &
-0.5d0 * mo_bi_ortho_tc_two_e(m,r,a,x) * tc_two_rdm(m,r,t,x) &
+0.5d0 * mo_bi_ortho_tc_two_e(r,t,m,x) * tc_two_rdm(r,a,m,x) &
-0.5d0 * mo_bi_ortho_tc_two_e(x,r,m,a) * tc_two_rdm(x,r,m,t)
enddo
enddo
enddo
if(dabs(res_r_inact).gt.1.d-12)then
if(dabs(res_r_inact_test - res_r_inact).gt.1.d-10)then
print*,'inact'
print*,'t,a',t,a
print*,res_r_inact_test , res_r_inact, dabs(res_r_inact_test - res_r_inact)
endif
endif
if(dabs(res_r_act).gt.1.d-12)then
if(dabs(res_r_act_test - res_r_act).gt.1.d-10)then
print*,'act'
print*,'t,a',t,a
print*,res_r_act_test , res_r_act, dabs(res_r_act_test - res_r_act)
endif
endif
endif
res_r(0) = res_r(1) + res_r(2) + res_r(3)
res_l(0) = res_l(1) + res_l(2) + res_l(3)
end

2
src/casscf_tc_bi/grad_old.irp.f
View File

@ -39,11 +39,13 @@
enddo
enddo
! print*,'old grad'
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)
! print*,indx,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)