LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
Anthony Scemama 2025-02-12 12:00:28 +01:00
commit e38ecc0a49
8 changed files with 366 additions and 249 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
src/casscf_cipsi/gradient.irp.f
View File

@ -14,8 +14,8 @@ END_PROVIDER
implicit none implicit none
n_c_a_prov = n_core_inact_orb * n_act_orb n_c_a_prov = n_core_inact_orb * n_act_orb
n_c_v_prov = n_core_inact_orb * n_virt_orb n_c_v_prov = n_core_inact_orb * n_virt_orb
n_a_v_prov = n_act_orb * n_virt_orb n_a_v_prov = n_act_orb * n_virt_orb
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] BEGIN_PROVIDER [integer, excit, (2,nMonoEx)]
&BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] &BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)]
@ -28,7 +28,7 @@ END_PROVIDER
BEGIN_DOC BEGIN_DOC
! a list of the orbitals involved in the excitation ! a list of the orbitals involved in the excitation
END_DOC END_DOC
implicit none implicit none
integer :: i,t,a,ii,tt,aa,indx,indx_tmp integer :: i,t,a,ii,tt,aa,indx,indx_tmp
indx=0 indx=0
@ -48,7 +48,7 @@ END_PROVIDER
mat_idx_c_a(ii,tt) = indx mat_idx_c_a(ii,tt) = indx
end do end do
end do end do
indx_tmp = 0 indx_tmp = 0
do ii=1,n_core_inact_orb do ii=1,n_core_inact_orb
i=list_core_inact(ii) i=list_core_inact(ii)
@ -61,11 +61,11 @@ END_PROVIDER
indx_tmp += 1 indx_tmp += 1
list_idx_c_v(1,indx_tmp) = indx list_idx_c_v(1,indx_tmp) = indx
list_idx_c_v(2,indx_tmp) = ii list_idx_c_v(2,indx_tmp) = ii
list_idx_c_v(3,indx_tmp) = aa list_idx_c_v(3,indx_tmp) = aa
mat_idx_c_v(ii,aa) = indx mat_idx_c_v(ii,aa) = indx
end do end do
end do end do
indx_tmp = 0 indx_tmp = 0
do tt=1,n_act_orb do tt=1,n_act_orb
t=list_act(tt) t=list_act(tt)
@ -82,7 +82,7 @@ END_PROVIDER
mat_idx_a_v(tt,aa) = indx mat_idx_a_v(tt,aa) = indx
end do end do
end do end do
if (bavard) then if (bavard) then
write(6,*) ' Filled the table of the Monoexcitations ' write(6,*) ' Filled the table of the Monoexcitations '
do indx=1,nMonoEx do indx=1,nMonoEx
@ -90,7 +90,7 @@ END_PROVIDER
,excit(2,indx),' ',excit_class(indx) ,excit(2,indx),' ',excit_class(indx)
end do end do
end if end if
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
@ -104,7 +104,7 @@ END_PROVIDER
implicit none implicit none
integer :: i,t,a,indx integer :: i,t,a,indx
real*8 :: gradvec_it,gradvec_ia,gradvec_ta real*8 :: gradvec_it,gradvec_ia,gradvec_ta
indx=0 indx=0
norm_grad_vec2_tab = 0.d0 norm_grad_vec2_tab = 0.d0
do i=1,n_core_inact_orb do i=1,n_core_inact_orb
@ -114,7 +114,7 @@ END_PROVIDER
norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx) norm_grad_vec2_tab(1) += gradvec2(indx)*gradvec2(indx)
end do end do
end do end do
do i=1,n_core_inact_orb do i=1,n_core_inact_orb
do a=1,n_virt_orb do a=1,n_virt_orb
indx+=1 indx+=1
@ -122,7 +122,7 @@ END_PROVIDER
norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx) norm_grad_vec2_tab(2) += gradvec2(indx)*gradvec2(indx)
end do end do
end do end do
do t=1,n_act_orb do t=1,n_act_orb
do a=1,n_virt_orb do a=1,n_virt_orb
indx+=1 indx+=1
@ -130,7 +130,7 @@ END_PROVIDER
norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx) norm_grad_vec2_tab(3) += gradvec2(indx)*gradvec2(indx)
end do end do
end do end do
norm_grad_vec2=0.d0 norm_grad_vec2=0.d0
do indx=1,nMonoEx do indx=1,nMonoEx
norm_grad_vec2+=gradvec2(indx)*gradvec2(indx) norm_grad_vec2+=gradvec2(indx)*gradvec2(indx)
@ -144,7 +144,7 @@ END_PROVIDER
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad_vec2 write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad_vec2
write(6,*) write(6,*)
endif endif
END_PROVIDER END_PROVIDER
real*8 function gradvec_it(i,t) real*8 function gradvec_it(i,t)
@ -154,23 +154,30 @@ real*8 function gradvec_it(i,t)
END_DOC END_DOC
implicit none implicit none
integer :: i,t integer :: i,t
integer :: ii,tt,v,vv,x,y integer :: ii,tt,v,vv,x,y
integer :: x3,y3 integer :: x3,y3
double precision :: bielec_PQxx_no double precision :: bielec_PQxx_no
ii=list_core_inact(i) ii=list_core_inact(i)
tt=list_act(t) tt=list_act(t)
gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
gradvec_it-=occnum(tt)*Fipq(ii,tt) gradvec_it-=occnum(tt)*Fipq(ii,tt)
do v=1,n_act_orb ! active do y=1,n_act_orb ! active
vv=list_act(v) ! y3=y+n_core_inact_orb ! list_act(y)
do x=1,n_act_orb ! active do x=1,n_act_orb ! active
x3=x+n_core_inact_orb ! list_act(x) ! x3=x+n_core_inact_orb ! list_act(x)
do y=1,n_act_orb ! active do v=1,n_act_orb ! active
y3=y+n_core_inact_orb ! list_act(y) vv=list_act(v)
! Gamma(2) a a a a 1/r12 i a a a ! Gamma(2) a a a a 1/r12 i a a a
gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) ! gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
integer :: ichol
double precision :: tmp
tmp = 0.d0
do ichol=1,cholesky_mo_num
tmp = tmp + cholesky_no_total_transp(ichol,vv,ii) * cholesky_no_total_transp(ichol,list_act(x),list_act(y))
enddo
gradvec_it = gradvec_it - 2.D0*P0tuvx_no(t,v,x,y)*tmp
end do end do
end do end do
end do end do
@ -183,12 +190,12 @@ real*8 function gradvec_ia(i,a)
END_DOC END_DOC
implicit none implicit none
integer :: i,a,ii,aa integer :: i,a,ii,aa
ii=list_core_inact(i) ii=list_core_inact(i)
aa=list_virt(a) aa=list_virt(a)
gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii))
gradvec_ia*=2.D0 gradvec_ia*=2.D0
end function gradvec_ia end function gradvec_ia
real*8 function gradvec_ta(t,a) real*8 function gradvec_ta(t,a)
@ -198,7 +205,7 @@ real*8 function gradvec_ta(t,a)
END_DOC END_DOC
implicit none implicit none
integer :: t,a,tt,aa,v,vv,x,y integer :: t,a,tt,aa,v,vv,x,y
tt=list_act(t) tt=list_act(t)
aa=list_virt(a) aa=list_virt(a)
gradvec_ta=0.D0 gradvec_ta=0.D0
@ -211,6 +218,6 @@ real*8 function gradvec_ta(t,a)
end do end do
end do end do
gradvec_ta*=2.D0 gradvec_ta*=2.D0
end function gradvec_ta end function gradvec_ta

