PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

AO <-> MO modifs

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-08 22:58:55 +01:00
parent 67f8fcd0da
commit 94cc75309e
23 changed files with 183 additions and 205 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/AOtoMO/AOtoMO_transform.f90
View File

@ -1,4 +1,4 @@
subroutine AOtoMO_transform(nBas,c,A)
subroutine AOtoMO_transform(nBas,c,A,B)
! Perform AO to MO transformation of a matrix A for given coefficients c
@ -8,11 +8,12 @@ subroutine AOtoMO_transform(nBas,c,A)
integer,intent(in) :: nBas
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: A(nBas,nBas)
! Output variables
double precision,intent(inout):: A(nBas,nBas)
double precision,intent(out) :: B(nBas,nBas)
A = matmul(transpose(c),matmul(A,c))
B = matmul(transpose(c),matmul(A,c))
end subroutine

4
src/AOtoMO/AOtoMO_transform_GHF.f90
View File

@ -10,11 +10,11 @@ subroutine AOtoMO_transform_GHF(nBas,nBas2,Ca,Cb,A,B)
integer,intent(in) :: nBas2
double precision,intent(in) :: Ca(nBas,nBas2)
double precision,intent(in) :: Cb(nBas,nBas2)
double precision,intent(inout):: A(nBas,nBas)
double precision,intent(in) :: A(nBas,nBas)
! Output variables
double precision,intent(inout):: B(nBas2,nBas2)
double precision,intent(out) :: B(nBas2,nBas2)
B = matmul(transpose(Ca),matmul(A,Ca)) &
+ matmul(transpose(Cb),matmul(A,Cb))

10
src/AOtoMO/MOtoAO_transform.f90
View File

@ -1,4 +1,4 @@
subroutine MOtoAO_transform(nBas,S,c,A)
subroutine MOtoAO_transform(nBas,S,c,B,A)
! Perform MO to AO transformation of a matrix A for a given metric S
! and coefficients c
@ -8,7 +8,9 @@ subroutine MOtoAO_transform(nBas,S,c,A)
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas),c(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: B(nBas,nBas)
! Local variables
@ -16,12 +18,12 @@ subroutine MOtoAO_transform(nBas,S,c,A)
! Output variables
double precision,intent(inout):: A(nBas,nBas)
double precision,intent(out) :: A(nBas,nBas)
! Memory allocation
allocate(Sc(nBas,nBas))
Sc = matmul(S,c)
A = matmul(Sc,matmul(A,transpose(Sc)))
A = matmul(Sc,matmul(B,transpose(Sc)))
end subroutine

10
src/GF/qsGF2.f90
View File

@ -73,7 +73,6 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
double precision,allocatable :: K(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: SigCm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
@ -104,7 +103,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
! Memory allocation
allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -157,10 +156,9 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
! Make correlation self-energy Hermitian and transform it back to AO basis
SigCp = 0.5d0*(SigC + transpose(SigC))
SigCm = 0.5d0*(SigC - transpose(SigC))
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO_transform(nBas,S,c,SigCp)
call MOtoAO_transform(nBas,S,c,SigC,SigCp)
! Solve the quasi-particle equation
@ -253,7 +251,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,tr
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis)
! Perform BSE calculation

10
src/GF/qsUGF2.f90
View File

@ -84,7 +84,6 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
double precision,allocatable :: K(:,:,:)
double precision,allocatable :: SigC(:,:,:)
double precision,allocatable :: SigCp(:,:,:)
double precision,allocatable :: SigCm(:,:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: error(:,:,:)
@ -120,7 +119,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
allocate(eGF2(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), &
SigCm(nBas,nBas,nspin),Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), &
Z(nBas,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), &
F_diis(nBasSq,max_diis,nspin))
! Initialization
@ -191,12 +190,11 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
! Make correlation self-energy Hermitian and transform it back to AO basis
do is=1,nspin
SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is)))
SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is)))
SigC(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is)))
end do
do is=1,nspin
call MOtoAO_transform(nBas,S,c(:,:,is),SigCp(:,:,is))
call MOtoAO_transform(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -326,7 +324,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_f
! Deallocate memory
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis)
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis)
! Perform BSE calculation

