LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

clean up

This commit is contained in:
Pierre-Francois Loos 2023-07-27 11:38:38 +02:00
parent 0a212628c3
commit d738e335ba
49 changed files with 49 additions and 229 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CC/BCCD.f90
View File

@ -233,4 +233,4 @@ subroutine BCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
write(*,*)
end subroutine BCCD
end subroutine

2
src/CC/CCD.f90
View File

@ -272,4 +272,4 @@ subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,e
if(do_DIP_EOM_CC_2h) call DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t)
end subroutine CCD
end subroutine

2
src/CC/CCD_correlation_energy.f90
View File

@ -33,4 +33,4 @@ subroutine CCD_correlation_energy(nC,nO,nV,nR,OOVV,t2,EcCCD)
EcCCD = 0.25d0*EcCCD
end subroutine CCD_correlation_energy
end subroutine

2
src/CC/CCGW.f90
View File

@ -349,4 +349,4 @@ subroutine CCGW(maxSCF,thresh,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,e)
enddo
write(*,*)'-------------------------------------------------------------------------------'
end subroutine CCGW
end subroutine

2
src/CC/CCSD.f90
View File

@ -302,4 +302,4 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,EN
end if
end subroutine CCSD
end subroutine

2
src/CC/CCSDT.f90
View File

@ -42,4 +42,4 @@ subroutine CCSDT(nC,nO,nV,nR,eO,eV,OOVV,VVVO,VOOO,t1,t2,EcCCT)
call form_T(nC,nO,nV,nR,delta_OOOVVV,ub,ubb,EcCCT)
end subroutine CCSDT
end subroutine

2
src/CC/CCSD_Ec_nc.f90
View File

@ -53,4 +53,4 @@ subroutine CCSD_Ec_nc(nO,nV,t1,t2,Fov,OOVV,EcCCSD)
end do
end do
end subroutine CCSD_Ec_nc
end subroutine

2
src/CC/CCSD_correlation_energy.f90
View File

@ -33,4 +33,4 @@ subroutine CCSD_correlation_energy(nC,nO,nV,nR,OOVV,tau,EcCCSD)
EcCCSD = 0.5d0*EcCCSD
end subroutine CCSD_correlation_energy
end subroutine

2
src/CC/DCD.f90
View File

@ -280,4 +280,4 @@ subroutine DCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
endif
end subroutine DCD
end subroutine

2
src/CC/DEA_EOM_CCD_2p.f90
View File

@ -113,4 +113,4 @@ subroutine DEA_EOM_CCD_2p(nC,nO,nV,nR,eV,OOVV,VVVV,t)
call print_excitation('DEA-EOM-CCD ',3,nVV,Om)
end subroutine DEA_EOM_CCD_2p
end subroutine

2
src/CC/DIP_EOM_CCD_2h.f90
View File

@ -113,4 +113,4 @@ subroutine DIP_EOM_CCD_2h(nC,nO,nV,nR,eO,OOVV,OOOO,t)
call print_excitation('DIP-EOM-CCD ',3,nOO,Om)
end subroutine DIP_EOM_CCD_2h
end subroutine

2
src/CC/EE_EOM_CCD_1h1p.f90
View File

@ -147,4 +147,4 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
end if
end subroutine EE_EOM_CCD_1h1p
end subroutine

2
src/CC/MP3_correlation_energy.f90
View File

@ -35,4 +35,4 @@ subroutine MP3_correlation_energy(nC,nO,nV,nR,OOVV,t2,v,delta_OOVV,EcMP3)
EcMP3 = 0.25d0*EcMP3
end subroutine MP3_correlation_energy
end subroutine

2
src/CC/antisymmetrize_ERI.f90
View File

@ -43,4 +43,4 @@ subroutine antisymmetrize_ERI(ispin,nBas,ERI,db_ERI)
endif
end subroutine antisymmetrize_ERI
end subroutine

2
src/CC/crCCD.f90
View File

