10
1
mirror of https://github.com/pfloos/quack synced 2024-12-23 12:55:25 +01:00

add SRG files

This commit is contained in:
Antoine Marie 2023-03-14 14:11:01 +01:00
parent 616d6747e6
commit 01b0a4d823
3 changed files with 112 additions and 51 deletions

View File

@ -68,6 +68,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
double precision :: EcGM double precision :: EcGM
double precision :: Conv double precision :: Conv
double precision :: rcond double precision :: rcond
double precision :: tao,tao1,tao2,tsrg,tsrg1,tsrg2,tlr,tlr1,tlr2,t1,t2,tt,tmo1,tmo2,tmo,tex,tex1,tex2
double precision,external :: trace_matrix double precision,external :: trace_matrix
double precision :: dipole(ncart) double precision :: dipole(ncart)
@ -155,34 +156,63 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
nSCF = nSCF + 1 nSCF = nSCF + 1
! Buid Coulomb matrix ! Buid Coulomb matrix
call wall_time(t1)
call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J)
! Compute exchange part of the self-energy ! Compute exchange part of the self-energy
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
call wall_time(t2)
tt=tt+t2-t1
! AO to MO transformation of two-electron integrals ! AO to MO transformation of two-electron integrals
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) call wall_time(tao1)
call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO)
call wall_time(tao2)
tao = tao + tao2 -tao1
! Compute linear response ! Compute linear response
call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & call wall_time(tlr1)
EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, &
EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call wall_time(tlr2)
tlr = tlr + tlr2 -tlr1
if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA)
! Compute correlation part of the self-energy ! Compute correlation part of the self-energy
call wall_time(tex1)
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA)
call wall_time(tex2)
tex=tex+tex2-tex1
call wall_time(tsrg1)
call self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC)
call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z)
call wall_time(tsrg2)
tsrg = tsrg + tsrg2 -tsrg1
! Make correlation self-energy Hermitian and transform it back to AO basis ! 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)
call wall_time(tmo2)
! Solve the quasi-particle equation tmo = tmo + tmo2 - tmo1
! Solve the quasi-particle equation
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigC(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigC(:,:)
@ -262,7 +292,14 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE
stop stop
endif endif
print *, "Wall time for Fock and exchange build", tt
print *, "Wall Time for AO to MO", tao
print *, "Wall Time for LR", tlr
print *, "Wall Time for excitation density", tex
print *, "Wall Time for SRG", tsrg
print *, "Wall time MO to AO Sigma", tmo
! Deallocate memory ! Deallocate memory

View File

@ -23,7 +23,8 @@ subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,
integer :: i,j,a,b integer :: i,j,a,b
integer :: p,q,r integer :: p,q,r
integer :: m integer :: m
double precision :: Dpim,Dqim,Dpam,Dqam double precision :: Dpim,Dqim,Dpam,Dqam
double precision :: t1,t2
! Output variables ! Output variables
@ -38,36 +39,59 @@ subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,
! SRG-GW self-energy ! ! SRG-GW self-energy !
!--------------------! !--------------------!
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
call wall_time(t1)
do m=1,nS !$OMP PARALLEL &
do i=nC+1,nO !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
do q=nC+1,nBas-nR !$OMP PRIVATE(m,i,q,p,Dpim,Dqim) &
do p=nC+1,nBas-nR !$OMP DEFAULT(NONE)
Dpim = e(p) - e(i) + Omega(m) !$OMP DO
Dqim = e(q) - e(i) + Omega(m) do q=nC+1,nBas-nR
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*(1-exp(-eta*Dpim**2)*exp(-eta*Dqim**2)) & do p=nC+1,nBas-nR
*(Dpim + Dqim)/(Dpim**2 + Dqim**2) do m=1,nS
do i=nC+1,nO
Dpim = e(p) - e(i) + Omega(m)
Dqim = e(q) - e(i) + Omega(m)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*(1d0-dexp(-eta*Dpim*Dpim)*dexp(-eta*Dqim*Dqim)) &
*(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim)
end do
end do end do
end do end do
end do
end do end do
!$OMP END DO
!$OMP END PARALLEL
call wall_time(t2)
print *, "first loop", (t2-t1)
! Virtual part of the correlation self-energy ! Virtual part of the correlation self-energy
do m=1,nS call wall_time(t1)
do a=nO+1,nBas-nR !$OMP PARALLEL &
do q=nC+1,nBas-nR !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nR,nBas,e,Omega) &
do p=nC+1,nBas-nR !$OMP PRIVATE(m,a,q,p,Dpam,Dqam) &
Dpam = e(p) - e(a) - Omega(m) !$OMP DEFAULT(NONE)
Dqam = e(q) - e(a) - Omega(m) !$OMP DO
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*(1-exp(-eta*Dpam**2)*exp(-eta*Dqam**2)) & do q=nC+1,nBas-nR
*(Dpam + Dqam)/(Dpam**2 + Dqam**2) do p=nC+1,nBas-nR
end do do m=1,nS
end do do a=nO+1,nBas-nR
Dpam = e(p) - e(a) - Omega(m)
Dqam = e(q) - e(a)- Omega(m)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*(1d0-exp(-eta*Dpam*Dpam)*exp(-eta*Dqam*Dqam)) &
*(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam)
end do
end do
end do end do
end do end do
!$OMP END DO
!$OMP END PARALLEL
call wall_time(t2)
print *, "second loop", (t2-t1)
! Galitskii-Migdal correlation energy ! Galitskii-Migdal correlation energy

