2022-02-04 20:11:34 +01:00
|
|
|
program determinants
|
2022-02-06 20:22:32 +01:00
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
|
integer*16 :: o, v
|
|
|
|
|
integer*16 :: e, k
|
|
|
|
|
integer*16 :: n, m
|
|
|
|
|
|
|
|
|
|
integer*16 :: s0d0
|
|
|
|
|
integer*16 :: s2d1
|
|
|
|
|
integer*16 :: s0d2, s2d2, s4d2
|
|
|
|
|
integer*16 :: s2d3, s4d3, s6d3
|
|
|
|
|
integer*16 :: s0d4, s2d4, s4d4, s6d4, s8d4
|
|
|
|
|
integer*16 :: s2d5, s4d5, s6d5
|
|
|
|
|
integer*16 :: s0d6, s2d6, s4d6
|
|
|
|
|
integer*16 :: s2d7
|
|
|
|
|
integer*16 :: s0d8
|
|
|
|
|
|
|
|
|
|
integer*16 :: rhf, fci
|
|
|
|
|
integer*16 :: cis0, cis2, cis4, cis6
|
|
|
|
|
integer*16 :: cis, cisd, cisdt, cisdtq
|
|
|
|
|
integer*16 :: cio1, cio1p5, cio2, cio2p5, cio3, cio3p5, cio4
|
|
|
|
|
|
|
|
|
|
integer :: nunit
|
|
|
|
|
character(len=8), parameter :: fmt1 = "(*(i16))"
|
|
|
|
|
character(len=8), parameter :: fmt2 = "(a,i16)"
|
|
|
|
|
|
|
|
|
|
character(len=32) :: filename
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
read(*,*) e, k
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
o = e/2
|
|
|
|
|
v = k - o
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
call binomial_coefficient(k,o,m)
|
|
|
|
|
! CIs0
|
|
|
|
|
cis0 = m
|
|
|
|
|
! CIs2
|
|
|
|
|
cis2 = cis0 + m*o*v
|
|
|
|
|
! CIs4
|
|
|
|
|
cis4 = cis2 + m*o*(o-1)*v*(v-1)/4
|
|
|
|
|
! CIs6
|
|
|
|
|
cis6 = cis4 + m*o*(o-1)*(o-2)*v*(v-1)*(v-2)/36
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
filename = 'det_CIs.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt2) '0 ', cis0
|
|
|
|
|
write(nunit,fmt2) '2 ', cis2
|
|
|
|
|
write(nunit,fmt2) '4 ', cis4
|
|
|
|
|
write(nunit,fmt2) '6 ', cis6
|
|
|
|
|
close(nunit)
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
s0d0 = 1
|
|
|
|
|
|
|
|
|
|
s0d2 = o*v
|
|
|
|
|
|
|
|
|
|
s0d4 = o*(o-1)*v*(v-1)/4
|
|
|
|
|
|
|
|
|
|
s0d6 = o*(o-1)*(o-2)*v*(v-1)*(v-2)/36
|
|
|
|
|
|
|
|
|
|
s0d8 = o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2)*(v-3)/400
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
s2d1 = 2*o*v
|
|
|
|
|
|
|
|
|
|
!s2d2 = o *v*(v-1) + &
|
|
|
|
|
! o*(o-1)*v
|
|
|
|
|
s2d2 = o*v*(o+v-2)
|
|
|
|
|
|
|
|
|
|
s2d3 = 2*o*(o-1)*v*(v-1)
|
|
|
|
|
|
|
|
|
|
!s2d4 = o*(o-1) *v*(v-1)*(v-2)/2 + &
|
|
|
|
|
! o*(o-1)*(o-2)*v*(v-1) /2
|
|
|
|
|
s2d4 = o*(o-1)*v*(v-1)*(o+v-4)/2
|
|
|
|
|
|
|
|
|
|
s2d5 = o*(o-1)*(o-2)*v*(v-1)*(v-2)/2
|
|
|
|
|
|
|
|
|
|
!s2d6 = o*(o-1)*(o-2) *v*(v-1)*(v-2)*(v-3)/12 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2) /12
|
|
|
|
|
s2d6 = o*(o-1)*(o-2)*v*(v-1)*(v-2)*(o+v-6)/12
|
|
|
|
|
|
|
|
|
|
s2d7 = o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2)*(v-3)/18
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
s4d2 = o*(o-1)*v*(v-1)*3/2
|
|
|
|
|
|
|
|
|
|
!s4d3 = o*(o-1) *v*(v-1)*(v-2) + &
|
|
|
|
|
! o*(o-1)*(o-2)*v*(v-1)
|
|
|
|
|
s4d3 = o*(o-1)*v*(v-1)*(o+v-4)
|
|
|
|
|
|
|
|
|
|
!s4d4 = o*(o-1) *v*(v-1)*(v-2)*(v-3)/8 + &
|
|
|
|
|
! o*(o-1)*(o-2) *v*(v-1)*(v-2) *3/2 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3)*v*(v-1) /8
|
|
|
|
|
s4d4 = o*(o-1)*v*(v-1) * ( (v-2)*(v-3) + (o-2)*(o-3) + 12*(o-2)*(v-2) )/8
|
|
|
|
|
|
|
|
|
|
!s4d5 = o*(o-1)*(o-2) *v*(v-1)*(v-2)*(v-3)/2 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2) /2
|
|
|
|
|
s4d5 = o*(o-1)*(o-2)*v*(v-1)*(v-2)*(o+v-6)/2
|
|
|
|
|
|
|
|
|
|
!s4d6 = o*(o-1)*(o-2) *v*(v-1)*(v-2)*(v-3)*(v-4)/24 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3) *v*(v-1)*(v-2)*(v-3) *3/8 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3)*(o-4)*v*(v-1)*(v-2) /24
|
|
|
|
|
s4d6 = o*(o-1)*(o-2)*v*(v-1)*(v-2) * ( (o-3)*(o-4) + (v-3)*(v-4) + 3*(o-3)*(v-3) )/24
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
s6d3 = o*(o-1)*(o-2)*v*(v-1)*(v-2)*5/9
|
|
|
|
|
|
|
|
|
|
!s6d4 = o*(o-1)*(o-2) *v*(v-1)*(v-2)*(v-3)*5/12 + &
|
|
|
|
|
! o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2) *5/12
|
|
|
|
|
s6d4 = o*(o-1)*(o-2)*v*(v-1)*(v-2)*(o+v-6)*5/12
|
|
|
|
|
|
|
|
|
|
s6d5 = o*(o-1)*(o-2) *v*(v-1)*(v-2)*(v-3)*(v-4)/12 + &
|
|
|
|
|
o*(o-1)*(o-2)*(o-3) *v*(v-1)*(v-2)*(v-3) *5/9 + &
|
|
|
|
|
o*(o-1)*(o-2)*(o-3)*(o-4)*v*(v-1)*(v-2) /12
|
|
|
|
|
!s6d5 = o*(o-1)*(o-2)*v*(v-1)*(v-2) * ( 3*(o-3)*(o-4) + 3*(v-3)*(v-4) + 20*(o-3)*(v-3) )/36
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
s8d4 = o*(o-1)*(o-2)*(o-3)*v*(v-1)*(v-2)*(v-3)*35/288
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
! HF
|
|
|
|
|
rhf = s0d0
|
|
|
|
|
! CIo1
|
|
|
|
|
cio1 = rhf + s0d2 + s2d1
|
|
|
|
|
! CIo1.5
|
|
|
|
|
cio1p5 = cio1 + s2d2
|
|
|
|
|
! CIo2
|
|
|
|
|
cio2 = cio1p5 + s0d4 + s2d3 + s4d2
|
|
|
|
|
! CIo2.5
|
|
|
|
|
cio2p5 = cio2 + s2d4 + s4d3
|
|
|
|
|
! CIo3
|
|
|
|
|
cio3 = cio2p5 + s0d6 + s2d5 + s4d4 + s6d3
|
|
|
|
|
! CIo3.5
|
|
|
|
|
cio3p5 = cio3 + s2d6 + s4d5 + s6d4
|
|
|
|
|
! CIo4
|
|
|
|
|
cio4 = cio3p5 + s0d8 + s2d7 + s4d6 + s6d5 + s8d4
|
|
|
|
|
|
|
|
|
|
filename = 'det_CIo.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt2) '0 ', rhf
|
|
|
|
|
write(nunit,fmt2) '1 ', cio1
|
|
|
|
|
write(nunit,fmt2) '1.5 ', cio1p5
|
|
|
|
|
write(nunit,fmt2) '2 ', cio2
|
|
|
|
|
write(nunit,fmt2) '2.5 ', cio2p5
|
|
|
|
|
write(nunit,fmt2) '3 ', cio3
|
|
|
|
|
write(nunit,fmt2) '3.5 ', cio3p5
|
|
|
|
|
write(nunit,fmt2) '4 ', cio4
|
|
|
|
|
close(nunit)
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
! HF
|
|
|
|
|
rhf = s0d0
|
|
|
|
|
! CIS
|
|
|
|
|
cis = rhf + s2d1
|
|
|
|
|
! CISD
|
|
|
|
|
cisd = cis + s0d2 + s2d2 + s4d2
|
|
|
|
|
! CISDT
|
|
|
|
|
cisdt = cisd + s2d3 + s4d3 + s6d3
|
|
|
|
|
! CISDTQ
|
|
|
|
|
cisdtq = cisdt + s0d4 + s2d4 + s4d4 + s8d4
|
|
|
|
|
|
|
|
|
|
filename = 'det_CI.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt2) '0 ', rhf
|
|
|
|
|
write(nunit,fmt2) '1 ', cis
|
|
|
|
|
write(nunit,fmt2) '2 ', cisd
|
|
|
|
|
write(nunit,fmt2) '3 ', cisdt
|
|
|
|
|
write(nunit,fmt2) '4 ', cisdtq
|
|
|
|
|
close(nunit)
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
call binomial_coefficient(k,o,m)
|
|
|
|
|
|
|
|
|
|
fci = m*m
|
|
|
|
|
|
|
|
|
|
filename = 'det_FCI.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt2) '0 ', fci
|
|
|
|
|
close(nunit)
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
filename = 'det.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt2) 'FCI ', fci
|
|
|
|
|
write(nunit,fmt2) 'RHF ', rhf
|
|
|
|
|
write(nunit,fmt2) 'CIS ', cis
|
|
|
|
|
write(nunit,fmt2) 'CISD ', cisd
|
|
|
|
|
write(nunit,fmt2) 'CISDT ', cisdt
|
|
|
|
|
write(nunit,fmt2) 'CISDTQ ', cisdtq
|
|
|
|
|
write(nunit,fmt2) 'CIs0 ', cis0
|
|
|
|
|
write(nunit,fmt2) 'CIs2 ', cis2
|
|
|
|
|
write(nunit,fmt2) 'CIs4 ', cis4
|
|
|
|
|
write(nunit,fmt2) 'CIo1 ', cio1
|
|
|
|
|
write(nunit,fmt2) 'CIo1.5 ', cio1p5
|
|
|
|
|
write(nunit,fmt2) 'CIo2 ', cio2
|
|
|
|
|
write(nunit,fmt2) 'CIo2.5 ', cio2p5
|
|
|
|
|
write(nunit,fmt2) 'CIo3 ', cio3
|
|
|
|
|
write(nunit,fmt2) 'CIo3.5 ', cio3p5
|
|
|
|
|
write(nunit,fmt2) 'CIo4 ', cio4
|
|
|
|
|
close(nunit)
|
|
|
|
|
|
|
|
|
|
!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
|
|
|
|
|
|
|
filename = 'table.dat'
|
|
|
|
|
open(newunit=nunit,file=filename,status='unknown')
|
|
|
|
|
write(nunit,fmt1) s0d0
|
|
|
|
|
write(nunit,fmt1) 0, s2d1
|
|
|
|
|
write(nunit,fmt1) s0d2, s2d2, s4d2
|
|
|
|
|
write(nunit,fmt1) 0, s2d3, s4d3, s6d3
|
|
|
|
|
write(nunit,fmt1) s0d4, s2d4, s4d4, s6d4, s8d4
|
|
|
|
|
write(nunit,fmt1) 0, s2d5, s4d5, s6d5
|
|
|
|
|
write(nunit,fmt1) s0d6, s2d6, s4d6
|
|
|
|
|
write(nunit,fmt1) 0, s2d7
|
|
|
|
|
write(nunit,fmt1) s0d8
|
|
|
|
|
close(nunit)
|
2022-02-04 20:11:34 +01:00
|
|
|
|
|
|
|
|
end program determinants
|
|
|
|
|
|
2022-02-06 20:22:32 +01:00
|
|
|
subroutine binomial_coefficient(n,k,c)
|
|
|
|
|
implicit none
|
|
|
|
|
! (n k) = n!/((n-k)!k!)
|
|
|
|
|
integer*16, intent(in) :: n, k
|
|
|
|
|
integer*16, intent(out) :: c
|
|
|
|
|
integer :: i
|
|
|
|
|
real*8 :: cr
|
|
|
|
|
cr = 1.0d0
|
|
|
|
|
do i=1,k
|
|
|
|
|
cr = cr * real( n+1-i) / real( i )
|
|
|
|
|
end do
|
|
|
|
|
c = int( cr )
|
|
|
|
|
end subroutine binomial_coefficient
|