@ -202,4 +202,4 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*)'----------------------------------------------------'
write(*,*)
end subroutine crCCD
end subroutine

2
src/CC/drCCD.f90
View File

@ -195,4 +195,4 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF
write(*,*)'----------------------------------------------------'
write(*,*)
end subroutine drCCD
end subroutine

2
src/CC/form_EOM_one_body.f90
View File

@ -103,4 +103,4 @@ subroutine form_EOM_one_body(nO,nV,foo,fov,fvv,t1,t2,OOOV,OOVV,OVVV,cFoo,cFov,cF
end do
end do
end subroutine form_EOM_one_body
end subroutine

2
src/CC/form_EOM_two_body.f90
View File

@ -282,4 +282,4 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv,
end do
end do
end subroutine form_EOM_two_body
end subroutine

2
src/CC/form_FB.f90
View File

@ -86,4 +86,4 @@ subroutine form_FB(nC,nO,nV,nR,foo,fvv,fov,OOOV,OOVV,OVVV,t,FooB,FvvB,FovB)
enddo
enddo
end subroutine form_FB
end subroutine

2
src/CC/form_T.f90
View File

@ -43,4 +43,4 @@ subroutine form_T(nC,nO,nV,nR,delta_OOOVVV,ub,ubb,EcCCT)
EcCCT = - EcCCT/36d0
end subroutine form_T
end subroutine

2
src/CC/form_X.f90
View File

@ -89,4 +89,4 @@ subroutine form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4)
enddo
enddo
end subroutine form_X
end subroutine

2
src/CC/form_abh.f90
View File

@ -102,4 +102,4 @@ subroutine form_abh(nC,nO,nV,nR,OOOO,OVOO,OOVV,VVVV,VOVV,OVVO,OVVV,t1,tau,aoooo,
end do
end do
end subroutine form_abh
end subroutine

2
src/CC/form_cF_nc.f90
View File

@ -104,4 +104,4 @@ subroutine form_cF_nc(nO,nV,t1,taus,Foo,Fov,Fvv,OOOV,OOVV,OVVV,cFoo,cFov,cFvv)
end do
end do
end subroutine form_cF_nc
end subroutine

2
src/CC/form_cW_nc.f90
View File

@ -111,4 +111,4 @@ subroutine form_cW_nc(nO,nV,t1,t2,tau,OOOO,OOOV,OOVO,OOVV,OVVO,OVVV,VOVV,VVVV,cW
end do
end subroutine form_cW_nc
end subroutine

2
src/CC/form_crossed_ring_r.f90
View File

@ -67,4 +67,4 @@ subroutine form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2)
end do
end do
end subroutine form_crossed_ring_r
end subroutine

2
src/CC/form_delta_OOOVVV.f90
View File

@ -34,4 +34,4 @@ subroutine form_delta_OOOVVV(nC,nO,nV,nR,eO,eV,delta)
enddo
enddo
end subroutine form_delta_OOOVVV
end subroutine

2
src/CC/form_delta_OOVV.f90
View File

@ -30,4 +30,4 @@ subroutine form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta)
enddo
enddo
end subroutine form_delta_OOVV
end subroutine

2
src/CC/form_delta_OV.f90
View File

@ -27,4 +27,4 @@ subroutine form_delta_OV(nC,nO,nV,nR,eO,eV,delta)
enddo
enddo
end subroutine form_delta_OV
end subroutine

2
src/CC/form_g.f90
View File

@ -50,4 +50,4 @@ subroutine form_g(nC,nO,nV,nR,hvv,hoo,VOVV,OOOV,t1,gvv,goo)
end do
end do
end subroutine form_g
end subroutine

2
src/CC/form_h.f90
View File

@ -76,4 +76,4 @@ subroutine form_h(nC,nO,nV,nR,eO,eV,OOVV,t1,tau,hvv,hoo,hvo)
! print*,'hvv',hvv
end subroutine form_h
end subroutine

2
src/CC/form_ladder_r.f90
View File

