From 01b0a4d823ec1187ccce77ea792d9ba11ae3e1bd Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Tue, 14 Mar 2023 14:11:01 +0100 Subject: [PATCH] add SRG files --- src/GW/SRG_qsGW.f90 | 53 +++++++++++++++--- src/GW/self_energy_correlation_SRG.f90 | 74 +++++++++++++++++--------- src/QuAcK/QuAcK.f90 | 36 ++++++------- 3 files changed, 112 insertions(+), 51 deletions(-) diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 76009a8..293217f 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -68,6 +68,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE double precision :: EcGM double precision :: Conv 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 :: dipole(ncart) @@ -155,34 +156,63 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE nSCF = nSCF + 1 ! Buid Coulomb matrix - + call wall_time(t1) call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy 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 - 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 - 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(tlr1) + + 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) ! 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 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 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 - + + call wall_time(tmo1) call MOtoAO_transform(nBas,S,c,SigC) - - ! Solve the quasi-particle equation + call wall_time(tmo2) + tmo = tmo + tmo2 - tmo1 + ! Solve the quasi-particle equation 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 - 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 diff --git a/src/GW/self_energy_correlation_SRG.f90 b/src/GW/self_energy_correlation_SRG.f90 index 1090662..2910269 100644 --- a/src/GW/self_energy_correlation_SRG.f90 +++ b/src/GW/self_energy_correlation_SRG.f90 @@ -23,8 +23,9 @@ subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM, integer :: i,j,a,b integer :: p,q,r integer :: m - double precision :: Dpim,Dqim,Dpam,Dqam - + double precision :: Dpim,Dqim,Dpam,Dqam + double precision :: t1,t2 + ! Output variables double precision,intent(out) :: EcGM @@ -38,36 +39,59 @@ subroutine self_energy_correlation_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM, ! SRG-GW self-energy ! !--------------------! -! Occupied part of the correlation self-energy + ! Occupied part of the correlation self-energy - - do m=1,nS - do i=nC+1,nO - do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - 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)*(1-exp(-eta*Dpim**2)*exp(-eta*Dqim**2)) & - *(Dpim + Dqim)/(Dpim**2 + Dqim**2) + call wall_time(t1) + + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & + !$OMP PRIVATE(m,i,q,p,Dpim,Dqim) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + 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 + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, "first loop", (t2-t1) ! Virtual part of the correlation self-energy - do m=1,nS - do a=nO+1,nBas-nR - do q=nC+1,nBas-nR - do p=nC+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)*(1-exp(-eta*Dpam**2)*exp(-eta*Dqam**2)) & - *(Dpam + Dqam)/(Dpam**2 + Dqam**2) - end do - end do + call wall_time(t1) + !$OMP PARALLEL & + !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nR,nBas,e,Omega) & + !$OMP PRIVATE(m,a,q,p,Dpam,Dqam) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do m=1,nS + 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 + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, "second loop", (t2-t1) + ! Galitskii-Migdal correlation energy diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 2c81dc4..756d59c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -155,7 +155,7 @@ program QuAcK doSph = .false. - call cpu_time(start_QuAcK) + call wall_time(start_QuAcK) ! Which calculations do you want to do? @@ -228,7 +228,7 @@ program QuAcK ! Read integrals - call cpu_time(start_int) + call wall_time(start_int) if(doSph) then @@ -241,11 +241,11 @@ program QuAcK end if - call cpu_time(end_int) + call wall_time(end_int) t_int = end_int - start_int 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(*,*) ! Compute orthogonalization matrix @@ -266,13 +266,13 @@ program QuAcK stop 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, & 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 - 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(*,*) end if @@ -350,7 +350,7 @@ program QuAcK ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ - call cpu_time(start_AOtoMO) + call wall_time(start_AOtoMO) write(*,*) write(*,*) 'AO to MO transformation... Please be patient' @@ -434,10 +434,10 @@ program QuAcK end if - call cpu_time(end_AOtoMO) + call wall_time(end_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(*,*) !------------------------------------------------------------------------ @@ -1032,7 +1032,7 @@ program QuAcK if(doqsGW) then - call cpu_time(start_qsGW) + call wall_time(start_qsGW) if(unrestricted) then @@ -1049,10 +1049,10 @@ program QuAcK end if - call cpu_time(end_qsGW) + call wall_time(end_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(*,*) end if @@ -1063,7 +1063,7 @@ program QuAcK if(doSRGqsGW) then - call cpu_time(start_qsGW) + call wall_time(start_qsGW) if(unrestricted) then @@ -1077,10 +1077,10 @@ program QuAcK end if - call cpu_time(end_qsGW) + call wall_time(end_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(*,*) end if @@ -1240,10 +1240,10 @@ program QuAcK ! End of QuAcK !------------------------------------------------------------------------ - call cpu_time(end_QuAcK) + call wall_time(end_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(*,*) end program QuAcK