View File

@ -155,7 +155,7 @@ program QuAcK
doSph = .false. doSph = .false.
call cpu_time(start_QuAcK) call wall_time(start_QuAcK)
! Which calculations do you want to do? ! Which calculations do you want to do?
@ -228,7 +228,7 @@ program QuAcK
! Read integrals ! Read integrals
call cpu_time(start_int) call wall_time(start_int)
if(doSph) then if(doSph) then
@ -241,11 +241,11 @@ program QuAcK
end if end if
call cpu_time(end_int) call wall_time(end_int)
t_int = end_int - start_int t_int = end_int - start_int
write(*,*) write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds'
write(*,*) write(*,*)
! Compute orthogonalization matrix ! Compute orthogonalization matrix
@ -266,13 +266,13 @@ program QuAcK
stop stop
end if end if
call cpu_time(start_HF) call wall_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*) write(*,*)
end if end if
@ -350,7 +350,7 @@ program QuAcK
! AO to MO integral transform for post-HF methods ! AO to MO integral transform for post-HF methods
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call cpu_time(start_AOtoMO) call wall_time(start_AOtoMO)
write(*,*) write(*,*)
write(*,*) 'AO to MO transformation... Please be patient' write(*,*) 'AO to MO transformation... Please be patient'
@ -434,10 +434,10 @@ program QuAcK
end if end if
call cpu_time(end_AOtoMO) call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*) write(*,*)
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -1032,7 +1032,7 @@ program QuAcK
if(doqsGW) then if(doqsGW) then
call cpu_time(start_qsGW) call wall_time(start_qsGW)
if(unrestricted) then if(unrestricted) then
@ -1049,10 +1049,10 @@ program QuAcK
end if end if
call cpu_time(end_qsGW) call wall_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW t_qsGW = end_qsGW - start_qsGW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_qsGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1063,7 +1063,7 @@ program QuAcK
if(doSRGqsGW) then if(doSRGqsGW) then
call cpu_time(start_qsGW) call wall_time(start_qsGW)
if(unrestricted) then if(unrestricted) then
@ -1077,10 +1077,10 @@ program QuAcK
end if end if
call cpu_time(end_qsGW) call wall_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW t_qsGW = end_qsGW - start_qsGW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_qsGW,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds'
write(*,*) write(*,*)
end if end if
@ -1240,10 +1240,10 @@ program QuAcK
! End of QuAcK ! End of QuAcK
!------------------------------------------------------------------------ !------------------------------------------------------------------------
call cpu_time(end_QuAcK) call wall_time(end_QuAcK)
t_QuAcK = end_QuAcK - start_QuAcK t_QuAcK = end_QuAcK - start_QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for QuAcK = ',t_QuAcK,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds'
write(*,*) write(*,*)
end program QuAcK end program QuAcK