10
src/GT/qsGTeh.f90
View File

@ -93,7 +93,6 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
double precision,allocatable :: K(:,:)
double precision,allocatable :: Sig(:,:)
double precision,allocatable :: Sigp(:,:)
double precision,allocatable :: Sigm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
@ -131,7 +130,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -189,10 +188,9 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Make correlation self-energy Hermitian and transform it back to AO basis
Sigp = 0.5d0*(Sig + transpose(Sig))
Sigm = 0.5d0*(Sig - transpose(Sig))
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO_transform(nBas,S,c,Sigp)
call MOtoAO_transform(nBas,S,c,Sig,Sigp)
! Solve the quasi-particle equation
@ -278,7 +276,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,error,error_diis,F_diis)
! Perform BSE calculation

10
src/GT/qsGTpp.f90
View File

@ -92,7 +92,6 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
double precision,allocatable :: K(:,:)
double precision,allocatable :: Sig(:,:)
double precision,allocatable :: Sigp(:,:)
double precision,allocatable :: Sigm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
@ -138,7 +137,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
! Memory allocation
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Sigm(nBas,nBas),Z(nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
@ -233,10 +232,9 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
! Make correlation self-energy Hermitian and transform it back to AO basis
Sigp = 0.5d0*(Sig + transpose(Sig))
Sigm = 0.5d0*(Sig - transpose(Sig))
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO_transform(nBas,S,c,Sigp)
call MOtoAO_transform(nBas,S,c,Sig,Sigp)
! Solve the quasi-particle equation
@ -322,7 +320,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,T
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Sigm,Z,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis)
! Perform BSE calculation

10
src/GT/qsUGTpp.f90
View File

@ -91,7 +91,6 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
double precision,allocatable :: K(:,:,:)
double precision,allocatable :: SigT(:,:,:)
double precision,allocatable :: SigTp(:,:,:)
double precision,allocatable :: SigTm(:,:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: eGT(:,:)
double precision,allocatable :: eOld(:,:)
@ -131,7 +130,7 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Memory allocation
allocate(SigT(nBas,nbas,nspin),SigTp(nBas,nbas,nspin),SigTm(nBas,nbas,nspin), &
allocate(SigT(nBas,nbas,nspin),SigTp(nBas,nbas,nspin), &
Z(nBas,nspin),eGT(nBas,nspin),eOld(nBas,nspin), &
error_diis(nBas,max_diis,nspin),e_diis(nBasSq,max_diis,nspin), &
F_diis(nBasSq,max_diis,nspin),error(nBas,nBas,nspin),&
@ -274,12 +273,11 @@ subroutine qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, &
! Make correlation self-energy Hermitian and transform it back to AO basis
do ispin=1,nspin
SigTp(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) + transpose(SigT(:,:,ispin)))
SigTm(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) - transpose(SigT(:,:,ispin)))
SigT(:,:,ispin) = 0.5d0*(SigT(:,:,ispin) + transpose(SigT(:,:,ispin)))
end do
do ispin=1,nspin
call MOtoAO_transform(nBas,S,c(:,:,ispin),SigTp(:,:,ispin))
call MOtoAO_transform(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin))
end do
! Solve the quasi-particle equation
@ -409,6 +407,6 @@ write(*,*) 'EcGM', EcGM(1)
Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, &
Om1bb,X1bb,Y1bb,Om2bb,X2bb,Y2bb,rho1bb,rho2bb)
deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,SigTm,Z,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,SigT,SigTp,Z,error,error_diis,F_diis)
end subroutine

36
src/GW/SRG_qsGW.f90
View File