248
src/casscf_cipsi/hessian.irp.f
View File

@ -11,13 +11,14 @@ real*8 function hessmat_itju(i,t,j,u)
integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj
real*8 :: term,t2 real*8 :: term,t2
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
ii=list_core_inact(i) ii=list_core_inact(i)
tt=list_act(t) tt=list_act(t)
if (i.eq.j) then if (i.eq.j) then
if (t.eq.u) then if (t.eq.u) then
! diagonal element ! diagonal element
term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & term = occnum(tt)*Fipq(ii,ii) + &
2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) &
-2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) -2.D0*(Fipq(ii,ii)+Fapq(ii,ii))
term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i)) term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i))
term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) & term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) &
@ -83,10 +84,10 @@ real*8 function hessmat_itju(i,t,j,u)
end do end do
end do end do
end if end if
term*=2.D0 term*=2.D0
hessmat_itju=term hessmat_itju=term
end function hessmat_itju end function hessmat_itju
real*8 function hessmat_itja(i,t,j,a) real*8 function hessmat_itja(i,t,j,a)
@ -97,7 +98,7 @@ real*8 function hessmat_itja(i,t,j,a)
integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y
real*8 :: term real*8 :: term
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
! it/ja ! it/ja
ii=list_core_inact(i) ii=list_core_inact(i)
tt=list_act(t) tt=list_act(t)
@ -120,7 +121,7 @@ real*8 function hessmat_itja(i,t,j,a)
end if end if
term*=2.D0 term*=2.D0
hessmat_itja=term hessmat_itja=term
end function hessmat_itja end function hessmat_itja
real*8 function hessmat_itua(i,t,u,a) real*8 function hessmat_itua(i,t,u,a)
@ -131,7 +132,7 @@ real*8 function hessmat_itua(i,t,u,a)
integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3
real*8 :: term real*8 :: term
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
ii=list_core_inact(i) ii=list_core_inact(i)
tt=list_act(t) tt=list_act(t)
t3=t+n_core_inact_orb t3=t+n_core_inact_orb
@ -162,7 +163,7 @@ real*8 function hessmat_itua(i,t,u,a)
end if end if
term*=2.D0 term*=2.D0
hessmat_itua=term hessmat_itua=term
end function hessmat_itua end function hessmat_itua
real*8 function hessmat_iajb(i,a,j,b) real*8 function hessmat_iajb(i,a,j,b)
@ -173,7 +174,7 @@ real*8 function hessmat_iajb(i,a,j,b)
integer :: i,a,j,b,ii,aa,jj,bb integer :: i,a,j,b,ii,aa,jj,bb
real*8 :: term real*8 :: term
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
ii=list_core_inact(i) ii=list_core_inact(i)
aa=list_virt(a) aa=list_virt(a)
if (i.eq.j) then if (i.eq.j) then
@ -199,7 +200,7 @@ real*8 function hessmat_iajb(i,a,j,b)
end if end if
term*=2.D0 term*=2.D0
hessmat_iajb=term hessmat_iajb=term
end function hessmat_iajb end function hessmat_iajb
real*8 function hessmat_iatb(i,a,t,b) real*8 function hessmat_iatb(i,a,t,b)
@ -210,7 +211,7 @@ real*8 function hessmat_iatb(i,a,t,b)
integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3
real*8 :: term real*8 :: term
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
ii=list_core_inact(i) ii=list_core_inact(i)
aa=list_virt(a) aa=list_virt(a)
tt=list_act(t) tt=list_act(t)
@ -231,7 +232,7 @@ real*8 function hessmat_iatb(i,a,t,b)
end if end if
term*=2.D0 term*=2.D0
hessmat_iatb=term hessmat_iatb=term
end function hessmat_iatb end function hessmat_iatb
real*8 function hessmat_taub(t,a,u,b) real*8 function hessmat_taub(t,a,u,b)
@ -240,83 +241,186 @@ real*8 function hessmat_taub(t,a,u,b)
END_DOC END_DOC
implicit none implicit none
integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y
integer :: v3,x3 integer :: v3,x3, ichol
real*8 :: term,t1,t2,t3 real*8 :: term,t1,t2,t3, tmp
double precision :: bielec_pqxx_no,bielec_pxxq_no double precision :: bielec_pqxx_no,bielec_pxxq_no
double precision, allocatable :: tmp1(:), tmp2(:,:)
allocate(tmp1(n_act_orb))
allocate(tmp2(n_act_orb,n_act_orb))
tt=list_act(t) tt=list_act(t)
aa=list_virt(a) aa=list_virt(a)
if (t == u) then if (t == u) then
if (a == b) then if (a == b) then
! ta/ta ! ta/ta
t1=occnum(tt)*Fipq(aa,aa) t1=occnum(tt)*Fipq(aa,aa) - occnum(tt)*Fipq(tt,tt)
t2=0.D0 t2=0.D0
t3=0.D0 ! do x=1,n_act_orb
t1-=occnum(tt)*Fipq(tt,tt) ! x3=x+n_core_inact_orb
! do v=1,n_act_orb
! v3=v+n_core_inact_orb
! tmp = 0.d0
! do ichol = 1, cholesky_mo_num
! tmp = tmp + cholesky_no_total_transp(ichol,aa,aa) * cholesky_no_total_transp(ichol,v3,x3)
! enddo
! t2 = t2 + 2.D0*P0tuvx_no(t,t,v,x)*tmp
! enddo
! enddo
do x=1,n_act_orb
x3=x+n_core_inact_orb
call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, &
cholesky_no_total_transp(1,aa,aa), 1, 0.d0, &
tmp1, 1)
do v=1,n_act_orb
t2 = t2 + P0tuvx_no(t,t,v,x)*tmp1(v)
enddo
enddo
! do v=1,n_act_orb
! v3=v+n_core_inact_orb
! do x=1,n_act_orb
! x3=x+n_core_inact_orb
! tmp = 0.d0
! do ichol = 1, cholesky_mo_num
! tmp = tmp + cholesky_no_total_transp(ichol,aa,x3) * cholesky_no_total_transp(ichol,v3,aa)
! enddo
! t2 = t2 + 2.d0*(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))*tmp
! end do
! end do
call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, &
cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, 0.d0, &
tmp2, n_act_orb)
do v=1,n_act_orb do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb do x=1,n_act_orb
xx=list_act(x) t2 = t2 + P0tuvx_no(t,x,v,t)*tmp2(x,v) + P0tuvx_no(t,x,t,v)*tmp2(x,v)
x3=x+n_core_inact_orb enddo
t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & enddo
+(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* &
bielec_pxxq_no(aa,x3,v3,aa)) t3=0.D0
do y=1,n_act_orb do x=1,n_act_orb
t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) xx=list_act(x)
do y=1,n_act_orb
do v=1,n_act_orb
t3 = t3 - P0tuvx_no(t,v,x,y)*bielecCI_no(v,t,y,xx)
end do end do
end do end do
end do end do
term=t1+t2+t3 term=t1+t2+t3*2.d0
else else
bb=list_virt(b) bb=list_virt(b)
! ta/tb b/=a ! ta/tb b/=a
term=occnum(tt)*Fipq(aa,bb) term=occnum(tt)*Fipq(aa,bb)
! 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)
! end do
! end do
do x=1,n_act_orb
x3=x+n_core_inact_orb
call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, &
cholesky_no_total_transp(1,aa,bb), 1, 0.d0, &
tmp1, 1)
do v=1,n_act_orb
term = term + P0tuvx_no(t,t,v,x)*tmp1(v)
enddo
enddo
! 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,x,v,t)+P0tuvx_no(t,x,t,v))*bielec_pxxq_no(aa,x3,v3,bb)
! end do
! end do
call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, &
cholesky_no_total_transp(1,n_core_inact_orb+1,bb), cholesky_mo_num, 0.d0, &
tmp2, n_act_orb)
do v=1,n_act_orb do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb do x=1,n_act_orb
xx=list_act(x) term = term + P0tuvx_no(t,x,v,t)*tmp2(x,v) + P0tuvx_no(t,x,t,v)*tmp2(x,v)
x3=x+n_core_inact_orb enddo
term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & enddo
+(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) &
*bielec_pxxq_no(aa,x3,v3,bb))
end do
end do
end if end if
else else
! ta/ub t/=u ! ta/ub t/=u
uu=list_act(u) uu=list_act(u)
bb=list_virt(b) bb=list_virt(b)
term=0.D0 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)
! end do
! end do
do x=1,n_act_orb
x3=x+n_core_inact_orb
call dgemv('T', cholesky_mo_num, n_act_orb, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,x3), cholesky_mo_num, &
cholesky_no_total_transp(1,aa,bb), 1, 0.d0, &
tmp1, 1)
do v=1,n_act_orb
term = term + P0tuvx_no(t,u,v,x)*tmp1(v)
enddo
enddo
! 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,x,v,u)+P0tuvx_no(t,x,u,v))*bielec_pxxq_no(aa,x3,v3,bb)
! end do
! end do
call dgemm('T','N', n_act_orb, n_act_orb, cholesky_mo_num, 2.d0, &
cholesky_no_total_transp(1,n_core_inact_orb+1,aa), cholesky_mo_num, &
cholesky_no_total_transp(1,n_core_inact_orb+1,bb), cholesky_mo_num, 0.d0, &
tmp2, n_act_orb)
do v=1,n_act_orb do v=1,n_act_orb
vv=list_act(v)
v3=v+n_core_inact_orb
do x=1,n_act_orb do x=1,n_act_orb
xx=list_act(x) term = term + P0tuvx_no(t,x,v,u)*tmp2(x,v)+P0tuvx_no(t,x,u,v)*tmp2(x,v)
x3=x+n_core_inact_orb enddo
term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & enddo
+(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) &
*bielec_pxxq_no(aa,x3,v3,bb))
end do
end do
if (a.eq.b) then if (a.eq.b) then
term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu))
do v=1,n_act_orb do v=1,n_act_orb
do y=1,n_act_orb do y=1,n_act_orb
do x=1,n_act_orb do x=1,n_act_orb
term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) term = 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) - P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt)
end do end do
end do end do
end do end do
end if end if
end if end if
term*=2.D0 term*=2.D0
hessmat_taub=term hessmat_taub=term
end function hessmat_taub end function hessmat_taub
BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
@ -326,7 +430,7 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
implicit none implicit none
integer :: i,t,a,indx,indx_shift integer :: i,t,a,indx,indx_shift
real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & !$OMP SHARED(hessdiag,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) &
!$OMP PRIVATE(i,indx,t,a,indx_shift) !$OMP PRIVATE(i,indx,t,a,indx_shift)
@ -339,9 +443,9 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
end do end do
end do end do
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
indx_shift = n_core_inact_orb*n_act_orb indx_shift = n_core_inact_orb*n_act_orb
!$OMP DO !$OMP DO
do a=1,n_virt_orb do a=1,n_virt_orb
do i=1,n_core_inact_orb do i=1,n_core_inact_orb
indx = a + (i-1)*n_virt_orb + indx_shift indx = a + (i-1)*n_virt_orb + indx_shift
@ -349,9 +453,9 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
end do end do
end do end do
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
indx_shift += n_core_inact_orb*n_virt_orb indx_shift += n_core_inact_orb*n_virt_orb
!$OMP DO !$OMP DO
do a=1,n_virt_orb do a=1,n_virt_orb
do t=1,n_act_orb do t=1,n_act_orb
indx = a + (t-1)*n_virt_orb + indx_shift indx = a + (t-1)*n_virt_orb + indx_shift
@ -360,7 +464,7 @@ BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)]
end do end do
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
END_PROVIDER END_PROVIDER
@ -377,7 +481,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
real*8 :: hessmat_taub real*8 :: hessmat_taub
! c-a c-v a-v ! c-a c-v a-v
! c-a | X X X ! c-a | X X X
! c-v | X X ! c-v | X X
! a-v | X ! a-v | X
provide all_mo_integrals provide all_mo_integrals
@ -390,12 +494,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP DO !$OMP DO
!!!! < Core-active| H |Core-active > !!!! < Core-active| H |Core-active >
! Core-active excitations ! Core-active excitations
do indx_tmp = 1, n_c_a_prov do indx_tmp = 1, n_c_a_prov
indx = list_idx_c_a(1,indx_tmp) indx = list_idx_c_a(1,indx_tmp)
i = list_idx_c_a(2,indx_tmp) i = list_idx_c_a(2,indx_tmp)
t = list_idx_c_a(3,indx_tmp) t = list_idx_c_a(3,indx_tmp)
! Core-active excitations ! Core-active excitations
do j = 1, n_core_inact_orb do j = 1, n_core_inact_orb
if (i.eq.j) then if (i.eq.j) then
ustart=t ustart=t
@ -418,12 +522,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP DO !$OMP DO
!!!! < Core-active| H |Core-VIRTUAL > !!!! < Core-active| H |Core-VIRTUAL >
! Core-active excitations ! Core-active excitations
do indx_tmp = 1, n_c_a_prov do indx_tmp = 1, n_c_a_prov
indx = list_idx_c_a(1,indx_tmp) indx = list_idx_c_a(1,indx_tmp)
i = list_idx_c_a(2,indx_tmp) i = list_idx_c_a(2,indx_tmp)
t = list_idx_c_a(3,indx_tmp) t = list_idx_c_a(3,indx_tmp)
! Core-VIRTUAL excitations ! Core-VIRTUAL excitations
do jndx_tmp = 1, n_c_v_prov do jndx_tmp = 1, n_c_v_prov
jndx = list_idx_c_v(1,jndx_tmp) jndx = list_idx_c_v(1,jndx_tmp)
j = list_idx_c_v(2,jndx_tmp) j = list_idx_c_v(2,jndx_tmp)
@ -441,12 +545,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP DO !$OMP DO
!!!! < Core-active| H |ACTIVE-VIRTUAL > !!!! < Core-active| H |ACTIVE-VIRTUAL >
! Core-active excitations ! Core-active excitations
do indx_tmp = 1, n_c_a_prov do indx_tmp = 1, n_c_a_prov
indx = list_idx_c_a(1,indx_tmp) indx = list_idx_c_a(1,indx_tmp)
i = list_idx_c_a(2,indx_tmp) i = list_idx_c_a(2,indx_tmp)
t = list_idx_c_a(3,indx_tmp) t = list_idx_c_a(3,indx_tmp)
! ACTIVE-VIRTUAL excitations ! ACTIVE-VIRTUAL excitations
do jndx_tmp = 1, n_a_v_prov do jndx_tmp = 1, n_a_v_prov
jndx = list_idx_a_v(1,jndx_tmp) jndx = list_idx_a_v(1,jndx_tmp)
u = list_idx_a_v(2,jndx_tmp) u = list_idx_a_v(2,jndx_tmp)
@ -466,12 +570,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx) !$OMP PRIVATE(indx_tmp,indx,i,a,j,b,bstart,jndx)
!$OMP DO !$OMP DO
!!!!! < Core-VIRTUAL | H |Core-VIRTUAL > !!!!! < Core-VIRTUAL | H |Core-VIRTUAL >
! Core-VIRTUAL excitations ! Core-VIRTUAL excitations
do indx_tmp = 1, n_c_v_prov do indx_tmp = 1, n_c_v_prov
indx = list_idx_c_v(1,indx_tmp) indx = list_idx_c_v(1,indx_tmp)
i = list_idx_c_v(2,indx_tmp) i = list_idx_c_v(2,indx_tmp)
a = list_idx_c_v(3,indx_tmp) a = list_idx_c_v(3,indx_tmp)
! Core-VIRTUAL excitations ! Core-VIRTUAL excitations
do j = 1, n_core_inact_orb do j = 1, n_core_inact_orb
if (i.eq.j) then if (i.eq.j) then
bstart=a bstart=a
@ -485,7 +589,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
enddo enddo
enddo enddo
enddo enddo
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
!$OMP END PARALLEL !$OMP END PARALLEL
endif endif
@ -496,12 +600,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP DO !$OMP DO
!!!! < Core-VIRTUAL | H |Active-VIRTUAL > !!!! < Core-VIRTUAL | H |Active-VIRTUAL >
! Core-VIRTUAL excitations ! Core-VIRTUAL excitations
do indx_tmp = 1, n_c_v_prov do indx_tmp = 1, n_c_v_prov
indx = list_idx_c_v(1,indx_tmp) indx = list_idx_c_v(1,indx_tmp)
i = list_idx_c_v(2,indx_tmp) i = list_idx_c_v(2,indx_tmp)
a = list_idx_c_v(3,indx_tmp) a = list_idx_c_v(3,indx_tmp)
! Active-VIRTUAL excitations ! Active-VIRTUAL excitations
do jndx_tmp = 1, n_a_v_prov do jndx_tmp = 1, n_a_v_prov
jndx = list_idx_a_v(1,jndx_tmp) jndx = list_idx_a_v(1,jndx_tmp)
t = list_idx_a_v(2,jndx_tmp) t = list_idx_a_v(2,jndx_tmp)
@ -520,12 +624,12 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP DO !$OMP DO
!!!! < Active-VIRTUAL | H |Active-VIRTUAL > !!!! < Active-VIRTUAL | H |Active-VIRTUAL >
! Active-VIRTUAL excitations ! Active-VIRTUAL excitations
do indx_tmp = 1, n_a_v_prov do indx_tmp = 1, n_a_v_prov
indx = list_idx_a_v(1,indx_tmp) indx = list_idx_a_v(1,indx_tmp)
t = list_idx_a_v(2,indx_tmp) t = list_idx_a_v(2,indx_tmp)
a = list_idx_a_v(3,indx_tmp) a = list_idx_a_v(3,indx_tmp)
! Active-VIRTUAL excitations ! Active-VIRTUAL excitations
do u=t,n_act_orb do u=t,n_act_orb
if (t.eq.u) then if (t.eq.u) then
bstart=a bstart=a
@ -542,4 +646,4 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)]
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
!$OMP END PARALLEL !$OMP END PARALLEL
END_PROVIDER END_PROVIDER

