LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-25 06:32:19 +02:00
Code Releases Activity

slow CIS(D) correction (in spinorbital basis)

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-10-01 23:05:23 +02:00
parent ac66ae8ea2
commit 69a0b5193d
12 changed files with 287 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.Be
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 3 1 0 0
1 2 2 0 0
# Znuc x y z
Be 0.0 0.0 0.0

2
examples/molecule.He
View File

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

6
input/methods
View File

@ -1,13 +1,13 @@
# RHF UHF MOM
F T F
T F F
# MP2* MP3 MP2-F12
F F F
# CCD CCSD CCSD(T)
F F F
# drCCD rCCD lCCD pCCD
F F F F
# CIS* CID CISD
T F F
# CIS* CIS(D) CID CISD
T T F F
# RPA* RPAx* ppRPA
F F F
# G0F2 evGF2 G0F3 evGF3

23
src/CI/CIS.f90
View File

@ -1,5 +1,4 @@
subroutine CIS(singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF)
! Perform configuration interaction single calculation`
@ -8,8 +7,9 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
! Input variables
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: doCIS_D
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
@ -20,6 +20,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
logical :: dump_matrix = .false.
logical :: dump_trans = .false.
integer :: ispin
integer :: maxS = 10
double precision :: lambda
double precision,allocatable :: A(:,:),Omega(:)
@ -41,7 +42,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
! Compute CIS matrix
if(singlet_manifold) then
if(singlet) then
ispin = 1
call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
@ -61,9 +62,14 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
write(*,*)
endif
! Compute CIS(D) correction
maxS = min(maxS,nS)
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS))
endif
if(triplet_manifold) then
if(triplet) then
ispin = 2
call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A)
@ -83,6 +89,11 @@ subroutine CIS(singlet_manifold,triplet_manifold, &
write(*,*)
endif
! Compute CIS(D) correction
maxS = min(maxS,nS)
if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS))
endif
end subroutine CIS

224
src/CI/D_correction.f90 Normal file
View File

@ -0,0 +1,224 @@
subroutine D_correction(ispin,nBasin,nCin,nOin,nVin,nRin,nSin,maxS,eHF,ERI,w,X)
! Compute the D correction of CIS(D)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBasin
integer,intent(in) :: nCin
integer,intent(in) :: nOin
integer,intent(in) :: nVin
integer,intent(in) :: nRin
integer,intent(in) :: nSin
integer,intent(in) :: maxS
double precision,intent(in) :: eHF(nBasin)
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
double precision,intent(in) :: w(maxS)
double precision,intent(in) :: X(nSin,maxS)
! Local variables
integer :: i,j,k
integer :: a,b,c
integer :: m,ia
integer :: nBas
integer :: nC
integer :: nO
integer :: nV
integer :: nR
double precision,allocatable :: seHF(:)
double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: dbERI(:,:,:,:)
double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:)
double precision,allocatable :: delta(:,:,:,:)
double precision,allocatable :: OOOV(:,:,:,:)
double precision,allocatable :: OOVV(:,:,:,:)
double precision,allocatable :: OVVV(:,:,:,:)
double precision,allocatable :: u(:,:,:,:)
double precision,allocatable :: v(:,:)
double precision,allocatable :: t(:,:,:,:)
double precision,allocatable :: rr(:,:),r(:,:)
double precision :: wD
double precision,external :: Kronecker_delta
! Spatial to spin orbitals
nBas = 2*nBasin
nC = 2*nCin
nO = 2*nOin
nV = 2*nVin
nR = 2*nRin
allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas))
call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF)
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
! Antysymmetrize ERIs
allocate(dbERI(nBas,nBas,nBas,nBas))
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
deallocate(sERI)
! Form energy denominator
allocate(eO(nO),eV(nV))
allocate(delta(nO,nO,nV,nV))
eO(:) = seHF(1:nO)
eV(:) = seHF(nO+1:nBas)
call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta)
deallocate(seHF,eO,eV)
! Create integral batches
allocate(OOOV(nO,nO,nO,nV),OOVV(nO,nO,nV,nV),OVVV(nO,nV,nV,nV))
OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas)
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas,nO+1:nBas)
deallocate(dbERI)
! Memory allocation
allocate(t(nO,nO,nV,nV),r(nO,nV),u(nO,nO,nV,nV),v(nO,nV))
! MP2 guess amplitudes
t(:,:,:,:) = -OOVV(:,:,:,:)/delta(:,:,:,:)
!------------------------------------------------------------------------
! Loop over single excitations
!------------------------------------------------------------------------
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*) ' CIS(D) correction '
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,'(2X,A5,1X,A20,1X,A20,1X,A20)') '#','CIS (eV)','CIS(D) (eV)','Correction (eV)'
write(*,*) '---------------------------------------------------------------------------------------------------'
do m=1,maxS
! Unfold r
allocate(rr(nOin,nVin))
ia = 0
do i=nCin+1,nOin
do a=1,nVin-nRin
ia = ia + 1
rr(i,a) = x(ia,m)
end do
end do
if(ispin == 1) then
do i=nC+1,nO
do a=1,nV-nR
r(i,a) = rr((i+1)/2,(a+1)/2)*Kronecker_delta(mod(i,2),mod(a,2))
end do
end do
elseif(ispin == 2) then
do i=nC+1,nO
do a=1,nV-nR
r(i,a) = rr((i+1)/2,(a+1)/2)*Kronecker_delta(mod(i,2),mod(a+1,2))
end do
end do
else
print*,'!!! CIS(D) must be for singlet or triplet !!!'
stop
end if
deallocate(rr)
! Compute u array
u(:,:,:,:) = 0d0
do i=nC+1,nO
do a=1,nV-nR
do j=nC+1,nO
do b=1,nV-nR
do c=1,nV-nR
u(i,j,a,b) = u(i,j,a,b) + OVVV(i,c,a,b)*r(j,c) - OVVV(j,c,a,b)*r(i,c)
end do
do k=nC+1,nO
u(i,j,a,b) = u(i,j,a,b) + OOOV(i,j,k,a)*r(k,b) - OOOV(i,j,k,b)*r(k,a)
end do
end do
end do
end do
end do
! Compute v array
v(:,:) = 0d0
do i=nC+1,nO
do j=nC+1,nO
do k=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
do c=1,nV-nR
v(i,a) = v(i,a) + OOVV(j,k,b,c)*(r(i,b)*t(j,k,c,a) + r(j,a)*t(i,k,c,b) + 2d0*r(j,b)*t(i,k,a,c))
end do
end do
end do
end do
end do
end do
v(:,:) = 0.5d0*v(:,:)
! Compute CIS(D) correction to CIS excitation energies
wD = 0d0
do i=nC+1,nO
do a=1,nV-nR
do j=nC+1,nO
do b=1,nV-nR
wD = wD - 0.25d0*u(i,j,a,b)**2/(delta(i,j,a,b) - w(m))
enddo
enddo
wD = wD + r(i,a)*v(i,a)
enddo
enddo
wD = 0.5d0*wD
! Flush results
write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6)') m,w(m)*HaToeV,(w(m)+wD)*HaToeV,wD*HaToeV
end do
write(*,*) '---------------------------------------------------------------------------------------------------'
write(*,*)
!------------------------------------------------------------------------
! End of loop over single excitations
!------------------------------------------------------------------------
end subroutine D_correction

10
src/CI/UCIS.f90
View File

@ -37,7 +37,9 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
integer :: nS_ab,nS_ba,nS_sf
double precision,allocatable :: A_sf(:,:)
double precision,allocatable :: Z_sf(:,:)
double precision,allocatable :: Omega_sf(:)
integer ,allocatable :: order(:)
! Hello world
@ -103,7 +105,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
nS_sf = nS_ab + nS_ba
allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf))
allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf),Z_sf(nS_sf,nS_sf))
call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, &
ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf)
@ -114,7 +116,11 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
write(*,*)
endif
call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
! call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
allocate(order(nS_sf))
call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf)
call quick_sort(Omega_sf,order(:),nS_sf)
call print_excitation('UCIS ',6,nS_sf,Omega_sf)
if(dump_trans) then

7
src/HF/print_UHF.f90
View File

@ -26,7 +26,11 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = e(nO(ispin),ispin)
LUMO(ispin) = e(nO(ispin)+1,ispin)
if(nO(ispin) < nBas) then
LUMO(ispin) = e(nO(ispin)+1,ispin)
else
LUMO(ispin) = 0d0
end if
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
@ -37,7 +41,6 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF)
! Dump results
write(*,*)
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A40)') ' Summary '

6
src/QuAcK/QuAcK.f90
View File

@ -9,7 +9,7 @@ program QuAcK
logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical :: doCIS,doCID,doCISD
logical :: doCIS,doCIS_D,doCID,doCISD
logical :: doRPA,doRPAx,doppRPA
logical :: doADC
logical :: doG0F2,doevGF2,doG0F3,doevGF3
@ -158,7 +158,7 @@ program QuAcK
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCID,doCISD, &
doCIS,doCIS_D,doCID,doCISD, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
@ -621,7 +621,7 @@ program QuAcK
else
call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF)
call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF)
end if
call cpu_time(end_CIS)

13
src/QuAcK/read_methods.f90
View File

@ -2,7 +2,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCID,doCISD, &
doCIS,doCIS_D,doCID,doCISD, &
doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
@ -19,7 +19,7 @@ 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,doCID,doCISD
logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD
logical,intent(out) :: doRPA,doRPAx,doppRPA
logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW
@ -54,6 +54,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
do_pCCD = .false.
doCIS = .false.
doCIS_D = .false.
doCID = .false.
doCISD = .false.
@ -111,10 +112,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
! Read excited state methods
read(1,*)
read(1,*) answer1,answer2,answer3
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doCIS = .true.
if(answer2 == 'T') doCID = .true.
if(answer3 == 'T') doCISD = .true.
if(answer2 == 'T') doCIS_D = .true.
if(answer3 == 'T') doCID = .true.
if(answer4 == 'T') doCISD = .true.
if(doCIS_D) doCIS = .true.
read(1,*)
read(1,*) answer1,answer2,answer3

18
src/RPA/unrestricted_linear_response_A_matrix.f90
View File

@ -121,6 +121,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
end do
end do
print*,'nSa,nSb,nSt',nSa,nSb,nSt
call matout(nSt,nSt,A_lr)
end if
!-----------------------------------------------
@ -141,9 +144,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do j=nC(1)+1,nO(1)
do b=nO(2)+1,nBas-nR(2)
jb = jb + 1
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)'
A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !&
! - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i)
end do
end do
@ -161,15 +164,18 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
do b=nO(1)+1,nBas-nR(1)
jb = jb + 1
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI_abab(j,a,i,b)
! print*,'ia,jb,A(ia,jb) = ',i,a,j,b,ia,jb,A_lr(ia,jb)
A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !&
! - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb)
end do
end do
end do
end do
print*,'nSa,nSb,nSt',nSa,nSb,nSt
call matout(nSt,nSt,A_lr)
end if

4
src/eDFT/eDFT.f90
View File

@ -107,8 +107,8 @@ program eDFT
! Allocate ensemble weights
allocate(wEns(maxEns),occnum(nBas,nspin,maxEns))
call read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice)
call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice)
!------------------------------------------------------------------------
! Read one- and two-electron integrals

6
src/eDFT/read_options.f90 → src/eDFT/read_options_dft.f90
View File

@ -1,5 +1,5 @@
subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice)
subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, &
maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice)
! Read DFT options
@ -201,4 +201,4 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_
close(unit=1)
end subroutine read_options
end subroutine read_options_dft