@ -77,10 +77,10 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
double precision,allocatable :: XmY_RPA(:,:)
double precision,allocatable :: rho_RPA(:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
double precision,allocatable :: XmY(:,:)
double precision,allocatable :: rho(:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGW(:)
@ -91,6 +91,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
@ -127,9 +128,9 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
! Memory allocation
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), &
rho_RPA(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -169,9 +170,8 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
call wall_time(tao1)
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
@ -187,40 +187,39 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call wall_time(tlr2)
tlr = tlr + tlr2 -tlr1
if(print_W) call print_excitation_energies('phRPA@SRG-qsGW',ispin,nS,OmRPA)
if(print_W) call print_excitation_energies('phRPA@SRG-qsGW',ispin,nS,Om)
! Compute correlation part of the self-energy
call wall_time(tex1)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
call wall_time(tex2)
tex=tex+tex2-tex1
call wall_time(tsrg1)
call SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z)
call SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call wall_time(tsrg2)
tsrg = tsrg + tsrg2 -tsrg1
! Make correlation self-energy Hermitian and transform it back to AO basis
call wall_time(tmo1)
call MOtoAO_transform(nBas,S,c,SigC)
call MOtoAO_transform(nBas,S,c,SigC,SigCp)
call wall_time(tmo2)
tmo = tmo + tmo2 - tmo1
! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigC(:,:)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
! Compute commutator and convergence criteria
@ -241,7 +240,8 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
SigC = matmul(transpose(c),matmul(SigC,c))
call AOtoMO_transform(nBas,c,SigCp,SigC)
! Compute new density matrix in the AO basis
@ -309,7 +309,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
! Perform BSE calculation

13
src/GW/print_qsGGW.f90
View File

@ -29,6 +29,7 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
! Local variables
logical :: dump_orb = .false.
integer :: p,ixyz,HOMO,LUMO
double precision :: Gap
double precision,external :: trace_matrix
@ -103,11 +104,13 @@ subroutine print_qsGGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
write(*,'(A50)') '-----------------------------------------'
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGGW orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
write(*,*)
if(dump_orb) then
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGGW orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGGW orbital energies'
write(*,'(A50)') '---------------------------------------'

13
src/GW/print_qsGW.f90
View File

@ -29,6 +29,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex
! Local variables
logical :: dump_orb = .false.
integer :: p,ixyz,HOMO,LUMO
double precision :: Gap
double precision,external :: trace_matrix
@ -103,11 +104,13 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex
write(*,'(A50)') '-----------------------------------------'
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGW MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
write(*,*)
if(dump_orb) then
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGW MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGW MO energies'
write(*,'(A50)') '---------------------------------------'

25
src/GW/print_qsUGW.f90
View File

@ -31,6 +31,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov, &
! Local variables
logical :: dump_orb = .false.
integer :: p
integer :: ispin,ixyz
double precision :: HOMO(nspin)
@ -158,17 +159,19 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov, &
! Print orbitals
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,1))
write(*,*)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,2))
write(*,*)
if(dump_orb) then
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-up orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,1))
write(*,*)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'qsUGW spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,cGW(:,:,2))
write(*,*)
end if
endif
end if
end subroutine

94
src/GW/qsGGW.f90
View File

