PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-23 11:17:38 +02:00
Code Releases Activity

star implementation of CISD

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-20 12:28:19 +02:00
parent 88c13a8ec9
commit 13b170fc2f
10 changed files with 292 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
input/basis
View File

@ -1,39 +1,9 @@
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
1 3
S 3
1 38.3600000 0.0238090
2 5.7700000 0.1548910
3 1.2400000 0.4699870
S 1
1 2.8360000 1.0000000
S 1
1 0.3782000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
1 0.2976000 1.0000000
P 1
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000
1 1.2750000 1.0000000

10
input/methods
View File

@ -6,13 +6,15 @@
F F F
# drCCD rCCD lCCD pCCD
F F F F
# CIS RPA RPAx ppRPA ADC
F F F F F
# CIS CID CISD
F F T
# RPA RPAx ppRPA
F F F
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0 evGW qsGW
T F F
F F F
# G0T0 evGT qsGT
T F F
F F F
# MCMP2
F

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 5 5 0 0
1 1 1 0 0
# Znuc x y z
Ne 0.0 0.0 0.0
He 0.0 0.0 0.0

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1
Ne 0.0000000000 0.0000000000 0.0000000000
He 0.0000000000 0.0000000000 0.0000000000

44
input/weight
View File

@ -1,39 +1,9 @@
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
1 3
S 3
1 38.3600000 0.0238090
2 5.7700000 0.1548910
3 1.2400000 0.4699870
S 1
1 2.8360000 1.0000000
S 1
1 0.3782000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
1 0.2976000 1.0000000
P 1
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000
1 1.2750000 1.0000000

204
src/QuAcK/CISD.f90 Normal file
View File

@ -0,0 +1,204 @@
subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
! Perform configuration interaction with singles and doubles
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
logical :: dump_trans = .false.
integer :: i,j,k,l
integer :: a,b,c,d
integer :: ia,jb,iajb,kcld
integer :: ishift,jshift
integer :: ispin
integer :: nS
integer :: nD
integer :: nSD
double precision,external :: Kronecker_delta
double precision,allocatable :: H(:,:),Omega(:)
! Hello world
write(*,*)
write(*,*)'******************************************************'
write(*,*)'| Configuration Interaction with Singles and Doubles |'
write(*,*)'******************************************************'
write(*,*)
! Compute CIS matrix
if(singlet_manifold) then
ispin = 1
! Dimensions
nS = (nO - nC)*(nV - nR)
nD = (nO - nC)*(nO - nC + 1)/2*(nV - nR)*(nV - nR + 1)/2
nSD = 1 + nS + nD
print*,'nS = ',nS
print*,'nD = ',nD
print*,'nSD = ',nSD
! Memory allocation
allocate(H(nSD,nSD),Omega(nSD))
! 0D block
ishift = 0
jshift = 1 + nS
iajb = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i,nO
do b=a,nV-nR
iajb = iajb + 1
H(ishift+1,jshift+iajb) = ERI(i,j,nO+a,nO+b)
end do
end do
end do
end do
! SS block
ishift = 1
jshift = 1
ia = 0
jb = 0
do i=nC+1,nO
do a=1,nV-nR
ia = ia + 1
do j=nC+1,nO
do b=1,nV-nR
jb = jb + 1
H(ishift+ia,jshift+jb) &
= Kronecker_delta(i,j)*Kronecker_delta(a,b)*(eHF(nO+a) - eHF(i)) &
+ ERI(nO+a,j,i,nO+b) - ERI(nO+a,j,nO+b,i)
end do
end do
end do
end do
! SD block
ishift = 1
jshift = 1 + nS
ia = 0
kcld = 0
do i=nC+1,nO
do a=1,nV-nR
ia = ia + 1
do k=nC+1,nO
do c=1,nV-nR
do l=k,nO
do d=c,nV-nR
kcld = kcld + 1
H(ishift+ia,jshift+kcld) &
= Kronecker_delta(i,k)*(ERI(nO+a,l,nO+c,nO+d) - ERI(nO+a,l,nO+d,nO+c)) &
- Kronecker_delta(i,l)*(ERI(nO+a,k,nO+c,nO+d) - ERI(nO+a,k,nO+d,nO+c)) &
- Kronecker_delta(a,c)*(ERI(k,l,i,nO+d) - ERI(k,l,nO+d,i)) &
+ Kronecker_delta(a,d)*(ERI(k,l,i,nO+c) - ERI(k,l,nO+c,i))
end do
end do
end do
end do
end do
end do
! DD block
ishift = 1 + nS
jshift = 1 + nS
iajb = 0
kcld = 0
do i=nC+1,nO
do a=1,nV-nR
do j=i,nO
do b=a,nV-nR
iajb = iajb + 1
do k=nC+1,nO
do c=1,nV-nR
do l=k,nO
do d=c,nV-nR
kcld = kcld + 1
! H(ishift+iajb,jshift+kcld) &
! = Kronecker_delta(i,k)*(ERI(a,l,c,d) - ERI(a,l,d,c)) &
! - Kronecker_delta(i,l)*(ERI(a,k,c,d) - ERI(a,k,d,c)) &
! - Kronecker_delta(a,c)*(ERI(k,l,i,d) - ERI(k,l,d,i)) &
! + Kronecker_delta(a,d)*(ERI(k,l,i,c) - ERI(k,l,c,i))
end do
end do
end do
end do
end do
end do
end do
end do
call diagonalize_matrix(nSD,H,Omega)
call print_excitation('CISD ',ispin,nS,Omega)
if(dump_trans) then
print*,'Singlet CISD transition vectors'
call matout(nSD,nSD,H)
write(*,*)
endif
endif
! if(triplet_manifold) then
! ispin = 2
!
! call diagonalize_matrix(nSD,H,Omega)
! call print_excitation('CISD ',ispin,nSD,Omega)
! if(dump_trans) then
! print*,'Triplet CIS transition vectors'
! call matout(nSD,nSD,H)
! write(*,*)
! endif
! endif
end subroutine CISD

