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

clean up + AC

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-01-14 14:44:01 +01:00
parent c6ec8f91e9
commit 51016fd5c1
12 changed files with 145 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
input/methods
View File

@ -2,14 +2,14 @@
T F F T F F
# MP2 MP3 MP2-F12 # MP2 MP3 MP2-F12
F F F F F F
# CCD CCSD CCSD(T) # CCD CCSD CCSD(T) ringCCD ladderCCD
F F T F F F T T
# CIS RPA TDHF ppRPA ADC # CIS RPA TDHF ppRPA ADC
F T T F F F T T T F
# GF2 GF3 # GF2 GF3
F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
F F F T F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

16
src/QuAcK/CCSD.f90
View File

@ -223,14 +223,14 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF)
end if end if
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*)' CCSD energy ' write(*,*)' CCSD energy '
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD
write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD write(*,'(1X,A20,1X,F10.6)')' Ec(CCSD) = ',EcCCSD
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)
! Deallocate memory ! Deallocate memory

30
src/QuAcK/Ec_AC.f90
View File

@ -1,4 +1,4 @@
subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC)
! Compute the correlation energy via the adiabatic connection formula ! Compute the correlation energy via the adiabatic connection formula
@ -12,6 +12,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
integer,intent(in) :: nBas,nC,nO,nV,nR,nS integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XpY(nS,nS)
double precision,intent(in) :: XmY(nS,nS)
! Local variables ! Local variables
@ -20,7 +21,10 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
double precision :: delta_spin double precision :: delta_spin
double precision :: delta_dRPA double precision :: delta_dRPA
double precision,allocatable :: P(:,:) double precision,allocatable :: P(:,:)
double precision,allocatable :: V(:,:) double precision,allocatable :: Ap(:,:)
double precision,allocatable :: Bp(:,:)
double precision,allocatable :: X(:,:)
double precision,allocatable :: Y(:,:)
double precision,external :: trace_matrix double precision,external :: trace_matrix
! Output variables ! Output variables
@ -40,7 +44,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
! Memory allocation ! Memory allocation
allocate(P(nS,nS),V(nS,nS)) allocate(P(nS,nS),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS))
! Compute P = (X+Y)(X+Y) - 1 ! Compute P = (X+Y)(X+Y) - 1
@ -50,7 +54,7 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
P(ia,ia) = P(ia,ia) - 1d0 P(ia,ia) = P(ia,ia) - 1d0
enddo enddo
! Compute Viajb = (ia|bj) ! Compute Aiajb = (ia|bj) and Biajb = (ia|jb)
ia = 0 ia = 0
do i=nC+1,nO do i=nC+1,nO
@ -61,16 +65,28 @@ subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC)
do b=nO+1,nBas-nR do b=nO+1,nBas-nR
jb = jb + 1 jb = jb + 1
V(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j)
Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,b,a)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
! Compute Tr(VP) ! Compute Tr(A x P)
EcAC = trace_matrix(nS,matmul(V,P)) ! EcAC = trace_matrix(nS,matmul(Ap,P))
! print*,'EcAC =',EcAC
X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
+ trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
- trace_matrix(nS,Ap)
! print*,'EcAC =',EcAC
end subroutine Ec_AC end subroutine Ec_AC

25
src/QuAcK/G0W0.f90
View File

@ -38,6 +38,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
double precision,allocatable :: Omega(:,:) double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:)
@ -79,7 +80,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Memory allocation ! Memory allocation
allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin), & allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin)) rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
AC = .true. AC = .true.
@ -88,7 +89,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
! Compute linear response ! Compute linear response
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
@ -129,11 +130,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if end if
@ -146,11 +147,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
EcBSE(ispin) = 0d0 EcBSE(ispin) = 0d0
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
end if end if
@ -192,13 +193,13 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
lambda = rAC(iAC) lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)
@ -230,13 +231,13 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, &
lambda = rAC(iAC) lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, &
rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin)

53
src/QuAcK/QuAcK.f90
View File