117
src/casscf_cipsi/mcscf_fock.irp.f
View File

@ -3,37 +3,53 @@ BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
! the inactive Fock matrix, in molecular orbitals ! the inactive Fock matrix, in molecular orbitals
END_DOC END_DOC
implicit none implicit none
integer :: p,q,k,kk,t,tt,u,uu integer :: i,p,q,k,kk,t,tt,u,uu
double precision :: bielec_pxxq_no, bielec_pqxx_no
do q=1,mo_num do q=1,mo_num
do p=1,mo_num do p=1,mo_num
Fipq(p,q)=one_ints_no(p,q) Fipq(p,q)=one_ints_no(p,q)
end do end do
end do end do
! the inactive Fock matrix ! the inactive Fock matrix
do k=1,n_core_inact_orb do k=1,n_core_inact_orb
kk=list_core_inact(k) kk=list_core_inact(k)
do q=1,mo_num ! do q=1,mo_num
do p=1,mo_num ! do p=1,mo_num
Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q) ! do i=1,cholesky_mo_num
end do ! Fipq(p,q) = Fipq(p,q) + 2.d0* cholesky_no_total_transp(i,p,q) * cholesky_no_total_transp(i,kk,kk)
end do ! enddo
! end do
! end do
call dgemv('T', cholesky_mo_num, mo_num*mo_num, 2.d0, &
cholesky_no_total_transp, cholesky_mo_num, &
cholesky_no_total_transp(1,kk,kk), 1, 1.d0, &
Fipq, 1)
! do q=1,mo_num
! do p=1,mo_num
! do i=1,cholesky_mo_num
! Fipq(p,q) = Fipq(p,q) - cholesky_no_total_transp(i,p,kk) * cholesky_no_total_transp(i,kk,q)
! enddo
! end do
! end do
call dgemm('T','N', mo_num, mo_num, cholesky_mo_num, -1.d0, &
cholesky_no_total_transp(1,1,kk), cholesky_mo_num, &
cholesky_no_total_transp(1,kk,1), cholesky_mo_num*mo_num, 1.d0, &
Fipq, mo_num)
end do end do
if (bavard) then if (bavard) then
integer :: i
write(6,*) write(6,*)
write(6,*) ' the diagonal of the inactive effective Fock matrix ' write(6,*) ' the diagonal of the inactive effective Fock matrix '
write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
write(6,*) write(6,*)
end if end if
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
BEGIN_DOC BEGIN_DOC
! the active active Fock matrix, in molecular orbitals ! the active active Fock matrix, in molecular orbitals
@ -45,27 +61,42 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
END_DOC END_DOC
implicit none implicit none
integer :: p,q,k,kk,t,tt,u,uu integer :: p,q,k,kk,t,tt,u,uu
double precision :: bielec_pxxq_no, bielec_pqxx_no
Fapq = 0.d0 Fapq = 0.d0
! the active Fock matrix, D0tu is diagonal ! the active Fock matrix, D0tu is diagonal
do t=1,n_act_orb do t=1,n_act_orb
tt=list_act(t) tt=list_act(t)
do q=1,mo_num ! do q=1,mo_num
do p=1,mo_num ! do p=1,mo_num
Fapq(p,q)+=occnum(tt) & ! do i=1,cholesky_mo_num
*(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q)) ! Fapq(p,q) = Fapq(p,q) + occnum(tt)* cholesky_no_total_transp(i,p,q) * cholesky_no_total_transp(i,tt,tt)
end do ! enddo
end do ! end do
! end do
call dgemv('T', cholesky_mo_num, mo_num*mo_num, occnum(tt), &
cholesky_no_total_transp, cholesky_mo_num, &
cholesky_no_total_transp(1,tt,tt), 1, 1.d0, &
Fapq, 1)
! do q=1,mo_num
! do p=1,mo_num
! do i=1,cholesky_mo_num
! Fapq(p,q) = Fapq(p,q) - 0.5d0*occnum(tt)*cholesky_no_total_transp(i,p,tt) * cholesky_no_total_transp(i,tt,q)
! enddo
! end do
! end do
call dgemm('T','N', mo_num, mo_num, cholesky_mo_num, -0.5d0*occnum(tt), &
cholesky_no_total_transp(1,1,tt), cholesky_mo_num, &
cholesky_no_total_transp(1,tt,1), cholesky_mo_num*mo_num, 1.d0, &
Fapq, mo_num)
end do end do
if (bavard) then if (bavard) then
integer :: i integer :: i
write(6,*) write(6,*)
write(6,*) ' the effective Fock matrix over MOs' write(6,*) ' the effective Fock matrix over MOs'
write(6,*) write(6,*)
write(6,*) write(6,*)
write(6,*) ' the diagonal of the inactive effective Fock matrix ' write(6,*) ' the diagonal of the inactive effective Fock matrix '
write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
@ -75,35 +106,35 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num)
write(6,*) write(6,*)
end if end if
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)] BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)]
&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)] &BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis ! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis
END_DOC END_DOC
SCF_density_matrix_ao_alpha = D0tu_alpha_ao SCF_density_matrix_ao_alpha = D0tu_alpha_ao
SCF_density_matrix_ao_beta = D0tu_beta_ao SCF_density_matrix_ao_beta = D0tu_beta_ao
soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta
mcscf_fock_beta_ao = fock_matrix_ao_beta mcscf_fock_beta_ao = fock_matrix_ao_beta
mcscf_fock_alpha_ao = fock_matrix_ao_alpha mcscf_fock_alpha_ao = fock_matrix_ao_alpha
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)] BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)] &BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis ! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis
END_DOC END_DOC
call ao_to_mo(mcscf_fock_alpha_ao,ao_num,mcscf_fock_alpha_mo,mo_num) call ao_to_mo(mcscf_fock_alpha_ao,ao_num,mcscf_fock_alpha_mo,mo_num)
call ao_to_mo(mcscf_fock_beta_ao,ao_num,mcscf_fock_beta_mo,mo_num) call ao_to_mo(mcscf_fock_beta_ao,ao_num,mcscf_fock_beta_mo,mo_num)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_mo, (mo_num,mo_num) ] BEGIN_PROVIDER [ double precision, mcscf_fock_mo, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, mcscf_fock_diag_mo, (mo_num)] &BEGIN_PROVIDER [ double precision, mcscf_fock_diag_mo, (mo_num)]
@ -118,13 +149,13 @@ END_PROVIDER
! |-----------------------| ! |-----------------------|
! | Fcv | F^a | Rvv | ! | Fcv | F^a | Rvv |
! !
! C: Core, O: Open, V: Virtual ! C: Core, O: Open, V: Virtual
! !
! Rcc = Acc Fcc^a + Bcc Fcc^b ! Rcc = Acc Fcc^a + Bcc Fcc^b
! Roo = Aoo Foo^a + Boo Foo^b ! Roo = Aoo Foo^a + Boo Foo^b
! Rvv = Avv Fvv^a + Bvv Fvv^b ! Rvv = Avv Fvv^a + Bvv Fvv^b
! Fcv = (F^a + F^b)/2 ! Fcv = (F^a + F^b)/2
! !
! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO) ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO)
! A,B: Coupling parameters ! A,B: Coupling parameters
! !
@ -133,7 +164,7 @@ END_PROVIDER
! cc oo vv ! cc oo vv
! A -0.5 0.5 1.5 ! A -0.5 0.5 1.5
! B 1.5 0.5 -0.5 ! B 1.5 0.5 -0.5
! !
END_DOC END_DOC
integer :: i,j,n integer :: i,j,n
if (elec_alpha_num == elec_beta_num) then if (elec_alpha_num == elec_beta_num) then
@ -194,4 +225,4 @@ END_PROVIDER
do i = 1, mo_num do i = 1, mo_num
mcscf_fock_diag_mo(i) = mcscf_fock_mo(i,i) mcscf_fock_diag_mo(i) = mcscf_fock_mo(i,i)
enddo enddo
END_PROVIDER END_PROVIDER