@ -70,4 +70,4 @@ subroutine form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2)
end do
end do
end subroutine form_ladder_r
end subroutine

2
src/CC/form_r1.f90
View File

@ -74,4 +74,4 @@ subroutine form_r1(nC,nO,nV,nR,OVVO,OVVV,OOOV,hvv,hoo,hvo,t1,t2,tau,r1)
end do
end do
end subroutine form_r1
end subroutine

2
src/CC/form_r1_nc.f90
View File

@ -79,4 +79,4 @@ subroutine form_r1_nc(nO,nV,t1,t2,delta_ov,Fov,cFoo,cFov,cFvv,OOVO,OVOV,OVVV,r1)
r1(:,:) = delta_ov(:,:)*t1(:,:) - r1(:,:)
end subroutine form_r1_nc
end subroutine

2
src/CC/form_r2.f90
View File

@ -136,4 +136,4 @@ subroutine form_r2(nC,nO,nV,nR,OOVV,OVOO,OVVV,OVVO,gvv,goo,aoooo,bvvvv,hovvo,t1,
end do
end do
end subroutine form_r2
end subroutine

2
src/CC/form_r2_nc.f90
View File

@ -106,4 +106,4 @@ subroutine form_r2_nc(nO,nV,t1,t2,tau,delta_oovv,cFoo,cFov,cFvv,cWoooo,cWvvvv,cW
r2(:,:,:,:) = delta_oovv(:,:,:,:)*t2(:,:,:,:) - r2(:,:,:,:)
end subroutine form_r2_nc
end subroutine

2
src/CC/form_ring_r.f90
View File

@ -67,4 +67,4 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2)
end do
end do
end subroutine form_ring_r
end subroutine

2
src/CC/form_tau.f90
View File

@ -31,4 +31,4 @@ subroutine form_tau(nC,nO,nV,nR,t1,t2,tau)
enddo
enddo
end subroutine form_tau
end subroutine

2
src/CC/form_tau_nc.f90
View File

@ -31,4 +31,4 @@ subroutine form_tau_nc(nO,nV,t1,t2,tau)
enddo
enddo
end subroutine form_tau_nc
end subroutine

2
src/CC/form_taus_nc.f90
View File

@ -31,4 +31,4 @@ subroutine form_taus_nc(nO,nV,t1,t2,taus)
enddo
enddo
end subroutine form_taus_nc
end subroutine

2
src/CC/form_u.f90
View File

@ -68,4 +68,4 @@ subroutine form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u)
enddo
enddo
end subroutine form_u
end subroutine

2
src/CC/form_ub.f90
View File

@ -45,4 +45,4 @@ subroutine form_ub(nC,nO,nV,nR,OOVV,t1,ub)
end do
end do
end subroutine form_ub
end subroutine

2
src/CC/form_ubb.f90
View File

@ -64,4 +64,4 @@ subroutine form_ubb(nC,nO,nV,nR,VVVO,VOOO,t2,ubb)
end do
end do
end subroutine form_ubb
end subroutine

2
src/CC/form_v.f90
View File

@ -76,4 +76,4 @@ subroutine form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v)
enddo
enddo
end subroutine form_v
end subroutine

2
src/CC/lCCD.f90
View File

@ -215,4 +215,4 @@ subroutine lCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
write(*,*)'----------------------------------------------------'
write(*,*)
end subroutine lCCD
end subroutine

2
src/CC/pCCD.f90
View File

@ -212,4 +212,4 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
endif
end subroutine pCCD
end subroutine

2
src/CC/rCCD.f90
View File

@ -223,4 +223,4 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,
if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t)
end subroutine rCCD
end subroutine

4
src/LR/print_unrestricted_transition_vectors.f90
View File

@ -49,7 +49,7 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
! Compute <S**2>
call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
! Print details about spin-conserved excitations
@ -167,4 +167,4 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
write(*,*)
end subroutine print_unrestricted_transition_vectors
end subroutine

180
src/LR/unrestricted_S2_expval.f90
View File