@ -7,6 +7,7 @@ program QuAcK
logical :: doRHF,doUHF,doMOM logical :: doRHF,doUHF,doMOM
logical :: doMP2,doMP3,doMP2F12 logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doCCSD,doCCSDT logical :: doCCD,doCCSD,doCCSDT
logical :: do_ring_CCD,do_ladder_CCD
logical :: doCIS,doRPA,doTDHF logical :: doCIS,doRPA,doTDHF
logical :: doppRPA,doADC logical :: doppRPA,doADC
logical :: doGF2,doGF3 logical :: doGF2,doGF3
@ -108,14 +109,15 @@ program QuAcK
! Which calculations do you want to do? ! Which calculations do you want to do?
call read_methods(doRHF,doUHF,doMOM, & call read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, & doCCD,doCCSD,doCCSDT, &
doCIS,doRPA,doTDHF, & do_ring_CCD,do_ladder_CCD, &
doppRPA,doADC, & doCIS,doRPA,doTDHF, &
doGF2,doGF3, & doppRPA,doADC, &
doG0W0,doevGW,doqsGW, & doGF2,doGF3, &
doG0T0,doevGT,doqsGT, & doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2) doMCMP2)
! Read options for methods ! Read options for methods
@ -351,10 +353,7 @@ program QuAcK
if(doCCD) then if(doCCD) then
call cpu_time(start_CCD) call cpu_time(start_CCD)
! call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
! call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
@ -380,6 +379,38 @@ program QuAcK
end if end if
!------------------------------------------------------------------------
! Perform ring CCD calculation
!------------------------------------------------------------------------
if(do_ring_CCD) then
call cpu_time(start_CCD)
call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ring CCD = ',t_CCD,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform ladder CCD calculation
!------------------------------------------------------------------------
if(do_ladder_CCD) then
call cpu_time(start_CCD)
call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Compute CIS excitations ! Compute CIS excitations
!------------------------------------------------------------------------ !------------------------------------------------------------------------

15
src/QuAcK/RPA.f90
View File

@ -29,6 +29,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
integer :: ispin integer :: ispin
double precision,allocatable :: Omega(:,:) double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision :: rho double precision :: rho
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
@ -65,7 +66,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
! Memory allocation ! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
AC = .true. AC = .true.
allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin)) allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
@ -77,7 +78,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
ispin = 1 ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
endif endif
@ -89,7 +90,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
ispin = 2 ispin = 2
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
endif endif
@ -131,9 +132,9 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
lambda = rAC(iAC) lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
@ -165,9 +166,9 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E
lambda = rAC(iAC) lambda = rAC(iAC)
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, &
EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)

23
src/QuAcK/TDHF.f90
View File

@ -29,6 +29,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
integer :: ispin integer :: ispin
double precision,allocatable :: Omega(:,:) double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision :: rho double precision :: rho
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
@ -65,7 +66,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
! Memory allocation ! Memory allocation
allocate(Omega(nS,nspin),XpY(nS,nS,nspin)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
AC = .true. AC = .true.
allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin)) allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin))
@ -77,7 +78,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
ispin = 1 ispin = 1
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif endif
@ -89,7 +90,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
ispin = 2 ispin = 2
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin))
endif endif
@ -129,15 +130,11 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
do iAC=1,nAC do iAC=1,nAC
lambda = rAC(iAC) lambda = rAC(iAC)
! call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)
@ -167,15 +164,11 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,
do iAC=1,nAC do iAC=1,nAC
lambda = rAC(iAC) lambda = rAC(iAC)
! call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, &
! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin)

11
src/QuAcK/ladder_CCD.f90
View File

@ -182,12 +182,13 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
endif endif
! Moller-Plesset energies
write(*,*) write(*,*)
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2 write(*,*)'----------------------------------------------------'
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3 write(*,*)' ladder-CCD energy '
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 write(*,*)'----------------------------------------------------'
write(*,'(1X,A30,1X,F15.10)')' E(ladder-CCD) = ',ECCD
write(*,'(1X,A30,1X,F10.6)')' Ec(ladder-CCSD) = ',EcCCD
write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)
end subroutine ladder_CCD end subroutine ladder_CCD