102
src/casscf_cipsi/natorb.irp.f
View File

@ -72,84 +72,27 @@ BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
BEGIN_DOC BEGIN_DOC
! 4-index transformation of 2part matrices ! 4-index transformation of 2part matrices
END_DOC END_DOC
integer :: i,j,k,l,p,q
real*8 :: d(n_act_orb)
! index per index double precision, allocatable :: tmp(:,:,:,:)
! first quarter allocate(tmp(n_act_orb,n_act_orb,n_act_orb,n_act_orb))
P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:)
do j=1,n_act_orb call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
do k=1,n_act_orb P0tuvx, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
do l=1,n_act_orb tmp, (n_act_orb*n_act_orb*n_act_orb))
do p=1,n_act_orb
d(p)=0.D0 call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
end do tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
do p=1,n_act_orb P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb))
do q=1,n_act_orb
d(p)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p) call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
end do P0tuvx_no, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
end do tmp, (n_act_orb*n_act_orb*n_act_orb))
do p=1,n_act_orb
P0tuvx_no(p,j,k,l)=d(p) call dgemm('T','N',(n_act_orb*n_act_orb*n_act_orb), n_act_orb, n_act_orb, 1.d0, &
end do tmp, n_act_orb, natorbsCI, n_act_orb, 0.d0, &
end do P0tuvx_no, (n_act_orb*n_act_orb*n_act_orb))
end do
end do deallocate(tmp)
! 2nd quarter
do j=1,n_act_orb
do k=1,n_act_orb
do l=1,n_act_orb
do p=1,n_act_orb
d(p)=0.D0
end do
do p=1,n_act_orb
do q=1,n_act_orb
d(p)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
P0tuvx_no(j,p,k,l)=d(p)
end do
end do
end do
end do
! 3rd quarter
do j=1,n_act_orb
do k=1,n_act_orb
do l=1,n_act_orb
do p=1,n_act_orb
d(p)=0.D0
end do
do p=1,n_act_orb
do q=1,n_act_orb
d(p)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
P0tuvx_no(j,k,p,l)=d(p)
end do
end do
end do
end do
! 4th quarter
do j=1,n_act_orb
do k=1,n_act_orb
do l=1,n_act_orb
do p=1,n_act_orb
d(p)=0.D0
end do
do p=1,n_act_orb
do q=1,n_act_orb
d(p)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
end do
end do
do p=1,n_act_orb
P0tuvx_no(j,k,l,p)=d(p)
end do
end do
end do
end do
END_PROVIDER END_PROVIDER
@ -160,6 +103,7 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
BEGIN_DOC BEGIN_DOC
! Transformed one-e integrals ! Transformed one-e integrals
END_DOC END_DOC
integer :: i,j, p, q integer :: i,j, p, q
real*8 :: d(n_act_orb) real*8 :: d(n_act_orb)
one_ints_no(:,:)=mo_one_e_integrals(:,:) one_ints_no(:,:)=mo_one_e_integrals(:,:)
@ -168,10 +112,8 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
do j=1,mo_num do j=1,mo_num
do p=1,n_act_orb do p=1,n_act_orb
d(p)=0.D0 d(p)=0.D0
end do
do p=1,n_act_orb
do q=1,n_act_orb do q=1,n_act_orb
d(p)+=one_ints_no(list_act(q),j)*natorbsCI(q,p) d(p) = d(p) + one_ints_no(list_act(q),j)*natorbsCI(q,p)
end do end do
end do end do
do p=1,n_act_orb do p=1,n_act_orb
@ -183,8 +125,6 @@ BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
do j=1,mo_num do j=1,mo_num
do p=1,n_act_orb do p=1,n_act_orb
d(p)=0.D0 d(p)=0.D0
end do
do p=1,n_act_orb
do q=1,n_act_orb do q=1,n_act_orb
d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p) d(p)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
end do end do