@ -1,180 +0,0 @@
subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
! Compute <S**2> for linear response excited states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: ispin
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
integer,intent(in) :: nSa
integer,intent(in) :: nSb
integer,intent(in) :: nSt
integer,intent(in) :: maxS
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt)
! Local variables
integer :: m
integer :: ia,i,a
double precision :: S2_exact
double precision :: S2_gs
double precision,allocatable :: Xa(:,:), Xb(:,:), Ya(:,:), Yb(:,:)
double precision,allocatable :: Xat(:,:),Xbt(:,:),Yat(:,:),Ybt(:,:)
double precision,allocatable :: OO(:,:), OV(:,:), VO(:,:), VV(:,:)
double precision,allocatable :: OOt(:,:),OVt(:,:),VOt(:,:),VVt(:,:)
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: S2(maxS)
! Memory allocation
allocate(OO(nO(1)-nC(1),nO(2)-nC(2)), OV(nO(1)-nC(1),nV(2)-nR(2)), VO(nV(1)-nR(1),nO(2)-nC(2)), VV(nV(1)-nR(1),nV(2)-nR(2)), &
OOt(nO(2)-nC(2),nO(1)-nC(1)),OVt(nV(2)-nR(2),nO(1)-nC(1)),VOt(nO(2)-nC(2),nV(1)-nR(1)),VVt(nV(2)-nR(2),nV(1)-nR(1)))
! Overlap matrix between spin-up and spin-down orbitals
OO(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nC(2)+1:nO(2) ,2)))
OV(:,:) = matmul(transpose(c(:,nC(1)+1:nO(1) ,1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2)))
VO(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nC(2)+1:nO(2) ,2)))
VV(:,:) = matmul(transpose(c(:,nO(1)+1:nBas-nR(1),1)),matmul(S,c(:,nO(2)+1:nBas-nR(2),2)))
OOt(:,:) = transpose(OO(:,:))
OVt(:,:) = transpose(OV(:,:))
VOt(:,:) = transpose(VO(:,:))
VVt(:,:) = transpose(VV(:,:))
!-------------------------!
! <S**2> for ground state !
!-------------------------!
S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
S2_gs = S2_exact + dble(nO(2)) - sum(OO(:,:)**2)
!------------------------------------------!
! <S**2> for spin-conserved-excited states !
!------------------------------------------!
if(ispin == 1) then
allocate(Xa(nO(1)-nC(1),nV(1)-nR(1)), Ya(nO(1)-nC(1),nV(1)-nR(1)), Xb(nO(2)-nC(2),nV(2)-nR(2)), Yb(nO(2)-nC(2),nV(2)-nR(2)), &
Xat(nV(1)-nR(1),nO(1)-nC(1)),Yat(nV(1)-nR(1),nO(1)-nC(1)),Xbt(nV(2)-nR(2),nO(2)-nC(2)),Ybt(nV(2)-nR(2),nO(2)-nC(2)))
do m=1,maxS
ia = 0
do i=nC(1)+1,nO(1)
do a=1,nV(1)-nR(1)
ia = ia + 1
Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia))
Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia))
end do
end do
ia = 0
do i=nC(2)+1,nO(2)
do a=1,nV(2)-nR(2)
ia = ia + 1
Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia))
Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia))
end do
end do
Xat(:,:) = transpose(Xa(:,:))
Xbt(:,:) = transpose(Xb(:,:))
Yat(:,:) = transpose(Ya(:,:))
Ybt(:,:) = transpose(Yb(:,:))
S2(m) = S2_gs &
+ trace_matrix(nV(1),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) &
+ trace_matrix(nV(2),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
- trace_matrix(nO(1),matmul(Xa,matmul(VO,matmul(VOt,Xat)))) &
- trace_matrix(nO(2),matmul(Xb,matmul(OVt,matmul(OV,Xbt)))) &
- 2d0*trace_matrix(nO(1),matmul(OO,matmul(Xb,matmul(VVt,Xat)))) &
- 2d0*trace_matrix(nV(2),matmul(OVt,matmul(Xa,matmul(VO,Yb)))) &
- 2d0*trace_matrix(nV(1),matmul(VO,matmul(Xb,matmul(OVt,Ya)))) &
- trace_matrix(nV(1),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) &
- trace_matrix(nV(2),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) &
+ trace_matrix(nO(1),matmul(Ya,matmul(VO,matmul(VOt,Yat)))) &
+ trace_matrix(nO(2),matmul(Yb,matmul(OVt,matmul(OV,Ybt)))) &
+ 2d0*trace_matrix(nO(1),matmul(Ya,matmul(VV,matmul(Ybt,OOt))))
end do
end if
!------------------------------------------!
! <S**2> for spin-conserved-excited states !
!------------------------------------------!
if(ispin == 2) then
allocate(Xa(nO(1)-nC(1),nV(2)-nR(2)), Ya(nO(1)-nC(1),nV(2)-nR(2)), Xb(nO(2)-nC(2),nV(1)-nR(1)), Yb(nO(2)-nC(2),nV(1)-nR(1)), &
Xat(nV(2)-nR(2),nO(1)-nC(1)),Yat(nV(2)-nR(2),nO(1)-nC(1)),Xbt(nV(1)-nR(1),nO(2)-nC(2)),Ybt(nV(1)-nR(1),nO(2)-nC(2)))
do m=1,maxS
ia = 0
do i=nC(1)+1,nO(1)
do a=1,nV(2)-nR(2)
ia = ia + 1
Xa(i,a) = 0.5d0*(XpY(m,ia) + XmY(m,ia))
Ya(i,a) = 0.5d0*(XpY(m,ia) - XmY(m,ia))
end do
end do
ia = 0
do i=nC(2)+1,nO(2)
do a=1,nV(1)-nR(1)
ia = ia + 1
Xb(i,a) = 0.5d0*(XpY(m,nSa+ia) + XmY(m,nSa+ia))
Yb(i,a) = 0.5d0*(XpY(m,nSa+ia) - XmY(m,nSa+ia))
end do
end do
Xat(:,:) = transpose(Xa(:,:))
Xbt(:,:) = transpose(Xb(:,:))
Yat(:,:) = transpose(Ya(:,:))
Ybt(:,:) = transpose(Yb(:,:))
S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0
S2(m) = S2(m) &
+ trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
- trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) &
+ trace_matrix(nO(2),matmul(Xb,VO))**2 &
+ trace_matrix(nV(2),matmul(Yat,matmul(OO,matmul(OOt,Ya)))) &
+ trace_matrix(nO(1),matmul(Ya,matmul(OVt,matmul(OV,Yat)))) &
+ trace_matrix(nO(1),matmul(Ya,OVt))**2 &
- 2d0*trace_matrix(nO(2),matmul(Xb,VO))*trace_matrix(nO(1),matmul(Ya,OVt)) &
+ trace_matrix(nV(2),matmul(Xat,matmul(OO,matmul(OOt,Xa)))) &
- trace_matrix(nO(1),matmul(Xa,matmul(OVt,matmul(OV,Xat)))) &
+ trace_matrix(nO(1),matmul(Xa,OVt))**2 &
+ trace_matrix(nV(1),matmul(Ybt,matmul(OOt,matmul(OO,Yb)))) &
- trace_matrix(nO(2),matmul(Yb,matmul(VO,matmul(VOt,Ybt)))) &
+ trace_matrix(nV(1),matmul(Ybt,VOt))**2 &
- 2d0*trace_matrix(nO(1),matmul(Xa,OVt))*trace_matrix(nO(2),matmul(Yb,VO))
end do
end if
end subroutine unrestricted_S2_expval

2
src/LR/unrestricted_oscillator_strength.f90
View File

@ -82,4 +82,4 @@ subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS
write(*,*) '---------------------------------------------------------------'
write(*,*)
end subroutine unrestricted_oscillator_strength
end subroutine