LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 12:56:38 +01:00
Code Releases Activity

openMP for GW

Browse Source
This commit is contained in:
Antoine Marie 2023-03-14 14:11:54 +01:00
parent 01b0a4d823
commit 3a80069009
2 changed files with 44 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/GW/self_energy_correlation.f90
View File

@ -78,29 +78,43 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
do jb=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(jb,i,q,p,eps) &
do p=nC+1,nBas-nR !$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do i=nC+1,nO
eps = e(p) - e(i) + Omega(jb) eps = e(p) - e(i) + Omega(jb)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2) SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*eps/(eps**2 + eta**2)
end do end do
end do
end do end do
end do end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Virtual part of the correlation self-energy ! Virtual part of the correlation self-energy
do jb=1,nS !$OMP PARALLEL &
do a=nO+1,nBas-nR !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Omega) &
do q=nC+1,nBas-nR !$OMP PRIVATE(jb,a,q,p,eps) &
do p=nC+1,nBas-nR !$OMP DEFAULT(NONE)
!$OMP DO
do q=nC+1,nBas-nR
do p=nC+1,nBas-nR
do jb=1,nS
do a=nO+1,nBas-nR
eps = e(p) - e(a) - Omega(jb) eps = e(p) - e(a) - Omega(jb)
SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2) SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*eps/(eps**2 + eta**2)
end do end do
end do
end do end do
end do end do
end do
!$OMP END DO
!$OMP END PARALLEL
! Galitskii-Migdal correlation energy ! Galitskii-Migdal correlation energy

24
src/LR/linear_response.f90
View File

@ -21,7 +21,9 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
double precision,intent(in) :: e(nBas) double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
integer :: i,j,k
double precision :: trace_matrix double precision :: trace_matrix
double precision,allocatable :: A(:,:) double precision,allocatable :: A(:,:)
@ -31,6 +33,7 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
double precision,allocatable :: AmBSq(:,:) double precision,allocatable :: AmBSq(:,:)
double precision,allocatable :: AmBIv(:,:) double precision,allocatable :: AmBIv(:,:)
double precision,allocatable :: Z(:,:) double precision,allocatable :: Z(:,:)
double precision,allocatable :: tmp(:,:)
! Output variables ! Output variables
@ -41,7 +44,7 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
! Memory allocation ! Memory allocation
allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS)) allocate(A(nS,nS),B(nS,nS),ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS))
! Build A and B matrices ! Build A and B matrices
@ -77,10 +80,11 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
! if(Omega(ia) < 0d0) Omega(ia) = 0d0 ! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do ! end do
call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) call ADAt(nS,AmB,1d0*dsqrt(Omega),AmBSq)
call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) call ADAt(nS,AmB,1d0/dsqrt(Omega),AmBIv)
Z = matmul(AmBSq,matmul(ApB,AmBSq)) call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
call diagonalize_matrix(nS,Z,Omega) call diagonalize_matrix(nS,Z,Omega)
@ -91,13 +95,13 @@ subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,E
! if(Omega(ia) < 0d0) Omega(ia) = 0d0 ! if(Omega(ia) < 0d0) Omega(ia) = 0d0
! end do ! end do
Omega = sqrt(Omega) Omega = dsqrt(Omega)
XpY = matmul(transpose(Z),AmBSq) call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1))
call DA(nS,1d0/sqrt(Omega),XpY) call DA(nS,1d0/dsqrt(Omega),XpY)
XmY = matmul(transpose(Z),AmBIv) call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
call DA(nS,1d0*sqrt(Omega),XmY) call DA(nS,1d0*dsqrt(Omega),XmY)
end if end if