41
src/dav_general_mat/dav_general.irp.f
View File

@ -82,7 +82,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
nproc_target = nproc nproc_target = nproc
double precision :: rss double precision :: rss
integer :: maxab integer :: maxab
maxab = sze maxab = sze
m=1 m=1
disk_based = .False. disk_based = .False.
@ -204,7 +204,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
u_in(i,k) = r1*dcos(r2) u_in(i,k) = r1*dcos(r2)
enddo enddo
enddo enddo
! Normalize all states ! Normalize all states
do k=1,N_st_diag do k=1,N_st_diag
call normalize(u_in(:,k),sze) call normalize(u_in(:,k),sze)
enddo enddo
@ -228,20 +228,13 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
shift = N_st_diag*(iter-1) shift = N_st_diag*(iter-1)
shift2 = N_st_diag*iter shift2 = N_st_diag*iter
if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i><i|H|u_k>
! Compute |W_k> = \sum_i |i><i|H|u_k> ! -----------------------------------
! -----------------------------------
! Gram-Smitt to orthogonalize all new guess with the previous vectors ! Gram-Smitt to orthogonalize all new guess with the previous vectors
call ortho_qr(U,size(U,1),sze,shift2) call ortho_qr(U,size(U,1),sze,shift2)
call ortho_qr(U,size(U,1),sze,shift2)
! call H_S2_u_0_nstates_openmp(W(:,shift+1),U(:,shift+1),N_st_diag,sze) call hpsi(W(:,shift+1),U(:,shift+1),N_st_diag,sze,h_mat)
call hpsi(W(:,shift+1),U(:,shift+1),N_st_diag,sze,h_mat)
else
! Already computed in update below
continue
endif
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l> ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
! ------------------------------------------- ! -------------------------------------------
@ -311,12 +304,12 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
do i=1,sze do i=1,sze
U(i,shift2+k) = & U(i,shift2+k) = &
(lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
/max(H_jj(i) - lambda (k),1.d-2) /max(dabs(H_jj(i) - lambda (k)),1.d-2) * dsign(1d0,H_jj(i) - lambda (k))
enddo enddo
if (k <= N_st) then if (k <= N_st) then
residual_norm(k) = u_dot_u(U(:,shift2+k),sze) residual_norm(k) = u_dot_u(U(:,shift2+k),sze)
to_print(1,k) = lambda(k) to_print(1,k) = lambda(k)
to_print(2,k) = residual_norm(k) to_print(2,k) = residual_norm(k)
endif endif
enddo enddo
@ -324,7 +317,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
if ((itertot>1).and.(iter == 1)) then if ((itertot>1).and.(iter == 1)) then
!don't print !don't print
continue continue
else else
write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,ES11.3))') iter-1, to_print(1:2,1:N_st) write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,ES11.3))') iter-1, to_print(1:2,1:N_st)
@ -333,11 +326,11 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
! Check convergence ! Check convergence
if (iter > 1) then if (iter > 1) then
converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
endif endif
do k=1,N_st do k=1,N_st
if (residual_norm(k) > 1.e8) then if (residual_norm(k) > 1.d8) then
print *, 'Davidson failed' print *, 'Davidson failed'
stop -1 stop -1
endif endif
@ -365,13 +358,15 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, &
U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
do k=1,N_st_diag do k=1,N_st_diag
do i=1,sze do i=1,sze
U(i,k) = u_in(i,k) U(i,k) = u_in(i,k)
enddo enddo
enddo enddo
call ortho_qr(U,size(U,1),sze,N_st_diag)
call ortho_qr(U,size(U,1),sze,N_st_diag) call ortho_qr(U,size(U,1),sze,N_st_diag)
do j=1,N_st_diag do j=1,N_st_diag
k=1 k=1
do while ((k<sze).and.(U(k,j) == 0.d0)) do while ((k<sze).and.(U(k,j) == 0.d0))
@ -412,7 +407,7 @@ subroutine hpsi(v,u,N_st,sze,h_mat)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Computes $v = H | u \rangle$ and ! Computes $v = H | u \rangle$ and
END_DOC END_DOC
integer, intent(in) :: N_st,sze integer, intent(in) :: N_st,sze
double precision, intent(in) :: u(sze,N_st),h_mat(sze,sze) double precision, intent(in) :: u(sze,N_st),h_mat(sze,sze)