@ -65,7 +65,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
double precision :: EK,EKaaaa,EKabba,EKbaab,EKbbbb
double precision :: EqsGW
double precision :: EcRPA
double precision :: EcBSE(nspin)
double precision :: EcBSE
double precision :: EcGM
double precision :: Conv
double precision :: rcond
@ -93,7 +93,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
double precision,allocatable :: C(:,:)
double precision,allocatable :: Cp(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eOld(:)
double precision,allocatable :: P(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: H(:,:)
@ -101,15 +100,16 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
double precision,allocatable :: X(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: err(:,:)
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Self-consistent qsGW calculation |'
write(*,*)'************************************************'
write(*,*)'********************************'
write(*,*)'| Generalized qsGW Calculation |'
write(*,*)'********************************'
write(*,*)
! Warning
@ -138,17 +138,15 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
! Memory allocation
allocate(P(nBas2,nBas2),Jaa(nBas,nBas),Jbb(nBas,nBas), &
allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),P(nBas2,nBas2),Jaa(nBas,nBas),Jbb(nBas,nBas), &
Kaa(nBas,nBas),Kab(nBas,nBas),Kba(nBas,nBas),Kbb(nBas,nBas), &
Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), &
Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), &
F(nBas2,nBas2),Fp(nBas2,nBas2),C(nBas2,nBas2),Cp(nBas2,nBas2), &
H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), &
err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis))
allocate(eGW(nBas2),eOld(nBas2),SigC(nBas2,nBas2),Z(nBas2), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas2,nBas2,nS))
err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis), &
eGW(nBas2),SigC(nBas2,nBas2),SigCp(nBas,nBas),Z(nBas2),Aph(nS,nS),Bph(nS,nS), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas2,nBas2,nS))
! Initialization
@ -158,7 +156,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
Conv = 1d0
P(:,:) = PHF(:,:)
eGW(:) = eHF(:)
eOld(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
err_diis(:,:) = 0d0
@ -182,6 +179,15 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
H( 1:nBas , 1:nBas ) = Hc(1:nBas,1:nBas)
H(nBas+1:nBas2,nBas+1:nBas2) = Hc(1:nBas,1:nBas)
! Construct super density matrix
P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO)))
Paa(:,:) = P( 1:nBas , 1:nBas )
Pab(:,:) = P( 1:nBas ,nBas+1:nBas2)
Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
!------------------------------------------------------------------------
! Main loop
!------------------------------------------------------------------------
@ -204,9 +210,23 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
call exchange_matrix_AO_basis(nBas,Pab,ERI_AO,Kba)
call exchange_matrix_AO_basis(nBas,Pbb,ERI_AO,Kbb)
! Build individual Fock matrices
Faa(:,:) = Hc(:,:) + Jaa(:,:) + Jbb(:,:) + Kaa(:,:)
Fab(:,:) = + Kab(:,:)
Fba(:,:) = + Kba(:,:)
Fbb(:,:) = Hc(:,:) + Jbb(:,:) + Jaa(:,:) + Kbb(:,:)
! Build super Fock matrix
F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas)
F( 1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas)
F(nBas+1:nBas2, 1:nBas ) = Fba(1:nBas,1:nBas)
F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas)
! AO to MO transformation of two-electron integrals
allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2))
allocate(ERI_tmp(nBas2,nBas2,nBas2,nBas2))
Ca(:,:) = C(1:nBas,1:nBas2)
Cb(:,:) = C(nBas+1:nBas2,1:nBas2)
@ -227,7 +247,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp)
ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
deallocate(Ca,Cb,ERI_tmp)
deallocate(ERI_tmp)
! Compute linear response
@ -249,30 +269,18 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
SigC = 0.5d0*(SigC + transpose(SigC))
! call MOtoAO_transform(nBas2,S,C,SigC)
! Build individual Fock matrices
Faa(:,:) = Hc(:,:) + Jaa(:,:) + Jbb(:,:) + Kaa(:,:)
Fab(:,:) = + Kab(:,:)
Fba(:,:) = + Kba(:,:)
Fbb(:,:) = Hc(:,:) + Jbb(:,:) + Jaa(:,:) + Kbb(:,:)
! Build super Fock matrix
F( 1:nBas , 1:nBas ) = Faa(1:nBas,1:nBas)
F( 1:nBas ,nBas+1:nBas2) = Fab(1:nBas,1:nBas)
F(nBas+1:nBas2, 1:nBas ) = Fba(1:nBas,1:nBas)
F(nBas+1:nBas2,nBas+1:nBas2) = Fbb(1:nBas,1:nBas)
call MOtoAO_transform_GHF(nBas2,nBas,S,Ca,Cb,SigC,SigCp)
! ... and add self-energy
F(:,:) = F(:,:) + SigC(:,:)
F(:,:) = F(:,:) + SigCp(:,:)
! Compute commutator and convergence criteria
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
if(nSCF > 1) Conv = maxval(abs(err))
! DIIS extrapolation
if(max_diis > 1) then
@ -280,6 +288,10 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis,F_diis,err,F)
end if
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X),matmul(F,X))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
Cp(:,:) = Fp(:,:)
@ -289,7 +301,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
C(:,:) = matmul(X,Cp)
SigC = matmul(transpose(c),matmul(SigC,c))
call AOtoMO_transform_GHF(nBas,nBas2,Ca,Cb,SigCp,SigC)
! Form super density matrix
@ -302,11 +314,6 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
! Save quasiparticles energy for next cycle
Conv = maxval(abs(err))
eOld(:) = eGW(:)
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
@ -350,7 +357,7 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
! Print results
call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsGW(nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
call print_qsGGW(nBas2,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
enddo
!------------------------------------------------------------------------
@ -377,19 +384,10 @@ subroutine qsGGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,do
call GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas2,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
if(exchange_kernel) then
EcBSE(1) = 0.5d0*EcBSE(1)
EcBSE(2) = 1.5d0*EcBSE(2)
end if
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW correlation energy =',EcBSE,' au'
write(*,'(2X,A50,F20.10,A4)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE,' au'
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

24
src/GW/qsGW.f90
View File

@ -85,7 +85,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eOld(:)
double precision,allocatable :: P(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
@ -93,7 +92,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
double precision,allocatable :: K(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: SigCm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error(:,:)
@ -130,9 +128,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Memory allocation
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization
@ -143,7 +141,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
Conv = 1d0
P(:,:) = PHF(:,:)
eGW(:) = eHF(:)
eOld(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
@ -169,9 +166,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! AO to MO transformation of two-electron integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
@ -194,10 +190,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Make correlation self-energy Hermitian and transform it back to AO basis
SigCp = 0.5d0*(SigC + transpose(SigC))
SigCm = 0.5d0*(SigC - transpose(SigC))
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO_transform(nBas,S,c,SigCp)
call MOtoAO_transform(nBas,S,c,SigC,SigCp)
! Solve the quasi-particle equation
@ -222,7 +217,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
SigCp = matmul(transpose(c),matmul(SigCp,c))
call AOtoMO_transform(nBas,c,SigCp,SigC)
! Compute new density matrix in the AO basis
@ -230,8 +225,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Save quasiparticles energy for next cycle
Conv = maxval(abs(error))
eOld(:) = eGW(:)
if(nSCF > 1) Conv = maxval(abs(error))
!------------------------------------------------------------------------
! Compute total energy
@ -283,7 +277,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
! Perform BSE calculation

70
src/GW/qsUGW.f90
View File

@ -69,7 +69,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision :: ET(nspin)
double precision :: EV(nspin)
double precision :: EJ(nsp)
double precision :: Ex(nspin)
double precision :: EK(nspin)
double precision :: EcRPA
double precision :: EcGM(nspin)
double precision :: EqsGW
@ -78,7 +78,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision :: Conv
double precision :: rcond(nspin)
double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:,:)
double precision,allocatable :: err_diis(:,:,:)
double precision,allocatable :: F_diis(:,:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: XpY(:,:)
@ -87,7 +87,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision,allocatable :: c(:,:,:)
double precision,allocatable :: cp(:,:,:)
double precision,allocatable :: eGW(:,:)
double precision,allocatable :: eOld(:,:)
double precision,allocatable :: P(:,:,:)
double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:,:)
@ -95,9 +94,8 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
double precision,allocatable :: K(:,:,:)
double precision,allocatable :: SigC(:,:,:)
double precision,allocatable :: SigCp(:,:,:)
double precision,allocatable :: SigCm(:,:,:)
double precision,allocatable :: Z(:,:)
double precision,allocatable :: error(:,:,:)
double precision,allocatable :: err(:,:,:)
! Hello world
@ -136,25 +134,24 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(eGW(nBas,nspin),eOld(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), &
SigCm(nBas,nBas,nspin),Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc), &
rho(nBas,nBas,nS_sc,nspin),error(nBas,nBas,nspin),error_diis(nBasSq,max_diis,nspin), &
allocate(eGW(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin),P(nBas,nBas,nspin),F(nBas,nBas,nspin), &
Fp(nBas,nBas,nspin),J(nBas,nBas,nspin),K(nBas,nBas,nspin),SigC(nBas,nBas,nspin),SigCp(nBas,nBas,nspin), &
Z(nBas,nspin),Om(nS_sc),XpY(nS_sc,nS_sc),XmY(nS_sc,nS_sc), &
rho(nBas,nBas,nS_sc,nspin),err(nBas,nBas,nspin),err_diis(nBasSq,max_diis,nspin), &
F_diis(nBasSq,max_diis,nspin))
! Initialization
nSCF = -1
n_diis = 0
ispin = 1
Conv = 1d0
P(:,:,:) = PHF(:,:,:)
eGW(:,:) = eHF(:,:)
eOld(:,:) = eHF(:,:)
c(:,:,:) = cHF(:,:,:)
F_diis(:,:,:) = 0d0
error_diis(:,:,:) = 0d0
rcond(:) = 0d0
nSCF = -1
n_diis = 0
ispin = 1
Conv = 1d0
P(:,:,:) = PHF(:,:,:)
eGW(:,:) = eHF(:,:)
c(:,:,:) = cHF(:,:,:)
F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0
rcond(:) = 0d0
!------------------------------------------------------------------------
! Main loop
@ -182,12 +179,9 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! AO to MO transformation of two-electron integrals
!--------------------------------------------------
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -228,12 +222,11 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Make correlation self-energy Hermitian and transform it back to AO basis
do is=1,nspin
SigCp(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is)))
SigCm(:,:,is) = 0.5d0*(SigC(:,:,is) - transpose(SigC(:,:,is)))
SigC(:,:,is) = 0.5d0*(SigC(:,:,is) + transpose(SigC(:,:,is)))
end do
do is=1,nspin
call MOtoAO_transform(nBas,S,c(:,:,is),SigCp(:,:,is))
call MOtoAO_transform(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -245,18 +238,18 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Check convergence
do is=1,nspin
error(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is))
err(:,:,is) = matmul(F(:,:,is),matmul(P(:,:,is),S(:,:))) - matmul(matmul(S(:,:),P(:,:,is)),F(:,:,is))
end do
if(nSCF > 1) conv = maxval(abs(error(:,:,:)))
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(minval(rcond(:)) > 1d-7) then
do is=1,nspin
if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,error_diis(:,1:n_diis,is), &
F_diis(:,1:n_diis,is),error(:,:,is),F(:,:,is))
if(nO(is) > 1) call DIIS_extrapolation(rcond(is),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,is), &
F_diis(:,1:n_diis,is),err(:,:,is),F(:,:,is))
end do
else
n_diis = 0
@ -293,11 +286,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
P(:,:,is) = matmul(c(:,1:nO(is),is),transpose(c(:,1:nO(is),is)))
end do
! Save quasiparticles energy for next cycle
Conv = maxval(abs(eGW(:,:) - eOld(:,:)))
eOld(:,:) = eGW(:,:)
!------------------------------------------------------------------------
! Compute total energy
!------------------------------------------------------------------------
@ -324,19 +312,19 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Exchange energy
do is=1,nspin
Ex(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is)))
EK(is) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,is),K(:,:,is)))
end do
! Total energy
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
EqsGW = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(EK(:))
!------------------------------------------------------------------------
! Print results
!------------------------------------------------------------------------
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,SigCp,Z,dipole)
enddo
!------------------------------------------------------------------------
@ -359,7 +347,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_W,
! Deallocate memory
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis)
! Perform BSE calculation

