mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-09 06:53:38 +01:00
Compare commits
No commits in common. "6531181316c131f30d54a3653d17b597c0a43f3b" and "328ab2dadf856732af82aba1fd5386a3ab3ee909" have entirely different histories.
6531181316
...
328ab2dadf
@ -1,151 +1,104 @@
|
|||||||
BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
|
! -*- F90 -*-
|
||||||
BEGIN_DOC
|
BEGIN_PROVIDER[real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
|
||||||
! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
|
&BEGIN_PROVIDER[real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
|
||||||
! indices are unshifted orbital numbers
|
BEGIN_DOC
|
||||||
END_DOC
|
! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
|
||||||
implicit none
|
! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
|
||||||
integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
|
! indices are unshifted orbital numbers
|
||||||
double precision, allocatable :: integrals_array(:,:)
|
! all integrals are read from files
|
||||||
real*8 :: mo_two_e_integral
|
END_DOC
|
||||||
|
implicit none
|
||||||
allocate(integrals_array(mo_num,mo_num))
|
integer :: i,j,p,q,indx,kk
|
||||||
|
real*8 :: hhh
|
||||||
bielec_PQxx = 0.d0
|
logical :: lread
|
||||||
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
do j=i,n_core_orb
|
|
||||||
jj=list_core(j)
|
|
||||||
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
|
|
||||||
do p=1,mo_num
|
|
||||||
do q=1,mo_num
|
|
||||||
bielec_PQxx(p,q,i,j)=integrals_array(p,q)
|
|
||||||
bielec_PQxx(p,q,j,i)=integrals_array(p,q)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do j=1,n_act_orb
|
|
||||||
jj=list_act(j)
|
|
||||||
j3=j+n_core_orb
|
|
||||||
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
|
|
||||||
do p=1,mo_num
|
|
||||||
do q=1,mo_num
|
|
||||||
bielec_PQxx(p,q,i,j3)=integrals_array(p,q)
|
|
||||||
bielec_PQxx(p,q,j3,i)=integrals_array(p,q)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
|
do i=1,n_core_orb+n_act_orb
|
||||||
|
do j=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,mo_num
|
||||||
|
bielec_PQxx(p,q,i,j)=0.D0
|
||||||
|
bielec_PxxQ(p,i,j,q)=0.D0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
! (ij|pq)
|
open(unit=12,form='formatted',status='old',file='bielec_PQxx.tmp')
|
||||||
do i=1,n_act_orb
|
lread=.true.
|
||||||
ii=list_act(i)
|
indx=0
|
||||||
i3=i+n_core_orb
|
do while (lread)
|
||||||
do j=i,n_act_orb
|
read(12,*,iostat=kk) p,q,i,j,hhh
|
||||||
jj=list_act(j)
|
if (kk.ne.0) then
|
||||||
j3=j+n_core_orb
|
lread=.false.
|
||||||
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map)
|
else
|
||||||
do p=1,mo_num
|
! stored with p.le.q, i.le.j
|
||||||
do q=1,mo_num
|
bielec_PQxx(p,q,i,j)=hhh
|
||||||
bielec_PQxx(p,q,i3,j3)=integrals_array(p,q)
|
bielec_PQxx(q,p,i,j)=hhh
|
||||||
bielec_PQxx(p,q,j3,i3)=integrals_array(p,q)
|
bielec_PQxx(q,p,j,i)=hhh
|
||||||
end do
|
bielec_PQxx(p,q,j,i)=hhh
|
||||||
end do
|
indx+=1
|
||||||
end do
|
end if
|
||||||
end do
|
end do
|
||||||
|
close(12)
|
||||||
write(6,*) ' provided integrals (PQ|xx) '
|
write(6,*) ' read ',indx,' integrals PQxx into core '
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='old',file='bielec_PxxQ.tmp')
|
||||||
|
lread=.true.
|
||||||
|
indx=0
|
||||||
|
do while (lread)
|
||||||
|
read(12,*,iostat=kk) p,i,j,q,hhh
|
||||||
|
if (kk.ne.0) then
|
||||||
|
lread=.false.
|
||||||
|
else
|
||||||
|
! stored with (ip).le.(jq)
|
||||||
|
bielec_PxxQ(p,i,j,q)=hhh
|
||||||
|
bielec_PxxQ(q,j,i,p)=hhh
|
||||||
|
indx+=1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
write(6,*) ' read ',indx,' integrals PxxQ into core '
|
||||||
|
close(12)
|
||||||
|
write(6,*) ' provided integrals (PQ|xx) and (Px|xQ) '
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER[real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
|
||||||
|
! index p runs over the whole basis, t,u,v only over the active orbitals
|
||||||
|
! inegrals read from file
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,k,p,t,u,v,kk,indx
|
||||||
|
real*8 :: hhh
|
||||||
|
logical :: lread
|
||||||
|
|
||||||
|
write(6,*) ' reading integrals bielecCI '
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
|
do i=1,n_act_orb
|
||||||
BEGIN_DOC
|
do j=1,n_act_orb
|
||||||
! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
|
do k=1,n_act_orb
|
||||||
! indices are unshifted orbital numbers
|
do p=1,mo_num
|
||||||
END_DOC
|
bielecCI(i,k,j,p)=0.D0
|
||||||
implicit none
|
end do
|
||||||
integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
|
end do
|
||||||
double precision, allocatable :: integrals_array(:,:)
|
end do
|
||||||
real*8 :: mo_two_e_integral
|
end do
|
||||||
|
|
||||||
allocate(integrals_array(mo_num,mo_num))
|
|
||||||
|
|
||||||
bielec_PxxQ = 0.d0
|
|
||||||
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
do j=i,n_core_orb
|
|
||||||
jj=list_core(j)
|
|
||||||
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map)
|
|
||||||
do p=1,mo_num
|
|
||||||
do q=1,mo_num
|
|
||||||
bielec_PxxQ(p,i,j,q)=integrals_array(p,q)
|
|
||||||
bielec_PxxQ(p,j,i,q)=integrals_array(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do j=1,n_act_orb
|
|
||||||
jj=list_act(j)
|
|
||||||
j3=j+n_core_orb
|
|
||||||
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map)
|
|
||||||
do p=1,mo_num
|
|
||||||
do q=1,mo_num
|
|
||||||
bielec_PxxQ(p,i,j3,q)=integrals_array(p,q)
|
|
||||||
bielec_PxxQ(p,j3,i,q)=integrals_array(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='old',file='bielecCI.tmp')
|
||||||
! (ip|qj)
|
lread=.true.
|
||||||
do i=1,n_act_orb
|
indx=0
|
||||||
ii=list_act(i)
|
do while (lread)
|
||||||
i3=i+n_core_orb
|
read(12,*,iostat=kk) i,j,k,p,hhh
|
||||||
do j=i,n_act_orb
|
if (kk.ne.0) then
|
||||||
jj=list_act(j)
|
lread=.false.
|
||||||
j3=j+n_core_orb
|
else
|
||||||
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map)
|
bielecCI(i,j,k,p)=hhh
|
||||||
do p=1,mo_num
|
bielecCI(j,i,k,p)=hhh
|
||||||
do q=1,mo_num
|
indx+=1
|
||||||
bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q)
|
end if
|
||||||
bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p)
|
end do
|
||||||
end do
|
write(6,*) ' read ',indx,' integrals (aa|aP) into core '
|
||||||
end do
|
close(12)
|
||||||
end do
|
write(6,*) ' provided integrals (tu|xP) '
|
||||||
end do
|
|
||||||
write(6,*) ' provided integrals (Px|xQ) '
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
|
||||||
BEGIN_DOC
|
|
||||||
! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
|
|
||||||
! index p runs over the whole basis, t,u,v only over the active orbitals
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i,j,k,p,t,u,v
|
|
||||||
double precision, allocatable :: integrals_array(:)
|
|
||||||
real*8 :: mo_two_e_integral
|
|
||||||
|
|
||||||
allocate(integrals_array(mo_num))
|
|
||||||
|
|
||||||
do i=1,n_act_orb
|
|
||||||
t=list_act(i)
|
|
||||||
do j=1,n_act_orb
|
|
||||||
u=list_act(j)
|
|
||||||
do k=1,n_act_orb
|
|
||||||
v=list_act(k)
|
|
||||||
! (tu|vp)
|
|
||||||
call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map)
|
|
||||||
do p=1,mo_num
|
|
||||||
bielecCI(i,k,j,p)=integrals_array(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' provided integrals (tu|xP) '
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|||||||
118
src/casscf/bielec_create.irp.f
Normal file
118
src/casscf/bielec_create.irp.f
Normal file
@ -0,0 +1,118 @@
|
|||||||
|
! -*- F90 -*-
|
||||||
|
BEGIN_PROVIDER[real*8, bielec_PQxxtmp, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
|
||||||
|
&BEGIN_PROVIDER[real*8, bielec_PxxQtmp, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
|
||||||
|
&BEGIN_PROVIDER[integer, num_PQxx_written]
|
||||||
|
&BEGIN_PROVIDER[integer, num_PxxQ_written]
|
||||||
|
BEGIN_DOC
|
||||||
|
! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active
|
||||||
|
! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active
|
||||||
|
! indices are unshifted orbital numbers
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,ii,jj,p,q,i3,j3,t3,v3
|
||||||
|
double precision, allocatable :: integrals_array1(:,:)
|
||||||
|
double precision, allocatable :: integrals_array2(:,:)
|
||||||
|
real*8 :: mo_two_e_integral
|
||||||
|
|
||||||
|
allocate(integrals_array1(mo_num,mo_num))
|
||||||
|
allocate(integrals_array2(mo_num,mo_num))
|
||||||
|
|
||||||
|
do i=1,n_core_orb+n_act_orb
|
||||||
|
do j=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,mo_num
|
||||||
|
bielec_PQxxtmp(p,q,i,j)=0.D0
|
||||||
|
bielec_PxxQtmp(p,i,j,q)=0.D0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do i=1,n_core_orb
|
||||||
|
ii=list_core(i)
|
||||||
|
do j=i,n_core_orb
|
||||||
|
jj=list_core(j)
|
||||||
|
! (ij|pq)
|
||||||
|
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
|
||||||
|
! (ip|qj)
|
||||||
|
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map)
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,mo_num
|
||||||
|
bielec_PQxxtmp(p,q,i,j)=integrals_array1(p,q)
|
||||||
|
bielec_PQxxtmp(p,q,j,i)=integrals_array1(p,q)
|
||||||
|
bielec_PxxQtmp(p,i,j,q)=integrals_array2(p,q)
|
||||||
|
bielec_PxxQtmp(p,j,i,q)=integrals_array2(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do j=1,n_act_orb
|
||||||
|
jj=list_act(j)
|
||||||
|
j3=j+n_core_orb
|
||||||
|
! (ij|pq)
|
||||||
|
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
|
||||||
|
! (ip|qj)
|
||||||
|
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map)
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,mo_num
|
||||||
|
bielec_PQxxtmp(p,q,i,j3)=integrals_array1(p,q)
|
||||||
|
bielec_PQxxtmp(p,q,j3,i)=integrals_array1(p,q)
|
||||||
|
bielec_PxxQtmp(p,i,j3,q)=integrals_array2(p,q)
|
||||||
|
bielec_PxxQtmp(p,j3,i,q)=integrals_array2(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do i=1,n_act_orb
|
||||||
|
ii=list_act(i)
|
||||||
|
i3=i+n_core_orb
|
||||||
|
do j=i,n_act_orb
|
||||||
|
jj=list_act(j)
|
||||||
|
j3=j+n_core_orb
|
||||||
|
! (ij|pq)
|
||||||
|
call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map)
|
||||||
|
! (ip|qj)
|
||||||
|
call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map)
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,mo_num
|
||||||
|
bielec_PQxxtmp(p,q,i3,j3)=integrals_array1(p,q)
|
||||||
|
bielec_PQxxtmp(p,q,j3,i3)=integrals_array1(p,q)
|
||||||
|
bielec_PxxQtmp(p,i3,j3,q)=integrals_array2(p,q)
|
||||||
|
bielec_PxxQtmp(p,j3,i3,q)=integrals_array2(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' provided integrals (PQ|xx) '
|
||||||
|
write(6,*) ' provided integrals (Px|xQ) '
|
||||||
|
!!$ write(6,*) ' 1 1 1 2 = ',bielec_PQxxtmp(2,2,2,3),bielec_PxxQtmp(2,2,2,3)
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER[real*8, bielecCItmp, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! bielecCI : integrals (tu|vp) with p arbitrary, tuv active
|
||||||
|
! index p runs over the whole basis, t,u,v only over the active orbitals
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,k,p,t,u,v
|
||||||
|
double precision, allocatable :: integrals_array1(:)
|
||||||
|
real*8 :: mo_two_e_integral
|
||||||
|
|
||||||
|
allocate(integrals_array1(mo_num))
|
||||||
|
|
||||||
|
do i=1,n_act_orb
|
||||||
|
t=list_act(i)
|
||||||
|
do j=1,n_act_orb
|
||||||
|
u=list_act(j)
|
||||||
|
do k=1,n_act_orb
|
||||||
|
v=list_act(k)
|
||||||
|
! (tu|vp)
|
||||||
|
call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array1,mo_integrals_map)
|
||||||
|
do p=1,mo_num
|
||||||
|
bielecCItmp(i,k,j,p)=integrals_array1(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' provided integrals (tu|xP) '
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
@ -1,273 +0,0 @@
|
|||||||
BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)]
|
|
||||||
BEGIN_DOC
|
|
||||||
! integral (pq|xx) in the basis of natural MOs
|
|
||||||
! indices are unshifted orbital numbers
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i,j,k,l,t,u,p,q,pp
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
|
|
||||||
bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:)
|
|
||||||
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,n_core_orb+n_act_orb
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PQxx_no(list_act(p),j,k,l)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2nd quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,n_core_orb+n_act_orb
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PQxx_no(j,list_act(p),k,l)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 3rd quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,mo_num
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PQxx_no(j,k,n_core_orb+q,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PQxx_no(j,k,n_core_orb+p,l)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 4th quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,mo_num
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PQxx_no(j,k,l,n_core_orb+q)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PQxx_no(j,k,l,n_core_orb+p)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' transformed PQxx'
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)]
|
|
||||||
BEGIN_DOC
|
|
||||||
! integral (px|xq) in the basis of natural MOs
|
|
||||||
! indices are unshifted orbital numbers
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i,j,k,l,t,u,p,q,pp
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
|
|
||||||
bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:)
|
|
||||||
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,n_core_orb+n_act_orb
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PxxQ_no(list_act(p),k,l,j)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2nd quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,n_core_orb+n_act_orb
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PxxQ_no(j,k,l,list_act(p))=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 3rd quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,mo_num
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PxxQ_no(j,n_core_orb+q,l,k)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PxxQ_no(j,n_core_orb+p,l,k)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 4th quarter
|
|
||||||
do j=1,mo_num
|
|
||||||
do k=1,mo_num
|
|
||||||
do l=1,n_core_orb+n_act_orb
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielec_PxxQ_no(j,l,n_core_orb+q,k)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielec_PxxQ_no(j,l,n_core_orb+p,k)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' transformed PxxQ '
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
|
||||||
BEGIN_DOC
|
|
||||||
! integrals (tu|vp) in the basis of natural MOs
|
|
||||||
! index p runs over the whole basis, t,u,v only over the active orbitals
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i,j,k,l,t,u,p,q,pp
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
|
|
||||||
bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:)
|
|
||||||
|
|
||||||
do j=1,n_act_orb
|
|
||||||
do k=1,n_act_orb
|
|
||||||
do l=1,mo_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielecCI_no(p,j,k,l)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 2nd quarter
|
|
||||||
do j=1,n_act_orb
|
|
||||||
do k=1,n_act_orb
|
|
||||||
do l=1,mo_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielecCI_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,mo_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielecCI_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
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
bielecCI_no(j,k,l,list_act(p))=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' transformed tuvP '
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
@ -19,13 +19,7 @@ subroutine run
|
|||||||
N_det = 1
|
N_det = 1
|
||||||
TOUCH N_det psi_det psi_coef
|
TOUCH N_det psi_det psi_coef
|
||||||
call run_cipsi
|
call run_cipsi
|
||||||
|
call driver_wdens
|
||||||
write(6,*) ' total energy = ',eone+etwo+ecore
|
|
||||||
mo_label = "MCSCF"
|
|
||||||
mo_label = "Natural"
|
|
||||||
mo_coef(:,:) = NatOrbsFCI(:,:)
|
|
||||||
call save_mos
|
|
||||||
|
|
||||||
call driver_optorb
|
call driver_optorb
|
||||||
energy_old = energy
|
energy_old = energy
|
||||||
energy = eone+etwo+ecore
|
energy = eone+etwo+ecore
|
||||||
|
|||||||
@ -1,216 +1,177 @@
|
|||||||
use bitmasks
|
! -*- F90 -*-
|
||||||
|
use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
|
BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
|
||||||
BEGIN_DOC
|
&BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
|
||||||
! the first-order density matrix in the basis of the starting MOs
|
BEGIN_DOC
|
||||||
! matrices are state averaged
|
! the first-order density matrix in the basis of the starting MOs
|
||||||
!
|
! the second-order density matrix in the basis of the starting MOs
|
||||||
! we use the spin-free generators of mono-excitations
|
! matrices are state averaged
|
||||||
! E_pq destroys q and creates p
|
!
|
||||||
! D_pq = <0|E_pq|0> = D_qp
|
! we use the spin-free generators of mono-excitations
|
||||||
!
|
! E_pq destroys q and creates p
|
||||||
END_DOC
|
! D_pq = <0|E_pq|0> = D_qp
|
||||||
implicit none
|
! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
|
||||||
integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
|
!
|
||||||
integer :: ierr
|
END_DOC
|
||||||
integer(bit_kind) :: det_mu(N_int,2)
|
implicit none
|
||||||
integer(bit_kind) :: det_mu_ex(N_int,2)
|
integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
|
||||||
integer(bit_kind) :: det_mu_ex1(N_int,2)
|
integer :: ierr
|
||||||
integer(bit_kind) :: det_mu_ex2(N_int,2)
|
integer(bit_kind), allocatable :: det_mu(:,:)
|
||||||
real*8 :: phase1,phase2,term
|
integer(bit_kind), allocatable :: det_mu_ex(:,:)
|
||||||
integer :: nu1,nu2
|
integer(bit_kind), allocatable :: det_mu_ex1(:,:)
|
||||||
integer :: ierr1,ierr2
|
integer(bit_kind), allocatable :: det_mu_ex11(:,:)
|
||||||
real*8 :: cI_mu(N_states)
|
integer(bit_kind), allocatable :: det_mu_ex12(:,:)
|
||||||
|
integer(bit_kind), allocatable :: det_mu_ex2(:,:)
|
||||||
write(6,*) ' providing density matrices D0 and P0 '
|
integer(bit_kind), allocatable :: det_mu_ex21(:,:)
|
||||||
|
integer(bit_kind), allocatable :: det_mu_ex22(:,:)
|
||||||
D0tu = 0.d0
|
real*8 :: phase1,phase11,phase12,phase2,phase21,phase22
|
||||||
|
integer :: nu1,nu2,nu11,nu12,nu21,nu22
|
||||||
! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
|
integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22
|
||||||
do mu=1,n_det
|
real*8 :: cI_mu(N_states),term
|
||||||
call det_extract(det_mu,mu,N_int)
|
allocate(det_mu(N_int,2))
|
||||||
do istate=1,n_states
|
allocate(det_mu_ex(N_int,2))
|
||||||
cI_mu(istate)=psi_coef(mu,istate)
|
allocate(det_mu_ex1(N_int,2))
|
||||||
end do
|
allocate(det_mu_ex11(N_int,2))
|
||||||
do t=1,n_act_orb
|
allocate(det_mu_ex12(N_int,2))
|
||||||
ipart=list_act(t)
|
allocate(det_mu_ex2(N_int,2))
|
||||||
do u=1,n_act_orb
|
allocate(det_mu_ex21(N_int,2))
|
||||||
ihole=list_act(u)
|
allocate(det_mu_ex22(N_int,2))
|
||||||
! apply E_tu
|
|
||||||
call det_copy(det_mu,det_mu_ex1,N_int)
|
write(6,*) ' providing density matrices D0 and P0 '
|
||||||
call det_copy(det_mu,det_mu_ex2,N_int)
|
|
||||||
call do_spinfree_mono_excitation(det_mu,det_mu_ex1 &
|
! set all to zero
|
||||||
,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
|
do t=1,n_act_orb
|
||||||
! det_mu_ex1 is in the list
|
do u=1,n_act_orb
|
||||||
if (nu1.ne.-1) then
|
D0tu(u,t)=0.D0
|
||||||
do istate=1,n_states
|
do v=1,n_act_orb
|
||||||
term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
|
do x=1,n_act_orb
|
||||||
D0tu(t,u)+=term
|
P0tuvx(x,v,u,t)=0.D0
|
||||||
end do
|
end do
|
||||||
end if
|
|
||||||
! det_mu_ex2 is in the list
|
|
||||||
if (nu2.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
|
|
||||||
D0tu(t,u)+=term
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
! we average by just dividing by the number of states
|
|
||||||
do x=1,n_act_orb
|
|
||||||
do v=1,n_act_orb
|
|
||||||
D0tu(v,x)*=1.0D0/dble(N_states)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
|
|
||||||
BEGIN_DOC
|
|
||||||
! the second-order density matrix in the basis of the starting MOs
|
|
||||||
! matrices are state averaged
|
|
||||||
!
|
|
||||||
! we use the spin-free generators of mono-excitations
|
|
||||||
! E_pq destroys q and creates p
|
|
||||||
! D_pq = <0|E_pq|0> = D_qp
|
|
||||||
! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0>
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart
|
|
||||||
integer :: ierr
|
|
||||||
real*8 :: phase1,phase11,phase12,phase2,phase21,phase22
|
|
||||||
integer :: nu1,nu2,nu11,nu12,nu21,nu22
|
|
||||||
integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22
|
|
||||||
real*8 :: cI_mu(N_states),term
|
|
||||||
integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex
|
|
||||||
integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12
|
|
||||||
integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22
|
|
||||||
|
|
||||||
write(6,*) ' providing density matrices D0 and P0 '
|
|
||||||
|
|
||||||
P0tuvx = 0.d0
|
|
||||||
|
|
||||||
! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
|
|
||||||
do mu=1,n_det
|
|
||||||
call det_extract(det_mu,mu,N_int)
|
|
||||||
do istate=1,n_states
|
|
||||||
cI_mu(istate)=psi_coef(mu,istate)
|
|
||||||
end do
|
|
||||||
do t=1,n_act_orb
|
|
||||||
ipart=list_act(t)
|
|
||||||
do u=1,n_act_orb
|
|
||||||
ihole=list_act(u)
|
|
||||||
! apply E_tu
|
|
||||||
call det_copy(det_mu,det_mu_ex1,N_int)
|
|
||||||
call det_copy(det_mu,det_mu_ex2,N_int)
|
|
||||||
call do_spinfree_mono_excitation(det_mu,det_mu_ex1 &
|
|
||||||
,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
|
|
||||||
! det_mu_ex1 is in the list
|
|
||||||
if (nu1.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
|
|
||||||
! and we fill P0_tvvu
|
|
||||||
do v=1,n_act_orb
|
|
||||||
P0tuvx(t,v,v,u)-=term
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
! det_mu_ex2 is in the list
|
|
||||||
if (nu2.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
|
|
||||||
do v=1,n_act_orb
|
|
||||||
P0tuvx(t,v,v,u)-=term
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! now we do the double excitation E_tu E_vx |0>
|
|
||||||
do mu=1,n_det
|
|
||||||
call det_extract(det_mu,mu,N_int)
|
|
||||||
do istate=1,n_states
|
|
||||||
cI_mu(istate)=psi_coef(mu,istate)
|
|
||||||
end do
|
|
||||||
do v=1,n_act_orb
|
|
||||||
ipart=list_act(v)
|
|
||||||
do x=1,n_act_orb
|
|
||||||
ihole=list_act(x)
|
|
||||||
! apply E_vx
|
|
||||||
call det_copy(det_mu,det_mu_ex1,N_int)
|
|
||||||
call det_copy(det_mu,det_mu_ex2,N_int)
|
|
||||||
call do_spinfree_mono_excitation(det_mu,det_mu_ex1 &
|
|
||||||
,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
|
|
||||||
! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0>
|
|
||||||
if (ierr1.eq.1) then
|
|
||||||
do t=1,n_act_orb
|
|
||||||
jpart=list_act(t)
|
|
||||||
do u=1,n_act_orb
|
|
||||||
jhole=list_act(u)
|
|
||||||
call det_copy(det_mu_ex1,det_mu_ex11,N_int)
|
|
||||||
call det_copy(det_mu_ex1,det_mu_ex12,N_int)
|
|
||||||
call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11&
|
|
||||||
,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12)
|
|
||||||
if (nu11.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)&
|
|
||||||
*phase11*phase1
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
if (nu12.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)&
|
|
||||||
*phase12*phase1
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
|
|
||||||
! we apply E_tu to the second resultant determinant
|
|
||||||
if (ierr2.eq.1) then
|
|
||||||
do t=1,n_act_orb
|
|
||||||
jpart=list_act(t)
|
|
||||||
do u=1,n_act_orb
|
|
||||||
jhole=list_act(u)
|
|
||||||
call det_copy(det_mu_ex2,det_mu_ex21,N_int)
|
|
||||||
call det_copy(det_mu_ex2,det_mu_ex22,N_int)
|
|
||||||
call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21&
|
|
||||||
,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22)
|
|
||||||
if (nu21.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)&
|
|
||||||
*phase21*phase2
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
if (nu22.ne.-1) then
|
|
||||||
do istate=1,n_states
|
|
||||||
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)&
|
|
||||||
*phase22*phase2
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
! we average by just dividing by the number of states
|
|
||||||
do x=1,n_act_orb
|
|
||||||
do v=1,n_act_orb
|
|
||||||
do u=1,n_act_orb
|
|
||||||
do t=1,n_act_orb
|
|
||||||
P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)
|
|
||||||
end do
|
end do
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
end do
|
! first loop: we apply E_tu, once for D_tu, once for -P_tvvu
|
||||||
|
do mu=1,n_det
|
||||||
|
call det_extract(det_mu,mu,N_int)
|
||||||
|
do istate=1,n_states
|
||||||
|
cI_mu(istate)=psi_coef(mu,istate)
|
||||||
|
end do
|
||||||
|
do t=1,n_act_orb
|
||||||
|
ipart=list_act(t)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
ihole=list_act(u)
|
||||||
|
! apply E_tu
|
||||||
|
call det_copy(det_mu,det_mu_ex1,N_int)
|
||||||
|
call det_copy(det_mu,det_mu_ex2,N_int)
|
||||||
|
call do_spinfree_mono_excitation(det_mu,det_mu_ex1 &
|
||||||
|
,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
|
||||||
|
! det_mu_ex1 is in the list
|
||||||
|
if (nu1.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
term=cI_mu(istate)*psi_coef(nu1,istate)*phase1
|
||||||
|
D0tu(t,u)+=term
|
||||||
|
! and we fill P0_tvvu
|
||||||
|
do v=1,n_act_orb
|
||||||
|
P0tuvx(t,v,v,u)-=term
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
! det_mu_ex2 is in the list
|
||||||
|
if (nu2.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
term=cI_mu(istate)*psi_coef(nu2,istate)*phase2
|
||||||
|
D0tu(t,u)+=term
|
||||||
|
do v=1,n_act_orb
|
||||||
|
P0tuvx(t,v,v,u)-=term
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! now we do the double excitation E_tu E_vx |0>
|
||||||
|
do mu=1,n_det
|
||||||
|
call det_extract(det_mu,mu,N_int)
|
||||||
|
do istate=1,n_states
|
||||||
|
cI_mu(istate)=psi_coef(mu,istate)
|
||||||
|
end do
|
||||||
|
do v=1,n_act_orb
|
||||||
|
ipart=list_act(v)
|
||||||
|
do x=1,n_act_orb
|
||||||
|
ihole=list_act(x)
|
||||||
|
! apply E_vx
|
||||||
|
call det_copy(det_mu,det_mu_ex1,N_int)
|
||||||
|
call det_copy(det_mu,det_mu_ex2,N_int)
|
||||||
|
call do_spinfree_mono_excitation(det_mu,det_mu_ex1 &
|
||||||
|
,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2)
|
||||||
|
! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0>
|
||||||
|
if (ierr1.eq.1) then
|
||||||
|
do t=1,n_act_orb
|
||||||
|
jpart=list_act(t)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
jhole=list_act(u)
|
||||||
|
call det_copy(det_mu_ex1,det_mu_ex11,N_int)
|
||||||
|
call det_copy(det_mu_ex1,det_mu_ex12,N_int)
|
||||||
|
call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11 &
|
||||||
|
,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12)
|
||||||
|
if (nu11.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate) &
|
||||||
|
*phase11*phase1
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
if (nu12.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate) &
|
||||||
|
*phase12*phase1
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
|
||||||
|
! we apply E_tu to the second resultant determinant
|
||||||
|
if (ierr2.eq.1) then
|
||||||
|
do t=1,n_act_orb
|
||||||
|
jpart=list_act(t)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
jhole=list_act(u)
|
||||||
|
call det_copy(det_mu_ex2,det_mu_ex21,N_int)
|
||||||
|
call det_copy(det_mu_ex2,det_mu_ex22,N_int)
|
||||||
|
call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21 &
|
||||||
|
,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22)
|
||||||
|
if (nu21.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate) &
|
||||||
|
*phase21*phase2
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
if (nu22.ne.-1) then
|
||||||
|
do istate=1,n_states
|
||||||
|
P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate) &
|
||||||
|
*phase22*phase2
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
! we average by just dividing by the number of states
|
||||||
|
do x=1,n_act_orb
|
||||||
|
do v=1,n_act_orb
|
||||||
|
D0tu(v,x)*=1.0D0/dble(N_states)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
do t=1,n_act_orb
|
||||||
|
P0tuvx(t,u,v,x)*=0.5D0/dble(N_states)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|||||||
@ -1,130 +1,131 @@
|
|||||||
use bitmasks
|
! -*- F90 -*-
|
||||||
|
use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
|
|
||||||
subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, &
|
||||||
ispin,phase,ierr)
|
ispin,phase,ierr)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! we create the mono-excitation, and determine, if possible,
|
! we create the mono-excitation, and determine, if possible,
|
||||||
! the phase and the number in the list of determinants
|
! the phase and the number in the list of determinants
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer(bit_kind) :: key1(N_int,2),key2(N_int,2)
|
integer(bit_kind) :: key1(N_int,2),key2(N_int,2)
|
||||||
integer(bit_kind), allocatable :: keytmp(:,:)
|
integer(bit_kind), allocatable :: keytmp(:,:)
|
||||||
integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin
|
integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin
|
||||||
real*8 :: phase
|
real*8 :: phase
|
||||||
logical :: found
|
logical :: found
|
||||||
allocate(keytmp(N_int,2))
|
allocate(keytmp(N_int,2))
|
||||||
|
|
||||||
nu=-1
|
nu=-1
|
||||||
phase=1.D0
|
phase=1.D0
|
||||||
ierr=0
|
ierr=0
|
||||||
call det_copy(key1,key2,N_int)
|
call det_copy(key1,key2,N_int)
|
||||||
! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin
|
! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin
|
||||||
! call print_det(key2,N_int)
|
! call print_det(key2,N_int)
|
||||||
call do_single_excitation(key2,ihole,ipart,ispin,ierr)
|
call do_single_excitation(key2,ihole,ipart,ispin,ierr)
|
||||||
! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin
|
! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin
|
||||||
! call print_det(key2,N_int)
|
! call print_det(key2,N_int)
|
||||||
! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr
|
! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr
|
||||||
if (ierr.eq.1) then
|
if (ierr.eq.1) then
|
||||||
! excitation is possible
|
! excitation is possible
|
||||||
! get the phase
|
! get the phase
|
||||||
call get_single_excitation(key1,key2,exc,phase,N_int)
|
call get_single_excitation(key1,key2,exc,phase,N_int)
|
||||||
! get the number in the list
|
! get the number in the list
|
||||||
found=.false.
|
found=.false.
|
||||||
nu=0
|
nu=0
|
||||||
do while (.not.found)
|
do while (.not.found)
|
||||||
nu+=1
|
nu+=1
|
||||||
if (nu.gt.N_det) then
|
if (nu.gt.N_det) then
|
||||||
! the determinant is possible, but not in the list
|
! the determinant is possible, but not in the list
|
||||||
found=.true.
|
found=.true.
|
||||||
nu=-1
|
nu=-1
|
||||||
else
|
else
|
||||||
call det_extract(keytmp,nu,N_int)
|
call det_extract(keytmp,nu,N_int)
|
||||||
integer :: i,ii
|
integer :: i,ii
|
||||||
found=.true.
|
found=.true.
|
||||||
do ii=1,2
|
do ii=1,2
|
||||||
do i=1,N_int
|
do i=1,N_int
|
||||||
if (keytmp(i,ii).ne.key2(i,ii)) then
|
if (keytmp(i,ii).ne.key2(i,ii)) then
|
||||||
found=.false.
|
found=.false.
|
||||||
end if
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end if
|
||||||
end do
|
end do
|
||||||
end do
|
! if (found) then
|
||||||
end if
|
! if (nu.eq.-1) then
|
||||||
end do
|
! write(6,*) ' image not found in the list, thus nu = ',nu
|
||||||
! if (found) then
|
! else
|
||||||
! if (nu.eq.-1) then
|
! write(6,*) ' found in the list as No ',nu,' phase = ',phase
|
||||||
! write(6,*) ' image not found in the list, thus nu = ',nu
|
! end if
|
||||||
! else
|
! end if
|
||||||
! write(6,*) ' found in the list as No ',nu,' phase = ',phase
|
end if
|
||||||
! end if
|
!
|
||||||
! end if
|
! we found the new string, the phase, and possibly the number in the list
|
||||||
end if
|
!
|
||||||
!
|
end subroutine do_signed_mono_excitation
|
||||||
! we found the new string, the phase, and possibly the number in the list
|
|
||||||
!
|
|
||||||
end subroutine do_signed_mono_excitation
|
|
||||||
|
|
||||||
subroutine det_extract(key,nu,Nint)
|
subroutine det_extract(key,nu,Nint)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! extract a determinant from the list of determinants
|
! extract a determinant from the list of determinants
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ispin,i,nu,Nint
|
integer :: ispin,i,nu,Nint
|
||||||
integer(bit_kind) :: key(Nint,2)
|
integer(bit_kind) :: key(Nint,2)
|
||||||
do ispin=1,2
|
do ispin=1,2
|
||||||
do i=1,Nint
|
do i=1,Nint
|
||||||
key(i,ispin)=psi_det(i,ispin,nu)
|
key(i,ispin)=psi_det(i,ispin,nu)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end subroutine det_extract
|
end subroutine det_extract
|
||||||
|
|
||||||
subroutine det_copy(key1,key2,Nint)
|
subroutine det_copy(key1,key2,Nint)
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! copy a determinant from key1 to key2
|
! copy a determinant from key1 to key2
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ispin,i,Nint
|
integer :: ispin,i,Nint
|
||||||
integer(bit_kind) :: key1(Nint,2),key2(Nint,2)
|
integer(bit_kind) :: key1(Nint,2),key2(Nint,2)
|
||||||
do ispin=1,2
|
do ispin=1,2
|
||||||
do i=1,Nint
|
do i=1,Nint
|
||||||
key2(i,ispin)=key1(i,ispin)
|
key2(i,ispin)=key1(i,ispin)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
end subroutine det_copy
|
end subroutine det_copy
|
||||||
|
|
||||||
|
subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 &
|
||||||
|
,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr)
|
||||||
|
BEGIN_DOC
|
||||||
|
! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q)
|
||||||
|
! we may create two determinants as result
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2)
|
||||||
|
integer(bit_kind) :: key_out2(N_int,2)
|
||||||
|
integer :: ihole,ipart,ierr,jerr,nu1,nu2
|
||||||
|
integer :: ispin
|
||||||
|
real*8 :: phase1,phase2
|
||||||
|
|
||||||
|
! write(6,*) ' applying E_',ipart,ihole,' on determinant '
|
||||||
|
! call print_det(key_in,N_int)
|
||||||
|
|
||||||
|
! spin alpha
|
||||||
|
ispin=1
|
||||||
|
call do_signed_mono_excitation(key_in,key_out1,nu1,ihole &
|
||||||
|
,ipart,ispin,phase1,ierr)
|
||||||
|
! if (ierr.eq.1) then
|
||||||
|
! write(6,*) ' 1 result is ',nu1,phase1
|
||||||
|
! call print_det(key_out1,N_int)
|
||||||
|
! end if
|
||||||
|
! spin beta
|
||||||
|
ispin=2
|
||||||
|
call do_signed_mono_excitation(key_in,key_out2,nu2,ihole &
|
||||||
|
,ipart,ispin,phase2,jerr)
|
||||||
|
! if (jerr.eq.1) then
|
||||||
|
! write(6,*) ' 2 result is ',nu2,phase2
|
||||||
|
! call print_det(key_out2,N_int)
|
||||||
|
! end if
|
||||||
|
|
||||||
subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 &
|
end subroutine do_spinfree_mono_excitation
|
||||||
,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr)
|
|
||||||
BEGIN_DOC
|
|
||||||
! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q)
|
|
||||||
! we may create two determinants as result
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2)
|
|
||||||
integer(bit_kind) :: key_out2(N_int,2)
|
|
||||||
integer :: ihole,ipart,ierr,jerr,nu1,nu2
|
|
||||||
integer :: ispin
|
|
||||||
real*8 :: phase1,phase2
|
|
||||||
|
|
||||||
! write(6,*) ' applying E_',ipart,ihole,' on determinant '
|
|
||||||
! call print_det(key_in,N_int)
|
|
||||||
|
|
||||||
! spin alpha
|
|
||||||
ispin=1
|
|
||||||
call do_signed_mono_excitation(key_in,key_out1,nu1,ihole &
|
|
||||||
,ipart,ispin,phase1,ierr)
|
|
||||||
! if (ierr.eq.1) then
|
|
||||||
! write(6,*) ' 1 result is ',nu1,phase1
|
|
||||||
! call print_det(key_out1,N_int)
|
|
||||||
! end if
|
|
||||||
! spin beta
|
|
||||||
ispin=2
|
|
||||||
call do_signed_mono_excitation(key_in,key_out2,nu2,ihole &
|
|
||||||
,ipart,ispin,phase2,jerr)
|
|
||||||
! if (jerr.eq.1) then
|
|
||||||
! write(6,*) ' 2 result is ',nu2,phase2
|
|
||||||
! call print_det(key_out2,N_int)
|
|
||||||
! end if
|
|
||||||
|
|
||||||
end subroutine do_spinfree_mono_excitation
|
|
||||||
|
|
||||||
|
|||||||
154
src/casscf/driver_wdens.irp.f
Normal file
154
src/casscf/driver_wdens.irp.f
Normal file
@ -0,0 +1,154 @@
|
|||||||
|
subroutine driver_wdens
|
||||||
|
implicit none
|
||||||
|
integer :: istate,p,q,r,s,indx,i,j
|
||||||
|
|
||||||
|
|
||||||
|
write(6,*) ' total energy = ',eone+etwo+ecore
|
||||||
|
write(6,*) ' generating natural orbitals '
|
||||||
|
write(6,*)
|
||||||
|
write(6,*)
|
||||||
|
call trf_to_natorb
|
||||||
|
|
||||||
|
write(6,*) ' all data available ! '
|
||||||
|
write(6,*) ' writing out files '
|
||||||
|
|
||||||
|
open(unit=12,file='D0tu.dat',form='formatted',status='unknown')
|
||||||
|
do p=1,n_act_orb
|
||||||
|
do q=1,n_act_orb
|
||||||
|
if (abs(D0tu(p,q)).gt.1.D-12) then
|
||||||
|
write(12,'(2i8,E20.12)') p,q,D0tu(p,q)
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
real*8 :: approx,np,nq,nr,ns
|
||||||
|
logical :: lpq,lrs,lps,lqr
|
||||||
|
open(unit=12,file='P0tuvx.dat',form='formatted',status='unknown')
|
||||||
|
do p=1,n_act_orb
|
||||||
|
np=D0tu(p,p)
|
||||||
|
do q=1,n_act_orb
|
||||||
|
lpq=p.eq.q
|
||||||
|
nq=D0tu(q,q)
|
||||||
|
do r=1,n_act_orb
|
||||||
|
lqr=q.eq.r
|
||||||
|
nr=D0tu(r,r)
|
||||||
|
do s=1,n_act_orb
|
||||||
|
lrs=r.eq.s
|
||||||
|
lps=p.eq.s
|
||||||
|
approx=0.D0
|
||||||
|
if (lpq.and.lrs) then
|
||||||
|
if (lqr) then
|
||||||
|
! pppp
|
||||||
|
approx=0.5D0*np*(np-1.D0)
|
||||||
|
else
|
||||||
|
! pprr
|
||||||
|
approx=0.5D0*np*nr
|
||||||
|
end if
|
||||||
|
else
|
||||||
|
if (lps.and.lqr.and..not.lpq) then
|
||||||
|
! pqqp
|
||||||
|
approx=-0.25D0*np*nq
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
if (abs(P0tuvx(p,q,r,s)).gt.1.D-12) then
|
||||||
|
write(12,'(4i4,2E20.12)') p,q,r,s,P0tuvx(p,q,r,s),approx
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='unknown',file='onetrf.tmp')
|
||||||
|
indx=0
|
||||||
|
do q=1,mo_num
|
||||||
|
do p=q,mo_num
|
||||||
|
if (abs(onetrf(p,q)).gt.1.D-12) then
|
||||||
|
write(12,'(2i6,E20.12)') p,q,onetrf(p,q)
|
||||||
|
indx+=1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' wrote ',indx,' mono-electronic integrals'
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='unknown',file='bielec_PQxx.tmp')
|
||||||
|
indx=0
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=p,mo_num
|
||||||
|
do r=1,n_core_orb+n_act_orb
|
||||||
|
do s=r,n_core_orb+n_act_orb
|
||||||
|
if (abs(bielec_PQxxtmp(p,q,r,s)).gt.1.D-12) then
|
||||||
|
write(12,'(4i8,E20.12)') p,q,r,s,bielec_PQxxtmp(p,q,r,s)
|
||||||
|
indx+=1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' wrote ',indx,' integrals (PQ|xx)'
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='unknown',file='bielec_PxxQ.tmp')
|
||||||
|
indx=0
|
||||||
|
do p=1,mo_num
|
||||||
|
do q=1,n_core_orb+n_act_orb
|
||||||
|
do r=q,n_core_orb+n_act_orb
|
||||||
|
integer ::s_start
|
||||||
|
if (q.eq.r) then
|
||||||
|
s_start=p
|
||||||
|
else
|
||||||
|
s_start=1
|
||||||
|
end if
|
||||||
|
do s=s_start,mo_num
|
||||||
|
if (abs(bielec_PxxQtmp(p,q,r,s)).gt.1.D-12) then
|
||||||
|
write(12,'(4i8,E20.12)') p,q,r,s,bielec_PxxQtmp(p,q,r,s)
|
||||||
|
indx+=1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' wrote ',indx,' integrals (Px|xQ)'
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
open(unit=12,form='formatted',status='unknown',file='bielecCI.tmp')
|
||||||
|
indx=0
|
||||||
|
do p=1,n_act_orb
|
||||||
|
do q=p,n_act_orb
|
||||||
|
do r=1,n_act_orb
|
||||||
|
do s=1,mo_num
|
||||||
|
if (abs(bielecCItmp(p,q,r,s)).gt.1.D-12) then
|
||||||
|
write(12,'(4i8,E20.12)') p,q,r,s,bielecCItmp(p,q,r,s)
|
||||||
|
indx+=1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' wrote ',indx,' integrals (tu|xP)'
|
||||||
|
close(12)
|
||||||
|
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) ' creating new orbitals '
|
||||||
|
do i=1,mo_num
|
||||||
|
write(6,*) ' Orbital No ',i
|
||||||
|
write(6,'(5F14.6)') (NatOrbsFCI(j,i),j=1,mo_num)
|
||||||
|
write(6,*)
|
||||||
|
end do
|
||||||
|
|
||||||
|
mo_label = "MCSCF"
|
||||||
|
mo_label = "Natural"
|
||||||
|
do i=1,mo_num
|
||||||
|
do j=1,ao_num
|
||||||
|
mo_coef(j,i)=NatOrbsFCI(j,i)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
call save_mos
|
||||||
|
|
||||||
|
write(6,*) ' ... done '
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
@ -1,249 +1,251 @@
|
|||||||
use bitmasks
|
! -*- F90 -*-
|
||||||
|
|
||||||
|
use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
|
|
||||||
BEGIN_PROVIDER [ integer, nMonoEx ]
|
BEGIN_PROVIDER [ integer, nMonoEx ]
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Number of single excitations
|
!
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
|
nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb
|
||||||
write(6,*) ' nMonoEx = ',nMonoEx
|
write(6,*) ' nMonoEx = ',nMonoEx
|
||||||
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)]
|
||||||
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
|
integer :: i,t,a,ii,tt,aa,indx
|
||||||
indx=0
|
indx=0
|
||||||
do ii=1,n_core_orb
|
do ii=1,n_core_orb
|
||||||
i=list_core(ii)
|
i=list_core(ii)
|
||||||
do tt=1,n_act_orb
|
do tt=1,n_act_orb
|
||||||
|
t=list_act(tt)
|
||||||
|
indx+=1
|
||||||
|
excit(1,indx)=i
|
||||||
|
excit(2,indx)=t
|
||||||
|
excit_class(indx)='c-a'
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do ii=1,n_core_orb
|
||||||
|
i=list_core(ii)
|
||||||
|
do aa=1,n_virt_orb
|
||||||
|
a=list_virt(aa)
|
||||||
|
indx+=1
|
||||||
|
excit(1,indx)=i
|
||||||
|
excit(2,indx)=a
|
||||||
|
excit_class(indx)='c-v'
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do tt=1,n_act_orb
|
||||||
t=list_act(tt)
|
t=list_act(tt)
|
||||||
indx+=1
|
do aa=1,n_virt_orb
|
||||||
excit(1,indx)=i
|
a=list_virt(aa)
|
||||||
excit(2,indx)=t
|
indx+=1
|
||||||
excit_class(indx)='c-a'
|
excit(1,indx)=t
|
||||||
end do
|
excit(2,indx)=a
|
||||||
end do
|
excit_class(indx)='a-v'
|
||||||
|
end do
|
||||||
do ii=1,n_core_orb
|
end do
|
||||||
i=list_core(ii)
|
|
||||||
do aa=1,n_virt_orb
|
if (bavard) then
|
||||||
a=list_virt(aa)
|
write(6,*) ' Filled the table of the Monoexcitations '
|
||||||
indx+=1
|
do indx=1,nMonoEx
|
||||||
excit(1,indx)=i
|
write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' &
|
||||||
excit(2,indx)=a
|
,excit(2,indx),' ',excit_class(indx)
|
||||||
excit_class(indx)='c-v'
|
end do
|
||||||
end do
|
end if
|
||||||
end do
|
|
||||||
|
|
||||||
do tt=1,n_act_orb
|
|
||||||
t=list_act(tt)
|
|
||||||
do aa=1,n_virt_orb
|
|
||||||
a=list_virt(aa)
|
|
||||||
indx+=1
|
|
||||||
excit(1,indx)=t
|
|
||||||
excit(2,indx)=a
|
|
||||||
excit_class(indx)='a-v'
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
if (bavard) then
|
|
||||||
write(6,*) ' Filled the table of the Monoexcitations '
|
|
||||||
do indx=1,nMonoEx
|
|
||||||
write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' &
|
|
||||||
,excit(2,indx),' ',excit_class(indx)
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
|
BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)]
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
|
! calculate the orbital gradient <Psi| H E_pq |Psi> by hand, i.e. for
|
||||||
! each determinant I we determine the string E_pq |I> (alpha and beta
|
! each determinant I we determine the string E_pq |I> (alpha and beta
|
||||||
! separately) and generate <Psi|H E_pq |I>
|
! separately) and generate <Psi|H E_pq |I>
|
||||||
! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
|
! sum_I c_I <Psi|H E_pq |I> is then the pq component of the orbital
|
||||||
! gradient
|
! gradient
|
||||||
! E_pq = a^+_pa_q + a^+_Pa_Q
|
! E_pq = a^+_pa_q + a^+_Pa_Q
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ii,tt,aa,indx,ihole,ipart,istate
|
integer :: ii,tt,aa,indx,ihole,ipart,istate
|
||||||
real*8 :: res
|
real*8 :: res
|
||||||
|
|
||||||
do indx=1,nMonoEx
|
do indx=1,nMonoEx
|
||||||
ihole=excit(1,indx)
|
ihole=excit(1,indx)
|
||||||
ipart=excit(2,indx)
|
ipart=excit(2,indx)
|
||||||
call calc_grad_elem(ihole,ipart,res)
|
call calc_grad_elem(ihole,ipart,res)
|
||||||
gradvec(indx)=res
|
gradvec(indx)=res
|
||||||
end do
|
end do
|
||||||
|
|
||||||
real*8 :: norm_grad
|
real*8 :: norm_grad
|
||||||
norm_grad=0.d0
|
norm_grad=0.d0
|
||||||
do indx=1,nMonoEx
|
do indx=1,nMonoEx
|
||||||
norm_grad+=gradvec(indx)*gradvec(indx)
|
norm_grad+=gradvec(indx)*gradvec(indx)
|
||||||
end do
|
end do
|
||||||
norm_grad=sqrt(norm_grad)
|
norm_grad=sqrt(norm_grad)
|
||||||
write(6,*)
|
write(6,*)
|
||||||
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad
|
||||||
write(6,*)
|
write(6,*)
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
subroutine calc_grad_elem(ihole,ipart,res)
|
subroutine calc_grad_elem(ihole,ipart,res)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! eq 18 of Siegbahn et al, Physica Scripta 1980
|
! eq 18 of Siegbahn et al, Physica Scripta 1980
|
||||||
! we calculate 2 <Psi| H E_pq | Psi>, q=hole, p=particle
|
! we calculate 2 <Psi| H E_pq | Psi>, q=hole, p=particle
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate
|
integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate
|
||||||
real*8 :: res
|
real*8 :: res
|
||||||
integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
|
integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:)
|
||||||
real*8 :: i_H_psi_array(N_states),phase
|
real*8 :: i_H_psi_array(N_states),phase
|
||||||
allocate(det_mu(N_int,2))
|
allocate(det_mu(N_int,2))
|
||||||
allocate(det_mu_ex(N_int,2))
|
allocate(det_mu_ex(N_int,2))
|
||||||
|
|
||||||
res=0.D0
|
res=0.D0
|
||||||
|
|
||||||
do mu=1,n_det
|
do mu=1,n_det
|
||||||
! get the string of the determinant
|
! get the string of the determinant
|
||||||
call det_extract(det_mu,mu,N_int)
|
call det_extract(det_mu,mu,N_int)
|
||||||
do ispin=1,2
|
do ispin=1,2
|
||||||
! do the monoexcitation on it
|
! do the monoexcitation on it
|
||||||
call det_copy(det_mu,det_mu_ex,N_int)
|
call det_copy(det_mu,det_mu_ex,N_int)
|
||||||
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
|
call do_signed_mono_excitation(det_mu,det_mu_ex,nu &
|
||||||
,ihole,ipart,ispin,phase,ierr)
|
,ihole,ipart,ispin,phase,ierr)
|
||||||
if (ierr.eq.1) then
|
if (ierr.eq.1) then
|
||||||
! write(6,*)
|
! write(6,*)
|
||||||
! write(6,*) ' mu = ',mu
|
! write(6,*) ' mu = ',mu
|
||||||
! call print_det(det_mu,N_int)
|
! call print_det(det_mu,N_int)
|
||||||
! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin
|
! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin
|
||||||
! call print_det(det_mu_ex,N_int)
|
! call print_det(det_mu_ex,N_int)
|
||||||
call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int &
|
call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int &
|
||||||
,N_det,N_det,N_states,i_H_psi_array)
|
,N_det,N_det,N_states,i_H_psi_array)
|
||||||
do istate=1,N_states
|
do istate=1,N_states
|
||||||
res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
|
res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase
|
||||||
|
end do
|
||||||
|
! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res
|
||||||
|
end if
|
||||||
end do
|
end do
|
||||||
! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res
|
end do
|
||||||
end if
|
|
||||||
end do
|
! state-averaged gradient
|
||||||
end do
|
res*=2.D0/dble(N_states)
|
||||||
|
|
||||||
! state-averaged gradient
|
end subroutine calc_grad_elem
|
||||||
res*=2.D0/dble(N_states)
|
|
||||||
|
|
||||||
end subroutine calc_grad_elem
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
|
BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)]
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! calculate the orbital gradient <Psi| H E_pq |Psi> from density
|
! calculate the orbital gradient <Psi| H E_pq |Psi> from density
|
||||||
! matrices and integrals; Siegbahn et al, Phys Scr 1980
|
! matrices and integrals; Siegbahn et al, Phys Scr 1980
|
||||||
! eqs 14 a,b,c
|
! eqs 14 a,b,c
|
||||||
END_DOC
|
END_DOC
|
||||||
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
|
||||||
real*8 :: norm_grad
|
real*8 :: norm_grad
|
||||||
|
|
||||||
indx=0
|
indx=0
|
||||||
do i=1,n_core_orb
|
do i=1,n_core_orb
|
||||||
do t=1,n_act_orb
|
do t=1,n_act_orb
|
||||||
indx+=1
|
indx+=1
|
||||||
gradvec2(indx)=gradvec_it(i,t)
|
gradvec2(indx)=gradvec_it(i,t)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
do i=1,n_core_orb
|
do i=1,n_core_orb
|
||||||
do a=1,n_virt_orb
|
do a=1,n_virt_orb
|
||||||
indx+=1
|
indx+=1
|
||||||
gradvec2(indx)=gradvec_ia(i,a)
|
gradvec2(indx)=gradvec_ia(i,a)
|
||||||
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
|
||||||
gradvec2(indx)=gradvec_ta(t,a)
|
gradvec2(indx)=gradvec_ta(t,a)
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
|
|
||||||
norm_grad=0.d0
|
norm_grad=0.d0
|
||||||
do indx=1,nMonoEx
|
do indx=1,nMonoEx
|
||||||
norm_grad+=gradvec2(indx)*gradvec2(indx)
|
norm_grad+=gradvec2(indx)*gradvec2(indx)
|
||||||
end do
|
end do
|
||||||
norm_grad=sqrt(norm_grad)
|
norm_grad=sqrt(norm_grad)
|
||||||
write(6,*)
|
write(6,*)
|
||||||
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad
|
||||||
write(6,*)
|
write(6,*)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
real*8 function gradvec_it(i,t)
|
real*8 function gradvec_it(i,t)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! the orbital gradient core -> active
|
! the orbital gradient core -> active
|
||||||
! we assume natural orbitals
|
! we assume natural orbitals
|
||||||
END_DOC
|
END_DOC
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,t
|
integer :: i,t
|
||||||
|
|
||||||
integer :: ii,tt,v,vv,x,y
|
|
||||||
integer :: x3,y3
|
|
||||||
|
|
||||||
ii=list_core(i)
|
|
||||||
tt=list_act(t)
|
|
||||||
gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
|
|
||||||
gradvec_it-=occnum(tt)*Fipq(ii,tt)
|
|
||||||
do v=1,n_act_orb
|
|
||||||
vv=list_act(v)
|
|
||||||
do x=1,n_act_orb
|
|
||||||
x3=x+n_core_orb
|
|
||||||
do y=1,n_act_orb
|
|
||||||
y3=y+n_core_orb
|
|
||||||
gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
gradvec_it*=2.D0
|
|
||||||
end function gradvec_it
|
|
||||||
|
|
||||||
real*8 function gradvec_ia(i,a)
|
integer :: ii,tt,v,vv,x,y
|
||||||
BEGIN_DOC
|
integer :: x3,y3
|
||||||
! the orbital gradient core -> virtual
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: i,a,ii,aa
|
|
||||||
|
|
||||||
ii=list_core(i)
|
|
||||||
aa=list_virt(a)
|
|
||||||
gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii))
|
|
||||||
gradvec_ia*=2.D0
|
|
||||||
|
|
||||||
end function gradvec_ia
|
|
||||||
|
|
||||||
real*8 function gradvec_ta(t,a)
|
ii=list_core(i)
|
||||||
BEGIN_DOC
|
tt=list_act(t)
|
||||||
! the orbital gradient active -> virtual
|
gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii))
|
||||||
! we assume natural orbitals
|
gradvec_it-=occnum(tt)*Fipq(ii,tt)
|
||||||
END_DOC
|
do v=1,n_act_orb
|
||||||
implicit none
|
vv=list_act(v)
|
||||||
integer :: t,a,tt,aa,v,vv,x,y
|
do x=1,n_act_orb
|
||||||
|
x3=x+n_core_orb
|
||||||
tt=list_act(t)
|
do y=1,n_act_orb
|
||||||
aa=list_virt(a)
|
y3=y+n_core_orb
|
||||||
gradvec_ta=0.D0
|
gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx(ii,vv,x3,y3)
|
||||||
gradvec_ta+=occnum(tt)*Fipq(aa,tt)
|
end do
|
||||||
do v=1,n_act_orb
|
end do
|
||||||
do x=1,n_act_orb
|
end do
|
||||||
do y=1,n_act_orb
|
gradvec_it*=2.D0
|
||||||
gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa)
|
end function gradvec_it
|
||||||
end do
|
|
||||||
end do
|
real*8 function gradvec_ia(i,a)
|
||||||
end do
|
BEGIN_DOC
|
||||||
gradvec_ta*=2.D0
|
! the orbital gradient core -> virtual
|
||||||
|
END_DOC
|
||||||
end function gradvec_ta
|
implicit none
|
||||||
|
integer :: i,a,ii,aa
|
||||||
|
|
||||||
|
ii=list_core(i)
|
||||||
|
aa=list_virt(a)
|
||||||
|
gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii))
|
||||||
|
gradvec_ia*=2.D0
|
||||||
|
|
||||||
|
end function gradvec_ia
|
||||||
|
|
||||||
|
real*8 function gradvec_ta(t,a)
|
||||||
|
BEGIN_DOC
|
||||||
|
! the orbital gradient active -> virtual
|
||||||
|
! we assume natural orbitals
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: t,a,tt,aa,v,vv,x,y
|
||||||
|
|
||||||
|
tt=list_act(t)
|
||||||
|
aa=list_virt(a)
|
||||||
|
gradvec_ta=0.D0
|
||||||
|
gradvec_ta+=occnum(tt)*Fipq(aa,tt)
|
||||||
|
do v=1,n_act_orb
|
||||||
|
do x=1,n_act_orb
|
||||||
|
do y=1,n_act_orb
|
||||||
|
gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
gradvec_ta*=2.D0
|
||||||
|
|
||||||
|
end function gradvec_ta
|
||||||
|
|
||||||
|
|||||||
File diff suppressed because it is too large
Load Diff
@ -1,80 +1,67 @@
|
|||||||
BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
|
! -*- F90 -*-
|
||||||
BEGIN_DOC
|
BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ]
|
||||||
! the inactive Fock matrix, in molecular orbitals
|
&BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
|
||||||
END_DOC
|
BEGIN_DOC
|
||||||
implicit none
|
! the inactive and the active Fock matrices, in molecular
|
||||||
integer :: p,q,k,kk,t,tt,u,uu
|
! orbitals
|
||||||
|
! we create them in MOs, quite expensive
|
||||||
do q=1,mo_num
|
!
|
||||||
do p=1,mo_num
|
! for an implementation in AOs we need first the natural orbitals
|
||||||
Fipq(p,q)=one_ints_no(p,q)
|
! for forming an active density matrix in AOs
|
||||||
end do
|
!
|
||||||
end do
|
END_DOC
|
||||||
|
implicit none
|
||||||
! the inactive Fock matrix
|
double precision, allocatable :: integrals_array1(:,:)
|
||||||
do k=1,n_core_orb
|
double precision, allocatable :: integrals_array2(:,:)
|
||||||
kk=list_core(k)
|
integer :: p,q,k,kk,t,tt,u,uu
|
||||||
do q=1,mo_num
|
allocate(integrals_array1(mo_num,mo_num))
|
||||||
|
allocate(integrals_array2(mo_num,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 q=1,mo_num
|
||||||
|
Fipq(p,q)=one_ints(p,q)
|
||||||
|
Fapq(p,q)=0.D0
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
end do
|
! the inactive Fock matrix
|
||||||
|
do k=1,n_core_orb
|
||||||
if (bavard) then
|
kk=list_core(k)
|
||||||
integer :: i
|
do p=1,mo_num
|
||||||
write(6,*)
|
do q=1,mo_num
|
||||||
write(6,*) ' the diagonal of the inactive effective Fock matrix '
|
Fipq(p,q)+=2.D0*bielec_pqxx(p,q,k,k) -bielec_pxxq(p,k,k,q)
|
||||||
write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
|
end do
|
||||||
write(6,*)
|
end do
|
||||||
end if
|
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
|
|
||||||
BEGIN_DOC
|
|
||||||
! the active active Fock matrix, in molecular orbitals
|
|
||||||
! we create them in MOs, quite expensive
|
|
||||||
!
|
|
||||||
! for an implementation in AOs we need first the natural orbitals
|
|
||||||
! for forming an active density matrix in AOs
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
implicit none
|
|
||||||
integer :: p,q,k,kk,t,tt,u,uu
|
|
||||||
|
|
||||||
Fapq = 0.d0
|
|
||||||
|
|
||||||
! the active Fock matrix, D0tu is diagonal
|
|
||||||
do t=1,n_act_orb
|
|
||||||
tt=list_act(t)
|
|
||||||
do q=1,mo_num
|
|
||||||
do p=1,mo_num
|
|
||||||
Fapq(p,q)+=occnum(tt) &
|
|
||||||
*(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q))
|
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
end do
|
! the active Fock matrix, D0tu is diagonal
|
||||||
|
do t=1,n_act_orb
|
||||||
if (bavard) then
|
tt=list_act(t)
|
||||||
integer :: i
|
do p=1,mo_num
|
||||||
write(6,*)
|
do q=1,mo_num
|
||||||
write(6,*) ' the effective Fock matrix over MOs'
|
Fapq(p,q)+=occnum(tt) &
|
||||||
write(6,*)
|
*(bielec_pqxx(p,q,tt,tt)-0.5D0*bielec_pxxq(p,tt,tt,q))
|
||||||
|
end do
|
||||||
write(6,*)
|
end do
|
||||||
write(6,*) ' the diagonal of the inactive effective Fock matrix '
|
end do
|
||||||
write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
|
|
||||||
write(6,*)
|
if (bavard) then
|
||||||
write(6,*)
|
integer :: i
|
||||||
write(6,*) ' the diagonal of the active Fock matrix '
|
write(6,*)
|
||||||
write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num)
|
write(6,*) ' the effective Fock matrix over MOs'
|
||||||
write(6,*)
|
write(6,*)
|
||||||
end if
|
|
||||||
|
write(6,*)
|
||||||
|
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,*)
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) ' the diagonal of the active Fock matrix '
|
||||||
|
write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num)
|
||||||
|
write(6,*)
|
||||||
|
end if
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -1,373 +1,548 @@
|
|||||||
BEGIN_PROVIDER [real*8, occnum, (mo_num)]
|
! -*- F90 -*-
|
||||||
implicit none
|
! diagonalize D0tu
|
||||||
BEGIN_DOC
|
! save the diagonal somewhere, in inverse order
|
||||||
! MO occupation numbers
|
! 4-index-transform the 2-particle density matrix over active orbitals
|
||||||
END_DOC
|
! correct the bielectronic integrals
|
||||||
|
! correct the monoelectronic integrals
|
||||||
integer :: i
|
! put integrals on file, as well orbitals, and the density matrices
|
||||||
occnum=0.D0
|
!
|
||||||
do i=1,n_core_orb
|
subroutine trf_to_natorb
|
||||||
occnum(list_core(i))=2.D0
|
implicit none
|
||||||
end do
|
integer :: i,j,k,l,t,u,p,q,pp
|
||||||
|
real*8 :: eigval(n_act_orb),natorbsCI(n_act_orb,n_act_orb)
|
||||||
do i=1,n_act_orb
|
real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
|
||||||
occnum(list_act(i))=occ_act(n_act_orb-i+1)
|
|
||||||
end do
|
|
||||||
|
|
||||||
write(6,*) ' occupation numbers '
|
call lapack_diag(eigval,natorbsCI,D0tu,n_act_orb,n_act_orb)
|
||||||
do i=1,mo_num
|
write(6,*) ' found occupation numbers as '
|
||||||
write(6,*) i,occnum(i)
|
do i=1,n_act_orb
|
||||||
end do
|
write(6,*) i,eigval(i)
|
||||||
|
end do
|
||||||
|
|
||||||
END_PROVIDER
|
if (bavard) then
|
||||||
|
!
|
||||||
|
|
||||||
BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ]
|
integer :: nmx
|
||||||
&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ]
|
real*8 :: xmx
|
||||||
implicit none
|
do i=1,n_act_orb
|
||||||
BEGIN_DOC
|
! largest element of the eigenvector should be positive
|
||||||
! Natural orbitals of CI
|
xmx=0.D0
|
||||||
END_DOC
|
nmx=0
|
||||||
integer :: i, j
|
do j=1,n_act_orb
|
||||||
|
if (abs(natOrbsCI(j,i)).gt.xmx) then
|
||||||
call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb)
|
nmx=j
|
||||||
|
xmx=abs(natOrbsCI(j,i))
|
||||||
write(6,*) ' found occupation numbers as '
|
end if
|
||||||
do i=1,n_act_orb
|
|
||||||
write(6,*) i,occ_act(i)
|
|
||||||
end do
|
|
||||||
|
|
||||||
if (bavard) then
|
|
||||||
!
|
|
||||||
|
|
||||||
integer :: nmx
|
|
||||||
real*8 :: xmx
|
|
||||||
do i=1,n_act_orb
|
|
||||||
! largest element of the eigenvector should be positive
|
|
||||||
xmx=0.D0
|
|
||||||
nmx=0
|
|
||||||
do j=1,n_act_orb
|
|
||||||
if (abs(natOrbsCI(j,i)).gt.xmx) then
|
|
||||||
nmx=j
|
|
||||||
xmx=abs(natOrbsCI(j,i))
|
|
||||||
end if
|
|
||||||
end do
|
|
||||||
xmx=sign(1.D0,natOrbsCI(nmx,i))
|
|
||||||
do j=1,n_act_orb
|
|
||||||
natOrbsCI(j,i)*=xmx
|
|
||||||
end do
|
|
||||||
|
|
||||||
write(6,*) ' Eigenvector No ',i
|
|
||||||
write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb)
|
|
||||||
end do
|
|
||||||
end if
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! 4-index transformation of 2part matrices
|
|
||||||
END_DOC
|
|
||||||
integer :: i,j,k,l,p,q,pp
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
|
|
||||||
! index per index
|
|
||||||
! first quarter
|
|
||||||
P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:)
|
|
||||||
|
|
||||||
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
|
end do
|
||||||
do p=1,n_act_orb
|
xmx=sign(1.D0,natOrbsCI(nmx,i))
|
||||||
pp=n_act_orb-p+1
|
do j=1,n_act_orb
|
||||||
do q=1,n_act_orb
|
natOrbsCI(j,i)*=xmx
|
||||||
d(pp)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p)
|
end do
|
||||||
|
|
||||||
|
|
||||||
|
write(6,*) ' Eigenvector No ',i
|
||||||
|
write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb)
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
|
||||||
|
do i=1,n_act_orb
|
||||||
|
do j=1,n_act_orb
|
||||||
|
D0tu(i,j)=0.D0
|
||||||
|
end do
|
||||||
|
! fill occupation numbers in descending order
|
||||||
|
D0tu(i,i)=eigval(n_act_orb-i+1)
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
! 4-index transformation of 2part matrices
|
||||||
|
!
|
||||||
|
! index per index
|
||||||
|
! first 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
|
end do
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
do p=1,n_act_orb
|
pp=n_act_orb-p+1
|
||||||
P0tuvx_no(p,j,k,l)=d(p)
|
do q=1,n_act_orb
|
||||||
end do
|
d(pp)+=P0tuvx(q,j,k,l)*natorbsCI(q,p)
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
end do
|
|
||||||
! 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
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p)
|
|
||||||
end do
|
end do
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
do p=1,n_act_orb
|
P0tuvx(p,j,k,l)=d(p)
|
||||||
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
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p)
|
|
||||||
end do
|
end do
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
do p=1,n_act_orb
|
end do
|
||||||
P0tuvx_no(j,k,p,l)=d(p)
|
! 2nd quarter
|
||||||
end do
|
do j=1,n_act_orb
|
||||||
end do
|
do k=1,n_act_orb
|
||||||
end do
|
do l=1,n_act_orb
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
! 4th quarter
|
d(p)=0.D0
|
||||||
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
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p)
|
|
||||||
end do
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=P0tuvx(j,q,k,l)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
P0tuvx(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
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=P0tuvx(j,k,q,l)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
P0tuvx(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
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=P0tuvx(j,k,l,q)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
P0tuvx(j,k,l,p)=d(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' transformed P0tuvx '
|
||||||
|
!
|
||||||
|
! one-electron integrals
|
||||||
|
!
|
||||||
|
do i=1,mo_num
|
||||||
|
do j=1,mo_num
|
||||||
|
onetrf(i,j)=mo_one_e_integrals(i,j)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 1st half-trf
|
||||||
|
do j=1,mo_num
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d(p)=0.D0
|
||||||
end do
|
end do
|
||||||
do p=1,n_act_orb
|
do p=1,n_act_orb
|
||||||
P0tuvx_no(j,k,l,p)=d(p)
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=onetrf(list_act(q),j)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
end do
|
onetrf(list_act(p),j)=d(p)
|
||||||
end do
|
end do
|
||||||
write(6,*) ' transformed P0tuvx '
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, one_ints_no, (mo_num,mo_num)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! Transformed one-e integrals
|
|
||||||
END_DOC
|
|
||||||
integer :: i,j, p, pp, q
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
one_ints_no(:,:)=mo_one_e_integrals(:,:)
|
|
||||||
|
|
||||||
! 1st half-trf
|
|
||||||
do j=1,mo_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=one_ints_no(list_act(q),j)*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
one_ints_no(list_act(p),j)=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
! 2nd half-trf
|
|
||||||
do j=1,mo_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=one_ints_no(j,list_act(q))*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
one_ints_no(j,list_act(p))=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' transformed one_ints '
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! FCI natural orbitals
|
|
||||||
END_DOC
|
|
||||||
integer :: i,j, p, pp, q
|
|
||||||
real*8 :: d(n_act_orb)
|
|
||||||
|
|
||||||
NatOrbsFCI(:,:)=mo_coef(:,:)
|
|
||||||
|
|
||||||
do j=1,ao_num
|
|
||||||
do p=1,n_act_orb
|
|
||||||
d(p)=0.D0
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
pp=n_act_orb-p+1
|
|
||||||
do q=1,n_act_orb
|
|
||||||
d(pp)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
do p=1,n_act_orb
|
|
||||||
NatOrbsFCI(j,list_act(p))=d(p)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
write(6,*) ' transformed orbitals '
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
subroutine trf_to_natorb()
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! save the diagonal somewhere, in inverse order
|
|
||||||
! 4-index-transform the 2-particle density matrix over active orbitals
|
|
||||||
! correct the bielectronic integrals
|
|
||||||
! correct the monoelectronic integrals
|
|
||||||
! put integrals on file, as well orbitals, and the density matrices
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
integer :: i,j,k,l,t,u,p,q,pp
|
|
||||||
real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb)
|
|
||||||
|
|
||||||
! we recalculate total energies
|
|
||||||
write(6,*)
|
|
||||||
write(6,*) ' recalculating energies after the transformation '
|
|
||||||
write(6,*)
|
|
||||||
write(6,*)
|
|
||||||
real*8 :: e_one_all
|
|
||||||
real*8 :: e_two_all
|
|
||||||
integer :: ii
|
|
||||||
integer :: jj
|
|
||||||
integer :: t3
|
|
||||||
integer :: tt
|
|
||||||
integer :: u3
|
|
||||||
integer :: uu
|
|
||||||
integer :: v
|
|
||||||
integer :: v3
|
|
||||||
integer :: vv
|
|
||||||
integer :: x
|
|
||||||
integer :: x3
|
|
||||||
integer :: xx
|
|
||||||
|
|
||||||
e_one_all=0.D0
|
|
||||||
e_two_all=0.D0
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
e_one_all+=2.D0*one_ints_no(ii,ii)
|
|
||||||
do j=1,n_core_orb
|
|
||||||
jj=list_core(j)
|
|
||||||
e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
|
|
||||||
end do
|
|
||||||
do t=1,n_act_orb
|
|
||||||
tt=list_act(t)
|
|
||||||
t3=t+n_core_orb
|
|
||||||
e_two_all += occnum(list_act(t)) * &
|
|
||||||
(2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3))
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
do t=1,n_act_orb
|
|
||||||
tt=list_act(t)
|
|
||||||
e_one_all += occnum(list_act(t))*one_ints_no(tt,tt)
|
|
||||||
do u=1,n_act_orb
|
|
||||||
uu=list_act(u)
|
|
||||||
do v=1,n_act_orb
|
|
||||||
v3=v+n_core_orb
|
|
||||||
do x=1,n_act_orb
|
|
||||||
x3=x+n_core_orb
|
|
||||||
e_two_all +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3)
|
|
||||||
end do
|
end do
|
||||||
end do
|
! 2nd half-trf
|
||||||
end do
|
do j=1,mo_num
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
write(6,*) ' e_one_all = ',e_one_all
|
d(p)=0.D0
|
||||||
write(6,*) ' e_two_all = ',e_two_all
|
|
||||||
ecore =nuclear_repulsion
|
|
||||||
ecore_bis=nuclear_repulsion
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
ecore +=2.D0*one_ints_no(ii,ii)
|
|
||||||
ecore_bis+=2.D0*one_ints_no(ii,ii)
|
|
||||||
do j=1,n_core_orb
|
|
||||||
jj=list_core(j)
|
|
||||||
ecore +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i)
|
|
||||||
ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
eone =0.D0
|
|
||||||
eone_bis=0.D0
|
|
||||||
etwo =0.D0
|
|
||||||
etwo_bis=0.D0
|
|
||||||
etwo_ter=0.D0
|
|
||||||
do t=1,n_act_orb
|
|
||||||
tt=list_act(t)
|
|
||||||
t3=t+n_core_orb
|
|
||||||
eone += occnum(list_act(t))*one_ints_no(tt,tt)
|
|
||||||
eone_bis += occnum(list_act(t))*one_ints_no(tt,tt)
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
eone += occnum(list_act(t)) * &
|
|
||||||
(2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3))
|
|
||||||
eone_bis += occnum(list_act(t)) * &
|
|
||||||
(2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt))
|
|
||||||
end do
|
|
||||||
do u=1,n_act_orb
|
|
||||||
uu=list_act(u)
|
|
||||||
u3=u+n_core_orb
|
|
||||||
do v=1,n_act_orb
|
|
||||||
vv=list_act(v)
|
|
||||||
v3=v+n_core_orb
|
|
||||||
do x=1,n_act_orb
|
|
||||||
xx=list_act(x)
|
|
||||||
x3=x+n_core_orb
|
|
||||||
real*8 :: h1,h2,h3
|
|
||||||
h1=bielec_PQxx_no(tt,uu,v3,x3)
|
|
||||||
h2=bielec_PxxQ_no(tt,u3,v3,xx)
|
|
||||||
h3=bielecCI_no(t,u,v,xx)
|
|
||||||
etwo +=P0tuvx_no(t,u,v,x)*h1
|
|
||||||
etwo_bis+=P0tuvx_no(t,u,v,x)*h2
|
|
||||||
etwo_ter+=P0tuvx_no(t,u,v,x)*h3
|
|
||||||
if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
|
|
||||||
write(6,9901) t,u,v,x,h1,h2,h3
|
|
||||||
9901 format('aie: ',4I4,3E20.12)
|
|
||||||
end if
|
|
||||||
end do
|
end do
|
||||||
end do
|
do p=1,n_act_orb
|
||||||
end do
|
pp=n_act_orb-p+1
|
||||||
end do
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=onetrf(j,list_act(q))*natorbsCI(q,p)
|
||||||
write(6,*) ' energy contributions '
|
end do
|
||||||
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
end do
|
||||||
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
do p=1,n_act_orb
|
||||||
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
onetrf(j,list_act(p))=d(p)
|
||||||
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
end do
|
||||||
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
end do
|
||||||
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
write(6,*) ' transformed onetrf '
|
||||||
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
!
|
||||||
write(6,*) ' ----------------------------------------- '
|
! Orbitals
|
||||||
write(6,*) ' sum of all = ',eone+etwo+ecore
|
!
|
||||||
write(6,*)
|
do j=1,ao_num
|
||||||
SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter
|
do i=1,mo_num
|
||||||
|
NatOrbsFCI(j,i)=mo_coef(j,i)
|
||||||
end subroutine trf_to_natorb
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
do j=1,ao_num
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
NatOrbsFCI(j,list_act(p))=d(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' transformed orbitals '
|
||||||
|
!
|
||||||
|
! now the bielectronic integrals
|
||||||
|
!
|
||||||
|
!!$ write(6,*) ' before the transformation '
|
||||||
|
!!$integer :: kk,ll,ii,jj
|
||||||
|
!!$real*8 :: h1,h2,h3
|
||||||
|
!!$ do i=1,n_act_orb
|
||||||
|
!!$ ii=list_act(i)
|
||||||
|
!!$ do j=1,n_act_orb
|
||||||
|
!!$ jj=list_act(j)
|
||||||
|
!!$ do k=1,n_act_orb
|
||||||
|
!!$ kk=list_act(k)
|
||||||
|
!!$ do l=1,n_act_orb
|
||||||
|
!!$ ll=list_act(l)
|
||||||
|
!!$ h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb)
|
||||||
|
!!$ h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll)
|
||||||
|
!!$ h3=bielecCItmp(i,j,k,ll)
|
||||||
|
!!$ if ((h1.ne.h2).or.(h1.ne.h3)) then
|
||||||
|
!!$ write(6,9901) i,j,k,l,h1,h2,h3
|
||||||
|
!!$9901 format(' aie ',4i4,3E20.12)
|
||||||
|
!!$9902 format('correct',4i4,3E20.12)
|
||||||
|
!!$ else
|
||||||
|
!!$ write(6,9902) i,j,k,l,h1,h2,h3
|
||||||
|
!!$ end if
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
|
||||||
|
do j=1,mo_num
|
||||||
|
do k=1,n_core_orb+n_act_orb
|
||||||
|
do l=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d1(p)=0.D0
|
||||||
|
d2(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d1(pp)+=bielec_PQxxtmp(list_act(q),j,k,l)*natorbsCI(q,p)
|
||||||
|
d2(pp)+=bielec_PxxQtmp(list_act(q),k,l,j)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielec_PQxxtmp(list_act(p),j,k,l)=d1(p)
|
||||||
|
bielec_PxxQtmp(list_act(p),k,l,j)=d2(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 2nd quarter
|
||||||
|
do j=1,mo_num
|
||||||
|
do k=1,n_core_orb+n_act_orb
|
||||||
|
do l=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d1(p)=0.D0
|
||||||
|
d2(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d1(pp)+=bielec_PQxxtmp(j,list_act(q),k,l)*natorbsCI(q,p)
|
||||||
|
d2(pp)+=bielec_PxxQtmp(j,k,l,list_act(q))*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielec_PQxxtmp(j,list_act(p),k,l)=d1(p)
|
||||||
|
bielec_PxxQtmp(j,k,l,list_act(p))=d2(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 3rd quarter
|
||||||
|
do j=1,mo_num
|
||||||
|
do k=1,mo_num
|
||||||
|
do l=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d1(p)=0.D0
|
||||||
|
d2(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d1(pp)+=bielec_PQxxtmp(j,k,n_core_orb+q,l)*natorbsCI(q,p)
|
||||||
|
d2(pp)+=bielec_PxxQtmp(j,n_core_orb+q,l,k)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielec_PQxxtmp(j,k,n_core_orb+p,l)=d1(p)
|
||||||
|
bielec_PxxQtmp(j,n_core_orb+p,l,k)=d2(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 4th quarter
|
||||||
|
do j=1,mo_num
|
||||||
|
do k=1,mo_num
|
||||||
|
do l=1,n_core_orb+n_act_orb
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d1(p)=0.D0
|
||||||
|
d2(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d1(pp)+=bielec_PQxxtmp(j,k,l,n_core_orb+q)*natorbsCI(q,p)
|
||||||
|
d2(pp)+=bielec_PxxQtmp(j,l,n_core_orb+q,k)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielec_PQxxtmp(j,k,l,n_core_orb+p)=d1(p)
|
||||||
|
bielec_PxxQtmp(j,l,n_core_orb+p,k)=d2(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' transformed PQxx and PxxQ '
|
||||||
|
!
|
||||||
|
! and finally the bielecCI integrals
|
||||||
|
!
|
||||||
|
do j=1,n_act_orb
|
||||||
|
do k=1,n_act_orb
|
||||||
|
do l=1,mo_num
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=bielecCItmp(q,j,k,l)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielecCItmp(p,j,k,l)=d(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
! 2nd quarter
|
||||||
|
do j=1,n_act_orb
|
||||||
|
do k=1,n_act_orb
|
||||||
|
do l=1,mo_num
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=bielecCItmp(j,q,k,l)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielecCItmp(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,mo_num
|
||||||
|
do p=1,n_act_orb
|
||||||
|
d(p)=0.D0
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=bielecCItmp(j,k,q,l)*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielecCItmp(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
|
||||||
|
pp=n_act_orb-p+1
|
||||||
|
do q=1,n_act_orb
|
||||||
|
d(pp)+=bielecCItmp(j,k,l,list_act(q))*natorbsCI(q,p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do p=1,n_act_orb
|
||||||
|
bielecCItmp(j,k,l,list_act(p))=d(p)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' transformed tuvP '
|
||||||
|
!
|
||||||
|
! that's all
|
||||||
|
!
|
||||||
|
!!$
|
||||||
|
!!$! test coherence of the bielectronic integals
|
||||||
|
!!$! PQxx = PxxQ = tuvP for some of the indices
|
||||||
|
!!$ write(6,*) ' after the transformation '
|
||||||
|
!!$ do i=1,n_act_orb
|
||||||
|
!!$ ii=list_act(i)
|
||||||
|
!!$ do j=1,n_act_orb
|
||||||
|
!!$ jj=list_act(j)
|
||||||
|
!!$ do k=1,n_act_orb
|
||||||
|
!!$ kk=list_act(k)
|
||||||
|
!!$ do l=1,n_act_orb
|
||||||
|
!!$ ll=list_act(l)
|
||||||
|
!!$ h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb)
|
||||||
|
!!$ h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll)
|
||||||
|
!!$ h3=bielecCItmp(i,j,k,ll)
|
||||||
|
!!$ if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
|
||||||
|
!!$ write(6,9901) i,j,k,l,h1,h1-h2,h1-h3
|
||||||
|
!!$ else
|
||||||
|
!!$ write(6,9902) i,j,k,l,h1,h2,h3
|
||||||
|
!!$ end if
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
!!$ end do
|
||||||
|
|
||||||
|
! we recalculate total energies
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) ' recalculating energies after the transformation '
|
||||||
|
write(6,*)
|
||||||
|
write(6,*)
|
||||||
|
real*8 :: e_one_all
|
||||||
|
real*8 :: e_two_all
|
||||||
|
integer :: ii
|
||||||
|
integer :: jj
|
||||||
|
integer :: t3
|
||||||
|
integer :: tt
|
||||||
|
integer :: u3
|
||||||
|
integer :: uu
|
||||||
|
integer :: v
|
||||||
|
integer :: v3
|
||||||
|
integer :: vv
|
||||||
|
integer :: x
|
||||||
|
integer :: x3
|
||||||
|
integer :: xx
|
||||||
|
|
||||||
|
e_one_all=0.D0
|
||||||
|
e_two_all=0.D0
|
||||||
|
do i=1,n_core_orb
|
||||||
|
ii=list_core(i)
|
||||||
|
e_one_all+=2.D0*onetrf(ii,ii)
|
||||||
|
do j=1,n_core_orb
|
||||||
|
jj=list_core(j)
|
||||||
|
e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
|
||||||
|
end do
|
||||||
|
do t=1,n_act_orb
|
||||||
|
tt=list_act(t)
|
||||||
|
t3=t+n_core_orb
|
||||||
|
do u=1,n_act_orb
|
||||||
|
uu=list_act(u)
|
||||||
|
u3=u+n_core_orb
|
||||||
|
e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
|
||||||
|
-bielec_PQxxtmp(tt,ii,i,u3))
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do t=1,n_act_orb
|
||||||
|
tt=list_act(t)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
uu=list_act(u)
|
||||||
|
e_one_all+=D0tu(t,u)*onetrf(tt,uu)
|
||||||
|
do v=1,n_act_orb
|
||||||
|
v3=v+n_core_orb
|
||||||
|
do x=1,n_act_orb
|
||||||
|
x3=x+n_core_orb
|
||||||
|
e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' e_one_all = ',e_one_all
|
||||||
|
write(6,*) ' e_two_all = ',e_two_all
|
||||||
|
ecore =nuclear_repulsion
|
||||||
|
ecore_bis=nuclear_repulsion
|
||||||
|
do i=1,n_core_orb
|
||||||
|
ii=list_core(i)
|
||||||
|
ecore +=2.D0*onetrf(ii,ii)
|
||||||
|
ecore_bis+=2.D0*onetrf(ii,ii)
|
||||||
|
do j=1,n_core_orb
|
||||||
|
jj=list_core(j)
|
||||||
|
ecore +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
|
||||||
|
ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
eone =0.D0
|
||||||
|
eone_bis=0.D0
|
||||||
|
etwo =0.D0
|
||||||
|
etwo_bis=0.D0
|
||||||
|
etwo_ter=0.D0
|
||||||
|
do t=1,n_act_orb
|
||||||
|
tt=list_act(t)
|
||||||
|
t3=t+n_core_orb
|
||||||
|
do u=1,n_act_orb
|
||||||
|
uu=list_act(u)
|
||||||
|
u3=u+n_core_orb
|
||||||
|
eone +=D0tu(t,u)*onetrf(tt,uu)
|
||||||
|
eone_bis+=D0tu(t,u)*onetrf(tt,uu)
|
||||||
|
do i=1,n_core_orb
|
||||||
|
ii=list_core(i)
|
||||||
|
eone +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
|
||||||
|
-bielec_PQxxtmp(tt,ii,i,u3))
|
||||||
|
eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) &
|
||||||
|
-bielec_PxxQtmp(tt,i,i,uu))
|
||||||
|
end do
|
||||||
|
do v=1,n_act_orb
|
||||||
|
vv=list_act(v)
|
||||||
|
v3=v+n_core_orb
|
||||||
|
do x=1,n_act_orb
|
||||||
|
xx=list_act(x)
|
||||||
|
x3=x+n_core_orb
|
||||||
|
real*8 :: h1,h2,h3
|
||||||
|
h1=bielec_PQxxtmp(tt,uu,v3,x3)
|
||||||
|
h2=bielec_PxxQtmp(tt,u3,v3,xx)
|
||||||
|
h3=bielecCItmp(t,u,v,xx)
|
||||||
|
etwo +=P0tuvx(t,u,v,x)*h1
|
||||||
|
etwo_bis+=P0tuvx(t,u,v,x)*h2
|
||||||
|
etwo_ter+=P0tuvx(t,u,v,x)*h3
|
||||||
|
if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then
|
||||||
|
write(6,9901) t,u,v,x,h1,h2,h3
|
||||||
|
9901 format('aie: ',4I4,3E20.12)
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
write(6,*) ' energy contributions '
|
||||||
|
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
||||||
|
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
||||||
|
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
||||||
|
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
||||||
|
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
||||||
|
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
||||||
|
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
||||||
|
write(6,*) ' ----------------------------------------- '
|
||||||
|
write(6,*) ' sum of all = ',eone+etwo+ecore
|
||||||
|
write(6,*)
|
||||||
|
|
||||||
|
end subroutine trf_to_natorb
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [real*8, onetrf, (mo_num,mo_num)]
|
||||||
|
&BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)]
|
||||||
|
END_PROVIDER
|
||||||
|
|||||||
65
src/casscf/natorb_casscf.irp.f
Normal file
65
src/casscf/natorb_casscf.irp.f
Normal file
@ -0,0 +1,65 @@
|
|||||||
|
! -*- F90 -*-
|
||||||
|
BEGIN_PROVIDER [real*8, occnum, (mo_num)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,kk,j
|
||||||
|
logical :: lread
|
||||||
|
real*8 :: rdum
|
||||||
|
do i=1,mo_num
|
||||||
|
occnum(i)=0.D0
|
||||||
|
end do
|
||||||
|
do i=1,n_core_orb
|
||||||
|
occnum(list_core(i))=2.D0
|
||||||
|
end do
|
||||||
|
|
||||||
|
open(unit=12,file='D0tu.dat',form='formatted',status='old')
|
||||||
|
lread=.true.
|
||||||
|
do while (lread)
|
||||||
|
read(12,*,iostat=kk) i,j,rdum
|
||||||
|
if (kk.ne.0) then
|
||||||
|
lread=.false.
|
||||||
|
else
|
||||||
|
if (i.eq.j) then
|
||||||
|
occnum(list_act(i))=rdum
|
||||||
|
else
|
||||||
|
write(6,*) ' WARNING - no natural orbitals !'
|
||||||
|
write(6,*) i,j,rdum
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
close(12)
|
||||||
|
write(6,*) ' read occupation numbers '
|
||||||
|
do i=1,mo_num
|
||||||
|
write(6,*) i,occnum(i)
|
||||||
|
end do
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,k,l,kk
|
||||||
|
real*8 :: rdum
|
||||||
|
logical :: lread
|
||||||
|
|
||||||
|
do i=1,n_act_orb
|
||||||
|
do j=1,n_act_orb
|
||||||
|
do k=1,n_act_orb
|
||||||
|
do l=1,n_act_orb
|
||||||
|
P0tuvx_no(l,k,j,i)=0.D0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
|
||||||
|
open(unit=12,file='P0tuvx.dat',form='formatted',status='old')
|
||||||
|
lread=.true.
|
||||||
|
do while (lread)
|
||||||
|
read(12,*,iostat=kk) i,j,k,l,rdum
|
||||||
|
if (kk.ne.0) then
|
||||||
|
lread=.false.
|
||||||
|
else
|
||||||
|
P0tuvx_no(i,j,k,l)=rdum
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
close(12)
|
||||||
|
write(6,*) ' read the 2-particle density matrix '
|
||||||
|
END_PROVIDER
|
||||||
26
src/casscf/one_ints.irp.f
Normal file
26
src/casscf/one_ints.irp.f
Normal file
@ -0,0 +1,26 @@
|
|||||||
|
! -*- F90 -*-
|
||||||
|
BEGIN_PROVIDER [real*8, one_ints, (mo_num,mo_num)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,j,kk
|
||||||
|
logical :: lread
|
||||||
|
real*8 :: rdum
|
||||||
|
do i=1,mo_num
|
||||||
|
do j=1,mo_num
|
||||||
|
one_ints(i,j)=0.D0
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
open(unit=12,file='onetrf.tmp',status='old',form='formatted')
|
||||||
|
lread=.true.
|
||||||
|
do while (lread)
|
||||||
|
read(12,*,iostat=kk) i,j,rdum
|
||||||
|
if (kk.ne.0) then
|
||||||
|
lread=.false.
|
||||||
|
else
|
||||||
|
one_ints(i,j)=rdum
|
||||||
|
one_ints(j,i)=rdum
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
close(12)
|
||||||
|
write(6,*) ' read MCSCF natural one-electron integrals '
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
@ -1,3 +1,4 @@
|
|||||||
|
! -*- F90 -*-
|
||||||
BEGIN_PROVIDER [real*8, etwo]
|
BEGIN_PROVIDER [real*8, etwo]
|
||||||
&BEGIN_PROVIDER [real*8, eone]
|
&BEGIN_PROVIDER [real*8, eone]
|
||||||
&BEGIN_PROVIDER [real*8, eone_bis]
|
&BEGIN_PROVIDER [real*8, eone_bis]
|
||||||
@ -5,117 +6,117 @@
|
|||||||
&BEGIN_PROVIDER [real*8, etwo_ter]
|
&BEGIN_PROVIDER [real*8, etwo_ter]
|
||||||
&BEGIN_PROVIDER [real*8, ecore]
|
&BEGIN_PROVIDER [real*8, ecore]
|
||||||
&BEGIN_PROVIDER [real*8, ecore_bis]
|
&BEGIN_PROVIDER [real*8, ecore_bis]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3
|
integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3
|
||||||
real*8 :: e_one_all,e_two_all
|
real*8 :: e_one_all,e_two_all
|
||||||
e_one_all=0.D0
|
e_one_all=0.D0
|
||||||
e_two_all=0.D0
|
e_two_all=0.D0
|
||||||
do i=1,n_core_orb
|
do i=1,n_core_orb
|
||||||
ii=list_core(i)
|
ii=list_core(i)
|
||||||
e_one_all+=2.D0*mo_one_e_integrals(ii,ii)
|
e_one_all+=2.D0*mo_one_e_integrals(ii,ii)
|
||||||
do j=1,n_core_orb
|
do j=1,n_core_orb
|
||||||
jj=list_core(j)
|
jj=list_core(j)
|
||||||
e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
|
e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
|
||||||
end do
|
end do
|
||||||
do t=1,n_act_orb
|
do t=1,n_act_orb
|
||||||
tt=list_act(t)
|
tt=list_act(t)
|
||||||
t3=t+n_core_orb
|
t3=t+n_core_orb
|
||||||
do u=1,n_act_orb
|
do u=1,n_act_orb
|
||||||
|
uu=list_act(u)
|
||||||
|
u3=u+n_core_orb
|
||||||
|
e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
|
||||||
|
-bielec_PQxxtmp(tt,ii,i,u3))
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
do t=1,n_act_orb
|
||||||
|
tt=list_act(t)
|
||||||
|
do u=1,n_act_orb
|
||||||
|
uu=list_act(u)
|
||||||
|
e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
||||||
|
do v=1,n_act_orb
|
||||||
|
v3=v+n_core_orb
|
||||||
|
do x=1,n_act_orb
|
||||||
|
x3=x+n_core_orb
|
||||||
|
e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(6,*) ' e_one_all = ',e_one_all
|
||||||
|
write(6,*) ' e_two_all = ',e_two_all
|
||||||
|
ecore =nuclear_repulsion
|
||||||
|
ecore_bis=nuclear_repulsion
|
||||||
|
do i=1,n_core_orb
|
||||||
|
ii=list_core(i)
|
||||||
|
ecore +=2.D0*mo_one_e_integrals(ii,ii)
|
||||||
|
ecore_bis+=2.D0*mo_one_e_integrals(ii,ii)
|
||||||
|
do j=1,n_core_orb
|
||||||
|
jj=list_core(j)
|
||||||
|
ecore +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i)
|
||||||
|
ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii)
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
eone =0.D0
|
||||||
|
eone_bis=0.D0
|
||||||
|
etwo =0.D0
|
||||||
|
etwo_bis=0.D0
|
||||||
|
etwo_ter=0.D0
|
||||||
|
do t=1,n_act_orb
|
||||||
|
tt=list_act(t)
|
||||||
|
t3=t+n_core_orb
|
||||||
|
do u=1,n_act_orb
|
||||||
uu=list_act(u)
|
uu=list_act(u)
|
||||||
u3=u+n_core_orb
|
u3=u+n_core_orb
|
||||||
e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) &
|
eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
||||||
-bielec_PQxx(tt,ii,i,u3))
|
eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
||||||
end do
|
do i=1,n_core_orb
|
||||||
end do
|
ii=list_core(i)
|
||||||
end do
|
eone +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) &
|
||||||
do t=1,n_act_orb
|
-bielec_PQxxtmp(tt,ii,i,u3))
|
||||||
tt=list_act(t)
|
eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) &
|
||||||
do u=1,n_act_orb
|
-bielec_PxxQtmp(tt,i,i,uu))
|
||||||
uu=list_act(u)
|
|
||||||
e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
|
||||||
do v=1,n_act_orb
|
|
||||||
v3=v+n_core_orb
|
|
||||||
do x=1,n_act_orb
|
|
||||||
x3=x+n_core_orb
|
|
||||||
e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3)
|
|
||||||
end do
|
end do
|
||||||
end do
|
do v=1,n_act_orb
|
||||||
end do
|
vv=list_act(v)
|
||||||
end do
|
v3=v+n_core_orb
|
||||||
write(6,*) ' e_one_all = ',e_one_all
|
do x=1,n_act_orb
|
||||||
write(6,*) ' e_two_all = ',e_two_all
|
|
||||||
ecore =nuclear_repulsion
|
|
||||||
ecore_bis=nuclear_repulsion
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
ecore +=2.D0*mo_one_e_integrals(ii,ii)
|
|
||||||
ecore_bis+=2.D0*mo_one_e_integrals(ii,ii)
|
|
||||||
do j=1,n_core_orb
|
|
||||||
jj=list_core(j)
|
|
||||||
ecore +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i)
|
|
||||||
ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii)
|
|
||||||
end do
|
|
||||||
end do
|
|
||||||
eone =0.D0
|
|
||||||
eone_bis=0.D0
|
|
||||||
etwo =0.D0
|
|
||||||
etwo_bis=0.D0
|
|
||||||
etwo_ter=0.D0
|
|
||||||
do t=1,n_act_orb
|
|
||||||
tt=list_act(t)
|
|
||||||
t3=t+n_core_orb
|
|
||||||
do u=1,n_act_orb
|
|
||||||
uu=list_act(u)
|
|
||||||
u3=u+n_core_orb
|
|
||||||
eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
|
||||||
eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu)
|
|
||||||
do i=1,n_core_orb
|
|
||||||
ii=list_core(i)
|
|
||||||
eone +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) &
|
|
||||||
-bielec_PQxx(tt,ii,i,u3))
|
|
||||||
eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii) &
|
|
||||||
-bielec_PxxQ(tt,i,i,uu))
|
|
||||||
end do
|
|
||||||
do v=1,n_act_orb
|
|
||||||
vv=list_act(v)
|
|
||||||
v3=v+n_core_orb
|
|
||||||
do x=1,n_act_orb
|
|
||||||
xx=list_act(x)
|
xx=list_act(x)
|
||||||
x3=x+n_core_orb
|
x3=x+n_core_orb
|
||||||
real*8 :: h1,h2,h3
|
real*8 :: h1,h2,h3
|
||||||
h1=bielec_PQxx(tt,uu,v3,x3)
|
h1=bielec_PQxxtmp(tt,uu,v3,x3)
|
||||||
h2=bielec_PxxQ(tt,u3,v3,xx)
|
h2=bielec_PxxQtmp(tt,u3,v3,xx)
|
||||||
h3=bielecCI(t,u,v,xx)
|
h3=bielecCItmp(t,u,v,xx)
|
||||||
etwo +=P0tuvx(t,u,v,x)*h1
|
etwo +=P0tuvx(t,u,v,x)*h1
|
||||||
etwo_bis+=P0tuvx(t,u,v,x)*h2
|
etwo_bis+=P0tuvx(t,u,v,x)*h2
|
||||||
etwo_ter+=P0tuvx(t,u,v,x)*h3
|
etwo_ter+=P0tuvx(t,u,v,x)*h3
|
||||||
if ((h1.ne.h2).or.(h1.ne.h3)) then
|
if ((h1.ne.h2).or.(h1.ne.h3)) then
|
||||||
write(6,9901) t,u,v,x,h1,h2,h3
|
write(6,9901) t,u,v,x,h1,h2,h3
|
||||||
9901 format('aie: ',4I4,3E20.12)
|
9901 format('aie: ',4I4,3E20.12)
|
||||||
end if
|
end if
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
|
end do
|
||||||
end do
|
end do
|
||||||
end do
|
|
||||||
end do
|
write(6,*) ' energy contributions '
|
||||||
|
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
||||||
write(6,*) ' energy contributions '
|
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
||||||
write(6,*) ' core energy = ',ecore,' using PQxx integrals '
|
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
||||||
write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals '
|
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
||||||
write(6,*) ' 1el energy = ',eone ,' using PQxx integrals '
|
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
||||||
write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals '
|
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
||||||
write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals '
|
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
||||||
write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals '
|
write(6,*) ' ----------------------------------------- '
|
||||||
write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals '
|
write(6,*) ' sum of all = ',eone+etwo+ecore
|
||||||
write(6,*) ' ----------------------------------------- '
|
write(6,*)
|
||||||
write(6,*) ' sum of all = ',eone+etwo+ecore
|
write(6,*) ' nuclear (qp) = ',nuclear_repulsion
|
||||||
write(6,*)
|
write(6,*) ' core energy (qp) = ',core_energy
|
||||||
write(6,*) ' nuclear (qp) = ',nuclear_repulsion
|
write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1)
|
||||||
write(6,*) ' core energy (qp) = ',core_energy
|
write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1)
|
||||||
write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1)
|
write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1)
|
||||||
write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1)
|
write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1)
|
||||||
write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1)
|
|
||||||
write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1)
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user