61
src/QuAcK/QuAcK.f90
View File

@ -8,8 +8,9 @@ program QuAcK
logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical :: doCIS,doRPA,doRPAx
logical :: doppRPA,doADC
logical :: doCIS,doCID,doCISD
logical :: doRPA,doRPAx,doppRPA
logical :: doADC
logical :: doG0F2,doevGF2,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW
logical :: doG0T0,doevGT,doqsGT
@ -61,6 +62,8 @@ program QuAcK
double precision :: start_CCD ,end_CCD ,t_CCD
double precision :: start_CCSD ,end_CCSD ,t_CCSD
double precision :: start_CIS ,end_CIS ,t_CIS
double precision :: start_CID ,end_CID ,t_CID
double precision :: start_CISD ,end_CISD ,t_CISD
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_RPAx ,end_RPAx ,t_RPAx
double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA
@ -128,8 +131,8 @@ program QuAcK
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doCIS,doCID,doCISD, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
@ -491,6 +494,38 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Compute CID excitations
!------------------------------------------------------------------------
if(doCID) then
call cpu_time(start_CID)
! call CID(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CID)
t_CID = end_CID - start_CID
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CID,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute CISD excitations
!------------------------------------------------------------------------
if(doCISD) then
call cpu_time(start_CISD)
call CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_CISD)
t_CISD = end_CISD - start_CISD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CISD,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute (direct) RPA excitations
!------------------------------------------------------------------------
@ -546,18 +581,18 @@ program QuAcK
! Compute ADC excitations
!------------------------------------------------------------------------
if(doADC) then
! if(doADC) then
call cpu_time(start_ADC)
call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO)
call cpu_time(end_ADC)
! call cpu_time(start_ADC)
! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
! nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO)
! call cpu_time(end_ADC)
t_ADC = end_ADC - start_ADC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ADC = ',t_ADC,' seconds'
write(*,*)
! t_ADC = end_ADC - start_ADC
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ADC = ',t_ADC,' seconds'
! write(*,*)
end if
! end if
!------------------------------------------------------------------------
! Compute G0F2 electronic binding energies

2
src/QuAcK/linear_response_A_matrix.f90
View File

@ -16,7 +16,7 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: Kronecker_delta
double precision,external :: Kronecker_delta
integer :: i,j,a,b,ia,jb

3
src/QuAcK/ppRPA.f90
View File

@ -56,9 +56,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
nOO = nO*(nO+1)/2
nVV = nV*(nV+1)/2
! print*,'nOO = ',nOO
! print*,'nVV = ',nVV
! Memory allocation
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &

25
src/QuAcK/read_methods.f90
View File

@ -2,8 +2,8 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doCIS,doCID,doCISD, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
@ -19,7 +19,8 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical,intent(out) :: doCIS,doRPA,doRPAx,doppRPA,doADC
logical,intent(out) :: doCIS,doCID,doCISD
logical,intent(out) :: doRPA,doRPAx,doppRPA
logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW
logical,intent(out) :: doG0T0,doevGT,doqsGT
@ -53,10 +54,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
do_pCCD = .false.
doCIS = .false.
doCID = .false.
doCISD = .false.
doRPA = .false.
doRPAx = .false.
doppRPA = .false.
doADC = .false.
doG0F2 = .false.
doevGF2 = .false.
@ -108,12 +111,16 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
! Read excited state methods
read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5
read(1,*) answer1,answer2,answer3
if(answer1 == 'T') doCIS = .true.
if(answer2 == 'T') doRPA = .true.
if(answer3 == 'T') doRPAx = .true.
if(answer4 == 'T') doppRPA = .true.
if(answer5 == 'T') doADC = .true.
if(answer2 == 'T') doCID = .true.
if(answer3 == 'T') doCISD = .true.
read(1,*)
read(1,*) answer1,answer2,answer3
if(answer1 == 'T') doRPA = .true.
if(answer2 == 'T') doRPAx = .true.
if(answer3 == 'T') doppRPA = .true.
! Read Green function methods