6
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -185,9 +185,9 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_extra, (ao_num,n_points_extra_final_grid)] BEGIN_PROVIDER[double precision, aos_in_r_array_extra, (ao_num,n_points_extra_final_grid)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point ! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid
END_DOC END_DOC
integer :: i,j integer :: i,j
double precision :: aos_array(ao_num), r(3) double precision :: aos_array(ao_num), r(3)
@ -214,7 +214,7 @@ END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)] BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)]
BEGIN_DOC BEGIN_DOC
! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point of the EXTRA grid
END_DOC END_DOC
implicit none implicit none

41
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -181,3 +181,44 @@
END_PROVIDER END_PROVIDER
!!!!!EXTRA GRID
BEGIN_PROVIDER[double precision, mos_in_r_array_extra_omp, (mo_num,n_points_extra_final_grid)]
implicit none
BEGIN_DOC
! mos_in_r_array_extra(i,j) = value of the ith mo on the jth grid point on the EXTRA GRID
END_DOC
integer :: i,j
double precision :: mos_array_extra(mo_num), r(3)
print*,'coucou'
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,mos_array_extra,j) &
!$OMP SHARED(mos_in_r_array_extra_omp,n_points_extra_final_grid,mo_num,final_grid_points_extra)
do i = 1, n_points_extra_final_grid
r(1) = final_grid_points_extra(1,i)
r(2) = final_grid_points_extra(2,i)
r(3) = final_grid_points_extra(3,i)
call give_all_mos_at_r(r,mos_array_extra)
do j = 1, mo_num
mos_in_r_array_extra_omp(j,i) = mos_array_extra(j)
enddo
enddo
!$OMP END PARALLEL DO
print*,'coucou fin'
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_in_r_array_extra_transp,(n_points_extra_final_grid,mo_num)]
implicit none
BEGIN_DOC
! mos_in_r_array_extra_transp(i,j) = value of the jth mo on the ith grid point
END_DOC
integer :: i,j
do i = 1, n_points_extra_final_grid
do j = 1, mo_num
mos_in_r_array_extra_transp(i,j) = mos_in_r_array_extra_omp(j,i)
enddo
enddo
END_PROVIDER

1
src/mo_two_e_ints/map_integrals.irp.f
View File

@ -86,7 +86,6 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_s
call set_multiple_levels_omp(.False.) call set_multiple_levels_omp(.False.)
!$OMP PARALLEL DO PRIVATE(k,l,ii) SCHEDULE(dynamic) !$OMP PARALLEL DO PRIVATE(k,l,ii) SCHEDULE(dynamic)
do l=mo_integrals_cache_min,mo_integrals_cache_max do l=mo_integrals_cache_min,mo_integrals_cache_max
do k=mo_integrals_cache_min,mo_integrals_cache_max do k=mo_integrals_cache_min,mo_integrals_cache_max