PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-05-06 23:24:58 +02:00
Code Releases Activity

bug fixes and rm error in linear solve

Browse Source
This commit is contained in:
Loris Burth 2025-04-17 09:56:12 +02:00
parent cdc1c25351
commit f84da77244
4 changed files with 47 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/AOtoMO/complex_MOtoAO.f90
View File

@ -17,7 +17,7 @@ subroutine complex_MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs)
allocate(SC(nBas,nOrb), BSC(nOrb,nBas),cS(nBas,nBas)) allocate(SC(nBas,nOrb), BSC(nOrb,nBas),cS(nBas,nBas))
cS(:,:) = (0d0,1d0)*S(:,:) cS(:,:) = (1d0,0d0)*S(:,:)
!SC = matmul(S, C) !SC = matmul(S, C)
!BSC = matmul(M_MOs, transpose(SC)) !BSC = matmul(M_MOs, transpose(SC))
!M_AOs = matmul(SC, BSC) !M_AOs = matmul(SC, BSC)

64
src/GW/complex_RGW_self_energy.f90
View File

@ -1,4 +1,4 @@
subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z) subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Sig,Z)
! Compute correlation part of the self-energy and the renormalization factor ! Compute correlation part of the self-energy and the renormalization factor
@ -15,8 +15,7 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
integer,intent(in) :: nV integer,intent(in) :: nV
integer,intent(in) :: nR integer,intent(in) :: nR
integer,intent(in) :: nS integer,intent(in) :: nS
double precision,intent(in) :: Re_e(nOrb) complex*16,intent(in) :: e(nOrb)
double precision,intent(in) :: Im_e(nOrb)
complex*16,intent(in) :: Om(nS) complex*16,intent(in) :: Om(nS)
complex*16,intent(in) :: rho(nOrb,nOrb,nS) complex*16,intent(in) :: rho(nOrb,nOrb,nS)
@ -26,29 +25,33 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
integer :: p,q,m integer :: p,q,m
double precision :: eps,eta_tilde double precision :: eps,eta_tilde
complex*16 :: num,tmp complex*16 :: num,tmp
double precision, allocatable :: Re_DS(:) double precision,allocatable :: Re_DS(:)
double precision, allocatable :: Im_DS(:) double precision,allocatable :: Im_DS(:)
double precision,allocatable :: Re_Sig(:,:)
double precision,allocatable :: Im_Sig(:,:)
double precision,allocatable :: Re_Z(:)
double precision,allocatable :: Im_Z(:)
! Output variables ! Output variables
complex*16,intent(out) :: EcGM complex*16,intent(out) :: EcGM
double precision,intent(out) :: Re_Sig(nOrb,nOrb) complex*16,intent(out) :: Sig(nOrb,nOrb)
double precision,intent(out) :: Im_Sig(nOrb,nOrb) complex*16,intent(out) :: Z(nOrb)
double precision,intent(out) :: Re_Z(nOrb)
double precision,intent(out) :: Im_Z(nOrb)
!----------------! !----------------!
! GW self-energy ! ! GW self-energy !
!----------------! !----------------!
allocate(Re_DS(nBas),Im_DS(nBas)) allocate(Re_DS(nBas),Im_DS(nBas),Re_Z(nOrb),Im_Z(nOrb),Re_Sig(nOrb,nOrb),Im_Sig(nOrb,nOrb))
Re_Sig(:,:) = 0d0 Re_Sig(:,:) = 0d0
Im_Sig(:,:) = 0d0 Im_Sig(:,:) = 0d0
Re_DS(:) = 0d0
Im_DS(:) = 0d0
! Occupied part of the correlation self-energy ! Occupied part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) & !$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,i,q,p,eps,num,eta_tilde,tmp) & !$OMP PRIVATE(m,i,q,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
@ -57,8 +60,8 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
do m=1,nS do m=1,nS
do i=nC+1,nO do i=nC+1,nO
eps = Re_e(p) - Re_e(i) + real(Om(m)) eps = real(e(p)) - real(e(i)) + real(Om(m))
eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m))
num = 2d0*rho(p,i,m)*rho(q,i,m) num = 2d0*rho(p,i,m)*rho(q,i,m)
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
eta_tilde/(eps**2+eta_tilde**2),kind=8) eta_tilde/(eps**2+eta_tilde**2),kind=8)
@ -75,7 +78,7 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
! Virtual part of the correlation self-energy ! Virtual part of the correlation self-energy
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) & !$OMP SHARED(Re_Sig,Im_Sig,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,a,q,p,eps,num,eta_tilde,tmp) & !$OMP PRIVATE(m,a,q,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
@ -84,8 +87,8 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
do m=1,nS do m=1,nS
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
eps = Re_e(p) - Re_e(a) - real(Om(m)) eps = real(e(p)) - real(e(a)) - real(Om(m))
eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m))
num = 2d0*rho(p,a,m)*rho(q,a,m) num = 2d0*rho(p,a,m)*rho(q,a,m)
tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
-eta_tilde/(eps**2 + eta_tilde**2),kind=8) -eta_tilde/(eps**2 + eta_tilde**2),kind=8)
@ -103,21 +106,18 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
! Renormalization factor ! ! Renormalization factor !
!------------------------! !------------------------!
Re_DS(:) = 0d0
Im_DS(:) = 0d0
! Occupied part of the renormalization factor ! Occupied part of the renormalization factor
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) & !$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) & !$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
do p=nC+1,nOrb-nR do p=nC+1,nOrb-nR
do m=1,nS do m=1,nS
do i=nC+1,nO do i=nC+1,nO
eps = Re_e(p) - Re_e(i) + real(Om(m)) eps = real(e(p)) - real(e(i)) + real(Om(m))
eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m))
num = 2d0*rho(p,i,m)*rho(p,i,m) num = 2d0*rho(p,i,m)*rho(p,i,m)
tmp = num*cmplx(-(eps**2-eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& tmp = num*cmplx(-(eps**2-eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
-2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8) -2*eta_tilde*eps/(eps**2 + eta_tilde**2)**2,kind=8)
@ -132,7 +132,7 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
! Virtual part of the renormalization factor ! Virtual part of the renormalization factor
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) & !$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,e,Om) &
!$OMP PRIVATE(m,a,p,eps,num,eta_tilde,tmp) & !$OMP PRIVATE(m,a,p,eps,num,eta_tilde,tmp) &
!$OMP DEFAULT(NONE) !$OMP DEFAULT(NONE)
!$OMP DO !$OMP DO
@ -140,8 +140,8 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
do m=1,nS do m=1,nS
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
eps = Re_e(p) - Re_e(a) - real(Om(m)) eps = real(e(p)) - real(e(a)) - real(Om(m))
eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m))
num = 2d0*rho(p,a,m)*rho(p,a,m) num = 2d0*rho(p,a,m)*rho(p,a,m)
tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,&
2*eta_tilde*eps/eps/(eps**2 + eta_tilde**2)**2,kind=8) 2*eta_tilde*eps/eps/(eps**2 + eta_tilde**2)**2,kind=8)
@ -156,7 +156,11 @@ subroutine complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho
! Compute renormalization factor from derivative ! Compute renormalization factor from derivative
Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) Re_Z(:) = (1d0-Re_DS(:))/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2) Im_Z(:) = Im_DS(:)/((1d0 - Re_DS(:))**2 + Im_DS(:)**2)
deallocate(Re_DS,Im_DS)
Z = cmplx(Re_Z,Im_Z,kind=8)
Sig = cmplx(Re_Sig,Im_Sig,kind=8)
deallocate(Re_DS,Im_DS,Re_Z,Im_Z,Re_Sig,Im_Sig)
!!-------------------------------------! !!-------------------------------------!
!! Galitskii-Migdal correlation energy ! !! Galitskii-Migdal correlation energy !

