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

Optimized Hessian

Browse Source
This commit is contained in:
Anthony Scemama 2019-07-06 02:17:07 +02:00
parent 970fd8837a
commit 3f69f95275
1 changed files with 52 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
src/casscf/hessian.irp.f
View File

@ -536,6 +536,9 @@ real*8 function hessmat_taub(t,a,u,b)
integer :: v3,x3
real*8 :: term,t1,t2,t3
double precision,allocatable :: P0tuvx_no_t(:,:,:)
double precision :: bielec_pqxx_no_2(n_act_orb,n_act_orb)
double precision :: bielec_pxxq_no_2(n_act_orb,n_act_orb)
tt=list_act(t)
aa=list_virt(a)
if (t == u) then
@ -545,59 +548,87 @@ real*8 function hessmat_taub(t,a,u,b)
t2=0.D0
t3=0.D0
t1-=occnum(tt)*Fipq(tt,tt)
do x=1,n_act_orb
xx=list_act(x)
x3=x+n_core_inact_orb
do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
t2+=P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3)
end do
end do
do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb
xx=list_act(x)
x3=x+n_core_inact_orb
t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) &
+(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* &
bielec_pxxq_no(aa,x3,v3,aa))
do y=1,n_act_orb
t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
t2+=(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* &
bielec_pxxq_no(aa,x3,v3,aa)
end do
end do
do y=1,n_act_orb
do x=1,n_act_orb
xx=list_act(x)
do v=1,n_act_orb
t3-=P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx)
end do
end do
end do
term=t1+t2+t3
term=t1+2.d0*(t2+t3)
else
bb=list_virt(b)
! ta/tb b/=a
term=occnum(tt)*Fipq(aa,bb)
term=0.5d0*occnum(tt)*Fipq(aa,bb)
do x=1,n_act_orb
xx=list_act(x)
x3=x+n_core_inact_orb
do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
term = term + P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3)
end do
end do
do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb
xx=list_act(x)
x3=x+n_core_inact_orb
term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) &
+(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) &
*bielec_pxxq_no(aa,x3,v3,bb))
term= term + (P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) &
*bielec_pxxq_no(aa,x3,v3,bb)
end do
end do
term += term
end if
else
! ta/ub t/=u
uu=list_act(u)
bb=list_virt(b)
term=0.D0
do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb
xx=list_act(x)
x3=x+n_core_inact_orb
term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) &
+(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) &
*bielec_pxxq_no(aa,x3,v3,bb))
allocate(P0tuvx_no_t(n_act_orb,n_act_orb,n_act_orb))
P0tuvx_no_t(:,:,:) = P0tuvx_no(t,:,:,:)
do x=1,n_act_orb
x3=x+n_core_inact_orb
do v=1,n_act_orb
v3=v+n_core_inact_orb
bielec_pqxx_no_2(v,x) = bielec_pqxx_no(aa,bb,v3,x3)
bielec_pxxq_no_2(v,x) = bielec_pxxq_no(aa,v3,x3,bb)
end do
end do
term=0.D0
do x=1,n_act_orb
do v=1,n_act_orb
term += P0tuvx_no_t(u,v,x)*bielec_pqxx_no_2(v,x)
term += bielec_pxxq_no_2(x,v) * (P0tuvx_no_t(x,v,u)+P0tuvx_no_t(x,u,v))
end do
end do
term = 6.d0*term
if (a.eq.b) then
term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu))
do v=1,n_act_orb
do y=1,n_act_orb
do x=1,n_act_orb
term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu)
term-=P0tuvx_no_t(v,x,y)*bielecCI_no(x,y,v,uu)
term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt)
end do
end do