diff --git a/src/AOtoMO/complex_MOtoAO.f90 b/src/AOtoMO/complex_MOtoAO.f90 index 8a537bd..9f1c874 100644 --- a/src/AOtoMO/complex_MOtoAO.f90 +++ b/src/AOtoMO/complex_MOtoAO.f90 @@ -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)) - cS(:,:) = (0d0,1d0)*S(:,:) + cS(:,:) = (1d0,0d0)*S(:,:) !SC = matmul(S, C) !BSC = matmul(M_MOs, transpose(SC)) !M_AOs = matmul(SC, BSC) diff --git a/src/GW/complex_RGW_self_energy.f90 b/src/GW/complex_RGW_self_energy.f90 index ba7630a..8dc0447 100644 --- a/src/GW/complex_RGW_self_energy.f90 +++ b/src/GW/complex_RGW_self_energy.f90 @@ -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 @@ -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) :: nR integer,intent(in) :: nS - double precision,intent(in) :: Re_e(nOrb) - double precision,intent(in) :: Im_e(nOrb) + complex*16,intent(in) :: e(nOrb) complex*16,intent(in) :: Om(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 double precision :: eps,eta_tilde complex*16 :: num,tmp - double precision, allocatable :: Re_DS(:) - double precision, allocatable :: Im_DS(:) + double precision,allocatable :: Re_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 complex*16,intent(out) :: EcGM - double precision,intent(out) :: Re_Sig(nOrb,nOrb) - double precision,intent(out) :: Im_Sig(nOrb,nOrb) - double precision,intent(out) :: Re_Z(nOrb) - double precision,intent(out) :: Im_Z(nOrb) + complex*16,intent(out) :: Sig(nOrb,nOrb) + complex*16,intent(out) :: Z(nOrb) !----------------! ! 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 Im_Sig(:,:) = 0d0 + Re_DS(:) = 0d0 + Im_DS(:) = 0d0 ! Occupied part of the correlation self-energy !$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 DEFAULT(NONE) !$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 i=nC+1,nO - eps = Re_e(p) - Re_e(i) + real(Om(m)) - eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) + eps = real(e(p)) - real(e(i)) + real(Om(m)) + eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(m)) num = 2d0*rho(p,i,m)*rho(q,i,m) tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& 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 !$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 DEFAULT(NONE) !$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 a=nO+1,nOrb-nR - eps = Re_e(p) - Re_e(a) - real(Om(m)) - eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) + eps = real(e(p)) - real(e(a)) - real(Om(m)) + eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(m)) num = 2d0*rho(p,a,m)*rho(q,a,m) tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& -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 ! !------------------------! - Re_DS(:) = 0d0 - Im_DS(:) = 0d0 - ! Occupied part of the renormalization factor -!$OMP PARALLEL & -!$OMP SHARED(Re_DS,Im_DS,rho,eta,nS,nC,nO,nOrb,nR,Re_e,Im_e,Om) & -!$OMP PRIVATE(m,i,p,eps,num,eta_tilde,tmp) & -!$OMP DEFAULT(NONE) -!$OMP DO + !$OMP PARALLEL & + !$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 DEFAULT(NONE) + !$OMP DO do p=nC+1,nOrb-nR do m=1,nS do i=nC+1,nO - eps = Re_e(p) - Re_e(i) + real(Om(m)) - eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) + eps = real(e(p)) - real(e(i)) + real(Om(m)) + eta_tilde = eta - aimag(e(p)) + aimag(e(i)) - aimag(Om(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,& -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 !$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 DEFAULT(NONE) !$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 a=nO+1,nOrb-nR - eps = Re_e(p) - Re_e(a) - real(Om(m)) - eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) + eps = real(e(p)) - real(e(a)) - real(Om(m)) + eta_tilde = eta + aimag(e(p)) - aimag(e(a)) - aimag(Om(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,& 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 Re_Z(:) = (1d0-Re_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 ! diff --git a/src/GW/complex_qsRGW.f90 b/src/GW/complex_qsRGW.f90 index 2e326b2..494fa7d 100644 --- a/src/GW/complex_qsRGW.f90 +++ b/src/GW/complex_qsRGW.f90 @@ -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) :: X(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) complex*16,intent(inout) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) 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" !call complex_RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) else - call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,real(eGW),aimag(eGW),Om,rho,& - EcGM,real(SigC),aimag(SigC),real(Z),aimag(Z)) + call complex_RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,& + EcGM,SigC,Z) end if - ! Make correlation self-energy Hermitian and transform it back to AO basis 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_MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1)) endif - + ! Compute commutator and convergence criteria 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)) cp(:,:) = Fp(:,:) call complex_diagonalize_matrix(nOrb,cp,eGW) + call complex_orthogonalize_matrix(nBas,cp) c = matmul(X,cp) else Fp = matmul(transpose(c),matmul(F,c)) cp(:,:) = Fp(:,:) call complex_diagonalize_matrix(nOrb,cp,eGW) + call complex_orthogonalize_matrix(nBas,cp) c = matmul(c,cp) endif call complex_complex_AOtoMO(nBas,nOrb,c,SigCp,SigC) - ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 95fe1df..971d9f5 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -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) - if (info /= 0) then - - print *, info - stop 'error in linear_solve (zsysv)!!' - - end if +! if (info /= 0) then +! +! print *, info +! stop 'error in linear_solve (zsysv)!!' +! +! end if deallocate(work,ipiv,rwork,AF) end subroutine