6
src/HF/GHF.f90
View File

@ -72,9 +72,9 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
! Hello world
write(*,*)
write(*,*)'***********************************************'
write(*,*)'* Generalized Hartree-Fock calculation *'
write(*,*)'***********************************************'
write(*,*)'****************************************'
write(*,*)'* Generalized Hartree-Fock Calculation *'
write(*,*)'****************************************'
write(*,*)
! Useful stuff

4
src/HF/ROHF_fock_matrix.f90
View File

@ -59,8 +59,8 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,Fa,Fb,F)
! Block-by-block Fock matrix
call AOtoMO_transform(nBas,c,Fa)
call AOtoMO_transform(nBas,c,Fb)
Fa = matmul(transpose(c),matmul(Fa,c))
Fb = matmul(transpose(c),matmul(Fb,c))
F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC )
F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO )

10
src/HF/UHF.f90
View File

@ -35,7 +35,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
integer :: nSCF
integer :: nBasSq
integer :: n_diis
double precision :: conv
double precision :: Conv
double precision :: rcond(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
@ -90,7 +90,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
! Initialization
nSCF = 0
conv = 1d0
Conv = 1d0
n_diis = 0
F_diis(:,:,:) = 0d0
@ -106,7 +106,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
'|','#','|','E(UHF)','|','Ex(UHF)','|','Conv','|'
write(*,*)'----------------------------------------------------------'
do while(conv > thresh .and. nSCF < maxSCF)
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
@ -136,7 +136,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin))
end do
if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
! DIIS extrapolation
@ -224,7 +224,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',conv,'|'
'|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|'
end do
write(*,*)'----------------------------------------------------------'