12
src/GW/complex_qsRGW.f90
View File

@ -53,7 +53,7 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: CAP_AO(nBas,nBas) double precision,intent(in) :: CAP_AO(nBas,nBas)
double precision,intent(inout):: CAP_MO(nBas,nBas) complex*16,intent(inout) :: CAP_MO(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
complex*16,intent(inout) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) complex*16,intent(inout) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
@ -218,10 +218,9 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
write(*,*) "SRG not implemented" write(*,*) "SRG not implemented"
!call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) !call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
else else
call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,real(eGW),aimag(eGW),Om,rho,& call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,&
EcGM,real(SigC),aimag(SigC),real(Z),aimag(Z)) EcGM,SigC,Z)
end if end if
! Make correlation self-energy Hermitian and transform it back to AO basis ! Make correlation self-energy Hermitian and transform it back to AO basis
SigC = 0.5d0*(SigC + transpose(SigC)) SigC = 0.5d0*(SigC + transpose(SigC))
@ -235,7 +234,7 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
call complex_complex_AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1)) call complex_complex_AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1))
call complex_MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1)) call complex_MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1))
endif endif
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
@ -281,16 +280,17 @@ subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d
Fp = matmul(transpose(X),matmul(F,X)) Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb,cp,eGW) call complex_diagonalize_matrix(nOrb,cp,eGW)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(X,cp) c = matmul(X,cp)
else else
Fp = matmul(transpose(c),matmul(F,c)) Fp = matmul(transpose(c),matmul(F,c))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call complex_diagonalize_matrix(nOrb,cp,eGW) call complex_diagonalize_matrix(nOrb,cp,eGW)
call complex_orthogonalize_matrix(nBas,cp)
c = matmul(c,cp) c = matmul(c,cp)
endif endif
call complex_complex_AOtoMO(nBas,nOrb,c,SigCp,SigC) call complex_complex_AOtoMO(nBas,nOrb,c,SigCp,SigC)
! Density matrix ! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))

12
src/utils/wrap_lapack.f90
View File

@ -410,12 +410,12 @@ subroutine complex_linear_solve(N,A,b,x,rcond)
call zsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,rwork,info) call zsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,rwork,info)
if (info /= 0) then ! if (info /= 0) then
!
print *, info ! print *, info
stop 'error in linear_solve (zsysv)!!' ! stop 'error in linear_solve (zsysv)!!'
!
end if ! end if
deallocate(work,ipiv,rwork,AF) deallocate(work,ipiv,rwork,AF)
end subroutine end subroutine