diff --git a/src/AOtoMO/complex_MOtoAO.f90 b/src/AOtoMO/complex_MOtoAO.f90 new file mode 100644 index 0000000..8a537bd --- /dev/null +++ b/src/AOtoMO/complex_MOtoAO.f90 @@ -0,0 +1,39 @@ +subroutine complex_MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs) + + ! Perform MO to AO transformation of a matrix M_AOs for a given metric S + ! and coefficients c + ! + ! M_AOs = S C M_MOs (S C).T + + implicit none + + integer, intent(in) :: nBas, nOrb + double precision, intent(in) :: S(nBas,nBas) + complex*16, intent(in) :: C(nBas,nOrb) + complex*16, intent(in) :: M_MOs(nOrb,nOrb) + complex*16, intent(out) :: M_AOs(nBas,nBas) + + complex*16, allocatable :: SC(:,:),BSC(:,:),cS(:,:) + + + allocate(SC(nBas,nOrb), BSC(nOrb,nBas),cS(nBas,nBas)) + cS(:,:) = (0d0,1d0)*S(:,:) + !SC = matmul(S, C) + !BSC = matmul(M_MOs, transpose(SC)) + !M_AOs = matmul(SC, BSC) + + call zgemm("N", "N", nBas, nOrb, nBas, 1.d0, & + cS(1,1), nBas, C(1,1), nBas, & + 0.d0, SC(1,1), nBas) + + call zgemm("N", "T", nOrb, nBas, nOrb, 1.d0, & + M_MOs(1,1), nOrb, SC(1,1), nBas, & + 0.d0, BSC(1,1), nOrb) + + call zgemm("N", "N", nBas, nBas, nOrb, 1.d0, & + SC(1,1), nBas, BSC(1,1), nOrb, & + 0.d0, M_AOs(1,1), nBas) + + deallocate(SC, BSC) + +end subroutine diff --git a/src/AOtoMO/complex_complex_AOtoMO.f90 b/src/AOtoMO/complex_complex_AOtoMO.f90 new file mode 100644 index 0000000..9acb1b9 --- /dev/null +++ b/src/AOtoMO/complex_complex_AOtoMO.f90 @@ -0,0 +1,24 @@ +subroutine complex_complex_AOtoMO(nBas, nOrb, C, M_AOs, M_MOs) + + ! Perform AO to MO transformation of a matrix M_AOs for given coefficients c + ! M_MOs = C.T M_AOs C + + implicit none + + integer, intent(in) :: nBas, nOrb + complex*16, intent(in) :: C(nBas,nOrb) + complex*16, intent(in) :: M_AOs(nBas,nBas) + + complex*16, intent(out) :: M_MOs(nOrb,nOrb) + + complex*16, allocatable :: AC(:,:) + complex*16, allocatable :: complex_C(:,:) + + allocate(AC(nBas,nOrb)) + + AC = matmul(M_AOs, C) + M_MOs = matmul(transpose(C), AC) + + deallocate(AC) + +end subroutine diff --git a/src/GW/complex_RGW.f90 b/src/GW/complex_RGW.f90 index 40b1225..72c0283 100644 --- a/src/GW/complex_RGW.f90 +++ b/src/GW/complex_RGW.f90 @@ -1,7 +1,7 @@ -subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, & +subroutine complex_RGW(dotest,docG0W0,doevGW,doqsGW,maxSCF,thresh,max_diis,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & - S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_AO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Restricted GW module @@ -12,7 +12,7 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, logical,intent(in) :: dotest - logical,intent(in) :: docG0W0,doevGW + logical,intent(in) :: docG0W0,doevGW,doqsGW integer,intent(in) :: maxSCF integer,intent(in) :: max_diis @@ -57,6 +57,7 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) complex*16,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: CAP_AO(nOrb,nOrb) complex*16,intent(in) :: CAP_MO(nOrb,nOrb) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) complex*16,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) @@ -100,4 +101,24 @@ subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, end if +!------------------------------------------------------------------------ +! Perform qsGW calculation +!------------------------------------------------------------------------ + + if(doqsGW) then + + call wall_time(start_GW) + call complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, & + TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,nNuc,ZNuc,rNuc, & + ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & + dipole_int_AO,dipole_int_MO,PHF,cHF,eHF, & + CAP_AO,CAP_MO) + call wall_time(end_GW) + + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' + write(*,*) + + end if + end subroutine diff --git a/src/GW/complex_RGW_SRG_self_energy_diag.f90 b/src/GW/complex_RGW_SRG_self_energy_diag.f90 new file mode 100644 index 0000000..c9e9b5f --- /dev/null +++ b/src/GW/complex_RGW_SRG_self_energy_diag.f90 @@ -0,0 +1,99 @@ +subroutine complex_RGW_SRG_self_energy_diag(flow,eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z) + +! Compute diagonal of the correlation part of the self-energy and the renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + double precision,intent(in) :: flow + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: Re_e(nBas) + double precision,intent(in) :: Im_e(nBas) + complex*16,intent(in) :: Om(nS) + complex*16,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,p,m + double precision :: eps + complex*16 :: num + double precision :: eta_tilde + double precision,allocatable :: Re_DS(:) + double precision,allocatable :: Im_DS(:) + complex*16 :: tmp + +! Output variables + + double precision,intent(out) :: Re_Sig(nBas) + double precision,intent(out) :: Im_Sig(nBas) + double precision,intent(out) :: Re_Z(nBas) + double precision,intent(out) :: Im_Z(nBas) + complex*16,intent(out) :: EcGM + +! Initialize + allocate(Re_DS(nBas),Im_DS(nBas)) + Re_Sig(:) = 0d0 + Im_Sig(:) = 0d0 + Re_DS(:) = 0d0 + Im_DS(:) = 0d0 + +!----------------! +! GW self-energy ! +!----------------! + +! Occupied part of the correlation self-energy + do p=nC+1,nBas-nR + do i=nC+1,nO + do m=1,nS + eps = Re_e(p) - Re_e(i) + real(Om(m)) + eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) + num = 2d0*rho(p,i,m)**2 + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& + eta_tilde/(eps**2+eta_tilde**2),kind=8) + Re_Sig(p) = Re_Sig(p) + real(tmp) + Im_Sig(p) = Im_Sig(p) + aimag(tmp) + 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) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + + end do + end do + end do + +! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do m=1,nS + + eps = Re_e(p) - Re_e(a) - real(Om(m)) + eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) + num = 2d0*rho(p,a,m)**2 + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& + -eta_tilde/(eps**2 + eta_tilde**2),kind=8) + Re_Sig(p) = Re_Sig(p) + real(tmp) + Im_Sig(p) = Im_Sig(p) + aimag(tmp) + 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) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + end do + end do + end do + + +! 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) +end subroutine diff --git a/src/GW/complex_RGW_self_energy.f90 b/src/GW/complex_RGW_self_energy.f90 new file mode 100644 index 0000000..ba7630a --- /dev/null +++ b/src/GW/complex_RGW_self_energy.f90 @@ -0,0 +1,178 @@ +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) + +! Compute correlation part of the self-energy and the renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + 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) :: Om(nS) + complex*16,intent(in) :: rho(nOrb,nOrb,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q,m + double precision :: eps,eta_tilde + complex*16 :: num,tmp + double precision, allocatable :: Re_DS(:) + double precision, allocatable :: Im_DS(:) + +! 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) + +!----------------! +! GW self-energy ! +!----------------! + allocate(Re_DS(nBas),Im_DS(nBas)) + + Re_Sig(:,:) = 0d0 + Im_Sig(:,:) = 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 PRIVATE(m,i,q,p,eps,num,eta_tilde,tmp) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + 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)) + 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) + Re_Sig(p,q) = Re_Sig(p,q) + real(tmp) + Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp) + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! 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 PRIVATE(m,a,q,p,eps,num,eta_tilde,tmp) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + 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)) + 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) + Re_Sig(p,q) = Re_Sig(p,q) + real(tmp) + Im_Sig(p,q) = Im_Sig(p,q) + aimag(tmp) + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +!------------------------! +! 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 + 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)) + 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) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! 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 PRIVATE(m,a,p,eps,num,eta_tilde,tmp) & + !$OMP DEFAULT(NONE) + !$OMP DO + do p=nC+1,nOrb-nR + 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)) + 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) + Re_DS(p) = Re_DS(p) + real(tmp) + Im_DS(p) = Im_DS(p) + aimag(tmp) + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! 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) + +!!-------------------------------------! +!! Galitskii-Migdal correlation energy ! +!!-------------------------------------! +! +! EcGM = 0d0 +! do m=1,nS +! do a=nO+1,nOrb-nR +! do i=nC+1,nO +! +! eps = e(a) - e(i) + Om(m) +! num = 4d0*rho(a,i,m)*rho(a,i,m) +! EcGM = EcGM - num*eps/(eps**2 + eta**2) +! +! end do +! end do +! end do +! +end subroutine diff --git a/src/GW/complex_evRGW.f90 b/src/GW/complex_evRGW.f90 index cc20933..f7ca17e 100644 --- a/src/GW/complex_evRGW.f90 +++ b/src/GW/complex_evRGW.f90 @@ -14,7 +14,7 @@ subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,d integer,intent(in) :: max_diis double precision,intent(in) :: thresh double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + complex*16,intent(in) :: ERHF logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS diff --git a/src/GW/complex_qsRGW.f90 b/src/GW/complex_qsRGW.f90 new file mode 100644 index 0000000..2e326b2 --- /dev/null +++ b/src/GW/complex_qsRGW.f90 @@ -0,0 +1,388 @@ +subroutine complex_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, & + TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,eta,doSRG,nNuc,ZNuc,rNuc, & + ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, & + ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF, & + CAP_AO,CAP_MO) + +! Perform a quasiparticle self-consistent GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: doppBSE + logical,intent(in) :: singlet + logical,intent(in) :: triplet + double precision,intent(in) :: eta + logical,intent(in) :: doSRG + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + complex*16,intent(in) :: ERHF + complex*16,intent(in) :: eHF(nOrb) + complex*16,intent(in) :: cHF(nBas,nOrb) + complex*16,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + 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) + 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) + complex*16,intent(inout) :: dipole_int_MO(nOrb,nOrb,ncart) + +! Local variables + + integer :: nSCF + integer :: nBas_Sq + integer :: ispin + integer :: ixyz + integer :: n_diis + complex*16 :: ET + complex*16 :: EV + complex*16 :: EJ + complex*16 :: EK + complex*16 :: EqsGW + complex*16 :: EW + complex*16 :: EcRPA + complex*16 :: EcBSE(nspin) + complex*16 :: EcGM + double precision :: Conv + double precision :: rcond + double precision,external :: complex_trace_matrix + complex*16 :: dipole(ncart) + + double precision :: flow + + logical :: dRPA_W = .true. + logical :: print_W = .false. + complex*16,allocatable :: err_diis(:,:) + complex*16,allocatable :: F_diis(:,:) + complex*16,allocatable :: Aph(:,:) + complex*16,allocatable :: Bph(:,:) + complex*16,allocatable :: Om(:) + complex*16,allocatable :: XpY(:,:) + complex*16,allocatable :: XmY(:,:) + complex*16,allocatable :: rho(:,:,:) + complex*16,allocatable :: c(:,:) + complex*16,allocatable :: cp(:,:) + complex*16,allocatable :: eGW(:) + complex*16,allocatable :: P(:,:) + complex*16,allocatable :: F(:,:) + complex*16,allocatable :: Fp(:,:) + complex*16,allocatable :: J(:,:) + complex*16,allocatable :: K(:,:) + complex*16,allocatable :: SigC(:,:) + complex*16,allocatable :: SigCp(:,:) + complex*16,allocatable :: Z(:) + complex*16,allocatable :: err(:,:) + +! Hello world + + write(*,*) + write(*,*)'*******************************' + write(*,*)'* Restricted qsGW Calculation *' + write(*,*)'*******************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs and CAP in MO basis will be overwritten in qsGW !!' + write(*,*) + +! Stuff + + nBas_Sq = nBas*nBas + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamical screening!' + write(*,*) + end if + +! SRG regularization + + flow = 500d0 + + if(doSRG) then + + write(*,*) '*** SRG regularized qsGW scheme ***' + write(*,*) + + end if + +! Memory allocation + + allocate(eGW(nOrb)) + allocate(Z(nOrb)) + + allocate(c(nBas,nOrb)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(SigC(nOrb,nOrb)) + + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(err(nBas,nBas)) + allocate(SigCp(nBas,nBas)) + + allocate(Aph(nS,nS)) + allocate(Bph(nS,nS)) + allocate(Om(nS)) + allocate(XpY(nS,nS)) + allocate(XmY(nS,nS)) + allocate(rho(nOrb,nOrb,nS)) + + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:) = PHF(:,:) + eGW(:) = eHF(:) + c(:,:) = cHF(:,:) + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + rcond = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Build Hartree-exchange matrix + + call complex_Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) + call complex_exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + + ! AO to MO transformation of two-electron integrals + + do ixyz=1,ncart + call complex_AOtoMO(nBas,nOrb,c,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz)) + end do + + call complex_AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO) + + ! Compute linear response + + call complex_phRLR_A(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call complex_phRLR_B(ispin,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + + call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) + + call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) + + if(doSRG) then + 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)) + end if + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + SigC = 0.5d0*(SigC + transpose(SigC)) + + call complex_MOtoAO(nBas,nOrb,S,c,SigC,SigCp) + + ! Solve the quasi-particle equation + + F(:,:) = cmplx(Hc(:,:),CAP_AO(:,:),kind=8) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) + if(nBas .ne. nOrb) then + 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) + + if(nSCF > 1) Conv = maxval(abs(err)) + + ! Kinetic energy + + ET = complex_trace_matrix(nBas,matmul(P,T)) + + ! Potential energy + + EV = complex_trace_matrix(nBas,matmul(P,V)) + + ! Hartree energy + + EJ = 0.5d0*complex_trace_matrix(nBas,matmul(P,J)) + + ! Exchange energy + + EK = 0.25d0*complex_trace_matrix(nBas,matmul(P,K)) + + ! CAP energy + + EW = complex_trace_matrix(nBas,matmul(P,(0d0,1d0)*CAP_AO)) + + ! Total energy + + EqsGW = ET + EV + EJ + EK + EW + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call complex_DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) + + end if + + ! Diagonalize Hamiltonian in AO basis + + if(nBas .eq. nOrb) then + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call complex_diagonalize_matrix(nOrb,cp,eGW) + c = matmul(X,cp) + else + Fp = matmul(transpose(c),matmul(F,c)) + cp(:,:) = Fp(:,:) + call complex_diagonalize_matrix(nOrb,cp,eGW) + 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))) + + ! Print results + + !call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_complex_qsRGW(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z, & + ENuc,ET,EV,EW,EJ,EK,EcGM,EcRPA,EqsGW,dipole) + + end do +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + +! if(nSCF == maxSCF+1) then +! +! write(*,*) +! write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' +! write(*,*)' Convergence failed ' +! write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' +! write(*,*) +! +! deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) +! stop +! +! end if +! +!! Deallocate memory +! +! deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) +! +!! Perform BSE calculation +! +! if(dophBSE) then +! +! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & +! nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@qsGW@RHF total energy = ',ENuc + EqsGW + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +!! Compute the BSE correlation energy via the adiabatic connection +! +! if(doACFDT) then +! +! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@qsGW@RHF total energy = ',ENuc + EqsGW + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +! end if +! +! if(doppBSE) then +! +! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@qsGW@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +!! Testing zone +! +! if(dotest) then +! +! call dump_test_value('R','qsGW correlation energy',EcRPA) +! call dump_test_value('R','qsGW HOMO energy',eGW(nO)) +! call dump_test_value('R','qsGW LUMO energy',eGW(nO+1)) +! +! end if +! +end subroutine diff --git a/src/GW/print_complex_qsRGW.f90 b/src/GW/print_complex_qsRGW.f90 new file mode 100644 index 0000000..2d435fc --- /dev/null +++ b/src/GW/print_complex_qsRGW.f90 @@ -0,0 +1,148 @@ + +! --- + +subroutine print_complex_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, & + Z, ENuc, ET, EV,EW, EJ, EK, EcGM, EcRPA, EqsGW, dipole) + +! Print useful information about qsRGW calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nO + integer,intent(in) :: nSCF + double precision,intent(in) :: ENuc + complex*16,intent(in) :: ET + complex*16,intent(in) :: EV + complex*16,intent(in) :: EW + complex*16,intent(in) :: EJ + complex*16,intent(in) :: EK + complex*16,intent(in) :: EcGM + complex*16,intent(in) :: EcRPA + double precision,intent(in) :: Conv + double precision,intent(in) :: thresh + complex*16,intent(in) :: eHF(nOrb) + complex*16,intent(in) :: eGW(nOrb) + complex*16,intent(in) :: c(nBas,nOrb) + complex*16,intent(in) :: SigC(nOrb,nOrb) + complex*16,intent(in) :: Z(nOrb) + complex*16,intent(in) :: EqsGW + complex*16,intent(in) :: dipole(ncart) + +! Local variables + + logical :: dump_orb = .false. + integer :: p,ixyz,HOMO,LUMO + complex*16 :: Gap + double precision,external :: complex_trace_matrix + +! Output variables + +! HOMO and LUMO + + HOMO = maxloc(real(eGW(1:nO)),1) + LUMO = minloc(real(eGW(nO+1:nBas)),1) + nO + Gap = eGW(LUMO)-eGW(HOMO) + +! Compute energies + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A20,I1,A1,I1,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation' + elseif(nSCF < 100) then + write(*,'(1X,A20,I2,A1,I2,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation' + else + write(*,'(1X,A20,I3,A1,I3,A17)')' Self-consistent qsG',nSCF,'W',nSCF,'@cRHF calculation' + end if + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','Re(e_HF (eV))','|','Re(Sig_GW) (eV)','|','Re(Z)','|','Re(e_GW) (eV)','|' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','Im(e_HF (eV))','|','Im(Sig_GW) (eV)','|','Im(Z)','|','Im(e_GW) (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nOrb + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',p,'|',real(eHF(p))*HaToeV,'|',real(SigC(p,p))*HaToeV,'|',real(Z(p)),'|',real(eGW(p))*HaToeV,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',p,'|',aimag(eHF(p))*HaToeV,'|',aimag(SigC(p,p))*HaToeV,'|',aimag(Z(p)),'|',aimag(eGW(p))*HaToeV,'|' + write(*,*)'-------------------------------------------------------------------------------' + if(p==nO) then + write(*,*)'-------------------------------------------------------------------------------' + end if + + end do + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO real energy = ',real(eGW(HOMO))*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO imag energy = ',aimag(eGW(HOMO))*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF LUMO real energy = ',real(eGW(LUMO))*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF LUMO imag energy = ',aimag(eGW(LUMO))*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO-LUMO gap = ',real(Gap)*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A3)') 'qsGW@RHF HOMO-LUMO gap = ',aimag(Gap)*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF total real energy = ',ENuc + real(EqsGW),' au' + write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF total imag energy = ',aimag(EqsGW),' au' + write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF exchange energy = ',real(EK),' au' + write(*,'(2X,A60,F15.6,A3)') ' qsGW@RHF exchange energy = ',aimag(EK),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Dump results for final iteration + + if(Conv < thresh) then + + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A33)') ' Summary ' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',real(ET) + real(EV) + real(EW),' au' + write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',aimag(ET) + aimag(EV) + aimag(EW),' au' + write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',real(ET),' au' + write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',aimag(ET),' au' + write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',real(EV),' au' + write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',aimag(EV),' au' + write(*,'(A33,1X,F16.10,A3)') ' CAP energy = ',real(EW),' au' + write(*,'(A33,1X,F16.10,A3)') ' CAP energy = ',aimag(EW),' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',real(EJ + EK),' au' + write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',aimag(EJ + EK),' au' + write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',real(EJ),' au' + write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',aimag(EJ),' au' + write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',real(EK),' au' + write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',aimag(EK),' au' + write(*,'(A33,1X,F16.10,A3)') ' Correlation energy = ',real(EcGM),' au' + write(*,'(A33,1X,F16.10,A3)') ' Correlation energy = ',aimag(EcGM),' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',real(EqsGW),' au' + write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',aimag(EqsGW),' au' + write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au' + write(*,'(A33,1X,F16.10,A3)') ' qsRGW energy = ',ENuc + real(EqsGW),' au' + write(*,'(A33,1X,F16.10,A3)') ' qsRGW energy = ',aimag(EqsGW),' au' + write(*,'(A50)') '---------------------------------------' + write(*,*) + + if(dump_orb) then + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' Restricted qsGW orbital coefficients' + write(*,'(A50)') '---------------------------------------' + call complex_matout(nBas, nOrb, c) + write(*,*) + end if + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' Restricted qsGW orbital energies (au) ' + write(*,'(A50)') '---------------------------------------' + call complex_vecout(nOrb, eGW) + write(*,*) + + end if + +end subroutine diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 092dbce..be30a43 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -137,16 +137,28 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,docRHF, allocate(complex_cHF(nBas,nOrb)) allocate(complex_FHF(nBas,nBas)) allocate(complex_dipole_int_MO(nOrb,nOrb,ncart)) + allocate(dipole_int_MO(0,0,0)) allocate(complex_ERI_MO(nOrb,nOrb,nOrb,nOrb)) - if (doCAP) allocate(complex_CAP_MO(nOrb,nOrb)) + allocate(CAP_MO(0,0)) + if (doCAP) then + allocate(complex_CAP_MO(nOrb,nOrb)) + else + allocate(complex_CAP_MO(0,0)) + end if else allocate(PHF(nBas,nBas)) allocate(eHF(nOrb)) allocate(cHF(nBas,nOrb)) allocate(FHF(nBas,nBas)) allocate(dipole_int_MO(nOrb,nOrb,ncart)) + allocate(complex_dipole_int_MO(0,0,0)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) - if (doCAP) allocate(CAP_MO(nOrb,nOrb)) + allocate(complex_CAP_MO(0,0)) + if (doCAP) then + allocate(CAP_MO(nOrb,nOrb)) + else + allocate(CAP_MO(0,0)) + end if end if allocate(ERI_AO(nBas,nBas,nBas,nBas)) @@ -418,10 +430,11 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 if(doGW .and. docRHF) then call wall_time(start_GW) - call complex_RGW(dotest,docG0W0,doevGW,maxSCF_GW,thresh_GW,max_diis_GW, & + call complex_RGW(dotest,docG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW, & doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, & - V,Hc,ERI_AO,complex_ERI_MO,complex_CAP_MO,dipole_int_AO,complex_dipole_int_MO,complex_PHF,complex_cHF,complex_eHF) + V,Hc,ERI_AO,complex_ERI_MO,CAP_AO,complex_CAP_MO,dipole_int_AO,& + complex_dipole_int_MO,complex_PHF,complex_cHF,complex_eHF) call wall_time(end_GW) t_GW = end_GW - start_GW