8
src/HF/print_GHF.f90
View File

@ -61,7 +61,7 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
Ca(:,:) = C( 1:nBas ,1:nBas2)
Cb(:,:) = C(nBas+1:nBas2,1:nBas2)
! Compute expectation values of S^2
! Compute expectation values of S^2 (WRONG!)
Sx2 = 0.25d0*trace_matrix(nBas,Paa+Pbb) + 0.25d0*trace_matrix(nBas,Pab+Pba)**2
do mu=1,nBas
@ -109,8 +109,8 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
write(*,'(A32,1X,F16.6,A3)') ' GHF LUMO energy: ',e(LUMO)*HaToeV,' eV'
write(*,'(A32,1X,F16.6,A3)') ' GHF HOMO-LUMO gap : ',Gap*HaToeV,' eV'
write(*,'(A50)') '-----------------------------------------'
write(*,'(A32,1X,F16.6)') ' <S**2> :',S2
write(*,'(A50)') '-----------------------------------------'
! write(*,'(A32,1X,F16.6)') ' <S**2> :',S2
! write(*,'(A50)') '-----------------------------------------'
write(*,'(A35)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
@ -124,8 +124,8 @@ subroutine print_GHF(nBas,nBas2,nO,e,C,P,ENuc,ET,EV,EJ,EK,EHF,dipole)
write(*,'(A32)') ' GHF orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas2,nBas2,c)
write(*,*)
end if
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' GHF orbital energies'
write(*,'(A50)') '---------------------------------------'

2
src/HF/print_RHF.f90
View File

@ -76,8 +76,8 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
write(*,'(A32)') ' RHF orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cHF)
write(*,*)
end if
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' RHF orbital energies'
write(*,'(A50)') '---------------------------------------'

2
src/HF/print_UHF.f90
View File

@ -119,8 +119,8 @@ subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
write(*,'(A50)') 'UHF spin-down orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c(:,:,2))
write(*,*)
end if
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' UHF spin-up orbital energies '
write(*,'(A50)') '---------------------------------------'

3
src/QuAcK/RQuAcK.f90
View File

@ -150,9 +150,8 @@ subroutine RQuAcK(doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Read and transform dipole-related integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
call AOtoMO_transform(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
! 4-index transform

7
src/QuAcK/UQuAcK.f90
View File

@ -157,12 +157,9 @@ subroutine UQuAcK(doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,do
! Read and transform dipole-related integrals
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block