9
src/QuAcK/linear_response.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY) subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY)
! Compute linear response ! Compute linear response
@ -22,7 +22,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
! Output variables ! Output variables
double precision,intent(out) :: EcRPA double precision,intent(out) :: EcRPA
double precision,intent(out) :: Omega(nS),XpY(nS,nS) double precision,intent(out) :: Omega(nS)
double precision,intent(out) :: XpY(nS,nS)
double precision,intent(out) :: XmY(nS,nS)
! Memory allocation ! Memory allocation
@ -83,6 +85,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
XpY = matmul(transpose(Z),AmBSq) XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(abs(Omega)),XpY) call DA(nS,1d0/sqrt(abs(Omega)),XpY)
XmY = matmul(transpose(Z),AmBSq)
call DA(nS,sqrt(abs(Omega)),XmY)
! print*,'X+Y' ! print*,'X+Y'
! call matout(nS,nS,XpY) ! call matout(nS,nS,XpY)

2
src/QuAcK/print_excitation.f90
View File

@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
! Input variables ! Input variables
character*4,intent(in) :: method character*6,intent(in) :: method
integer,intent(in) :: ispin,nS integer,intent(in) :: ispin,nS
double precision,intent(in) :: Omega(nS) double precision,intent(in) :: Omega(nS)

31
src/QuAcK/read_methods.f90
View File

@ -1,11 +1,12 @@
subroutine read_methods(doRHF,doUHF,doMOM, & subroutine read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, & doCCD,doCCSD,doCCSDT, &
doCIS,doRPA,doTDHF, & do_ring_CCD,do_ladder_CCD, &
doppRPA,doADC, & doCIS,doRPA,doTDHF, &
doGF2,doGF3, & doppRPA,doADC, &
doG0W0,doevGW,doqsGW, & doGF2,doGF3, &
doG0T0,doevGT,doqsGT, & doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2) doMCMP2)
! Read desired methods ! Read desired methods
@ -17,6 +18,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
logical,intent(out) :: doRHF,doUHF,doMOM logical,intent(out) :: doRHF,doUHF,doMOM
logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCCD,doCCSD,doCCSDT logical,intent(out) :: doCCD,doCCSD,doCCSDT
logical,intent(out) :: do_ring_CCD,do_ladder_CCD
logical,intent(out) :: doCIS,doRPA,doTDHF,doppRPA,doADC logical,intent(out) :: doCIS,doRPA,doTDHF,doppRPA,doADC
logical,intent(out) :: doGF2,doGF3 logical,intent(out) :: doGF2,doGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0W0,doevGW,doqsGW
@ -45,6 +47,9 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
doCCSD = .false. doCCSD = .false.
doCCSDT = .false. doCCSDT = .false.
do_ring_CCD = .false.
do_ladder_CCD = .false.
doCIS = .false. doCIS = .false.
doRPA = .false. doRPA = .false.
doTDHF = .false. doTDHF = .false.
@ -83,10 +88,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
! Read CC methods ! Read CC methods
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3 read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') doCCD = .true. if(answer1 == 'T') doCCD = .true.
if(answer2 == 'T') doCCSD = .true. if(answer2 == 'T') doCCSD = .true.
if(answer3 == 'T') doCCSDT = .true. if(answer3 == 'T') doCCSDT = .true.
if(answer4 == 'T') do_ring_CCD = .true.
if(answer5 == 'T') do_ladder_CCD = .true.
! Read excited state methods ! Read excited state methods

11
src/QuAcK/ring_CCD.f90
View File

@ -182,12 +182,13 @@ subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
endif endif
! Moller-Plesset energies
write(*,*) write(*,*)
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2 write(*,*)'----------------------------------------------------'
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3 write(*,*)' ring-CCD energy '
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 write(*,*)'----------------------------------------------------'
write(*,'(1X,A30,1X,F15.10)')' E(ring-CCD) = ',ECCD
write(*,'(1X,A30,1X,F15.10)')' Ec(ring-CCD) = ',EcCCD
write(*,*)'----------------------------------------------------'
write(*,*) write(*,*)
end subroutine ring_CCD end subroutine ring_CCD