diff --git a/src/GW/complex_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90 new file mode 100644 index 0000000..17fcaab --- /dev/null +++ b/src/GW/complex_cRG0W0.f90 @@ -0,0 +1,245 @@ +subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & + linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF) + +! Perform a fully complex G0W0 calculation with CAP + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: doppBSE + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: doSRG + + 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) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: CAP(nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) + double precision,intent(in) :: eHF(nOrb) + +! Local variables + + logical :: print_W = .false. + logical :: plot_self = .false. + logical :: dRPA_W + integer :: isp_W + integer :: p + double precision :: flow + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcGM + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Re_SigC(:) + double precision,allocatable :: Im_SigC(:) + double precision,allocatable :: Re_Z(:) + double precision,allocatable :: Im_Z(:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) + + + double precision,allocatable :: Re_eGWlin(:) + double precision, allocatable :: Im_eGWlin(:) + double precision,allocatable :: Re_eGW(:) + double precision,allocatable :: Im_eGW(:) + double precision, allocatable :: e_cap(:) + +! Hello world + + write(*,*) + write(*,*)'***************************************' + write(*,*)'* Restricted complex G0W0 Calculation *' + write(*,*)'***************************************' + write(*,*) + +! Spin manifold and TDA for dynamical screening + + isp_W = 1 + dRPA_W = .true. + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamical screening!' + write(*,*) + end if + +! SRG regularization + + flow = 500d0 + + if(doSRG) then + ! Not implemented + write(*,*) '*** SRG regularized G0W0 scheme ***' + write(*,*) '!!! No SRG with cRG0W0 !!!' + write(*,*) + + end if + +! Memory allocation + + allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), & + Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),e_cap(nOrb)) + do p = 1, nOrb + e_cap(p) = CAP(p,p) + end do +!-------------------! +! Compute screening ! +!-------------------! + + call complex_phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + if(.not.TDA_W) call complex_phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + !if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) + +!--------------------------! +! Compute spectral weights ! +!--------------------------! + + call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + +!------------------------! +! Compute GW self-energy ! +!------------------------! + +!!!! STOPPED HERE PROCEED WITH IMPLEMENTING COMPLEX_CRGW_SELF_ENERGY_DIAG +!!!! + call complex_cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap) + +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! + + ! Linearized or graphical solution? + Re_eGWlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) + Im_eGWlin(:) = e_cap(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + Re_eGW(:) = Re_eGWlin(:) + Im_eGW(:) = Im_eGWlin(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search *** ' + write(*,*) + + call cRGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,eHF,e_cap,Om,rho,Re_eGWlin,Im_eGWlin,eHF,e_cap,Re_eGW,Im_eGW,Re_Z,Im_Z) + end if + +! Plot self-energy, renormalization factor, and spectral function + + if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) + +! Cumulant expansion + +! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z) + +! Compute the RPA correlation energy + + call phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,Re_eGW,ERI,Aph) + if(.not.TDA_W) call phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + +!--------------! +! Dump results ! +!--------------! + + call print_cRG0W0(nOrb,nO,eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM,CAP) +!---------------------------! +! Perform phBSE 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,dipole_int,eHF,Re_eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem +! +! if(doACFDT) then +! +! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +! end if +! +!!---------------------------! +!! Perform ppBSE calculation ! +!!---------------------------! +! +! if(doppBSE) then +! +! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,Re_eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +!! Testing zone +! +! if(dotest) then +! +! call dump_test_value('R','G0W0 correlation energy',EcRPA) +! call dump_test_value('R','G0W0 HOMO energy',Re_eGW(nO)) +! call dump_test_value('R','G0W0 LUMO energy',Re_eGW(nO+1)) +! +! end if +! +end subroutine diff --git a/src/GW/complex_cRGW_excitation_density.f90 b/src/GW/complex_cRGW_excitation_density.f90 new file mode 100644 index 0000000..e2d5d95 --- /dev/null +++ b/src/GW/complex_cRGW_excitation_density.f90 @@ -0,0 +1,80 @@ +subroutine RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + +! Compute excitation densities + + implicit none + +! Input variables + + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nR + integer,intent(in) :: nS + complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + complex*16,intent(in) :: XpY(nS,nS) + +! Local variables + + integer :: ia,jb,p,q,j,b + complex*16, allocatable :: tmp(:,:,:) + +! Output variables + + complex*16,intent(out) :: rho(nOrb,nOrb,nS) + + if(nOrb .lt. 256) then + + allocate(tmp(nOrb,nOrb,nS)) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(p, q, j, b, jb) & + !$OMP SHARED(nOrb, nC, nO, nR, ERI, tmp) + !$OMP DO COLLAPSE(2) + do p = 1, nOrb + do q = 1, nOrb + jb = 0 + do j = nC+1, nO + do b = nO+1, nOrb-nR + jb = jb + 1 + tmp(p,q,jb) = ERI(p,j,q,b) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call zgemm("N", "T", nOrb*nOrb, nS, nS, 1.d0, & + tmp(1,1,1), nOrb*nOrb, XpY(1,1), nS, & + 0.d0, rho(1,1,1), nOrb*nOrb) + + deallocate(tmp) + + else + + rho(:,:,:) = 0d0 + !$OMP PARALLEL & + !$OMP SHARED(nC,nOrb,nR,nO,nS,rho,ERI,XpY) & + !$OMP PRIVATE(q,p,jb,ia) & + !$OMP DEFAULT(NONE) + !$OMP DO + do q=nC+1,nOrb-nR + do p=nC+1,nOrb-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + do ia=1,nS + rho(p,q,ia) = rho(p,q,ia) + ERI(p,j,q,b)*XpY(ia,jb) + end do + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + endif + +end subroutine diff --git a/src/LR/complex_phRLR.f90 b/src/LR/complex_phRLR.f90 new file mode 100644 index 0000000..2c2e85b --- /dev/null +++ b/src/LR/complex_phRLR.f90 @@ -0,0 +1,65 @@ +subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + integer,intent(in) :: nS + complex*16,intent(in) :: Aph(nS,nS) + complex*16,intent(in) :: Bph(nS,nS) + + ! Local variables + + complex*16 :: complex_trace_matrix + complex*16 :: t1, t2 + complex*16,allocatable :: RPA_matrix(:,:) + complex*16,allocatable :: Z(:,:) + complex*16,allocatable :: tmp(:,:) + complex*16,allocatable :: OmOmminus(:) + +! Output variables + + complex*16,intent(out) :: EcRPA + complex*16,intent(out) :: Om(nS) + complex*16,intent(out) :: XpY(nS,nS) + complex*16,intent(out) :: XmY(nS,nS) + + + +! Tamm-Dancoff approximation + + if(TDA) then + + XpY(:,:) = Aph(:,:) + call complex_diagonalize_matrix(nS,XpY,Om) + XpY(:,:) = transpose(XpY(:,:)) + XmY(:,:) = XpY(:,:) + + else + + allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS)) + + call complex_diagonalize_matrix(2*nS,RPA_matrix,Om) + Om(:) = OmOmminus(1:nS) + call complex_vecout(nS,Om) + call complex_vecout(nS,Om(nS+1:2*nS)) + if (maxval(abs(Om(:)+Om(nS+1:2*nS))) > 1e10) & + call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.') + if(minval(abs(Om(:))) < 0d0) & + call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') + XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,nS+1:2*nS) + XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,nS+1:2*nS) + + deallocate(RPA_matrix) + + end if + + ! Compute the RPA correlation energy + + EcRPA = 0.5d0*(sum(Om) - complex_trace_matrix(nS,Aph)) + +end subroutine diff --git a/src/LR/complex_phRLR_A.f90 b/src/LR/complex_phRLR_A.f90 new file mode 100644 index 0000000..eb3f6ee --- /dev/null +++ b/src/LR/complex_phRLR_A.f90 @@ -0,0 +1,109 @@ +subroutine complex_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + +! Compute resonant block of the ph channel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: lambda + complex*16,intent(in) :: e(nBas) + complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_dRPA + + integer :: i,j,a,b,ia,jb + integer :: nn,jb0 + logical :: i_eq_j + complex*16 :: ct1,ct2 + +! Output variables + + complex*16,intent(out) :: Aph(nS,nS) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build A matrix for single manifold + + if(ispin == 1) then + + nn = nBas - nR - nO + ct1 = 2d0 * lambda + ct2 = - (1d0 - delta_dRPA) * lambda + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i, a, j, b, i_eq_j, ia, jb0, jb) & + !$OMP SHARED (nC, nO, nR, nBas, nn, ct1, ct2, e, ERI, Aph) + !$OMP DO COLLAPSE(2) + do i = nC+1, nO + do a = nO+1, nBas-nR + ia = a - nO + (i - nC - 1) * nn + + do j = nC+1, nO + i_eq_j = i == j + jb0 = (j - nC - 1) * nn - nO + + do b = nO+1, nBas-nR + jb = b + jb0 + + Aph(ia,jb) = ct1 * ERI(b,i,j,a) + ct2 * ERI(b,j,a,i) + if(i_eq_j) then + if(a == b) Aph(ia,jb) = Aph(ia,jb) + e(a) - e(i) + endif + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + end if + +! Build A matrix for triplet manifold + + if(ispin == 2) then + + nn = nBas - nR - nO + ct2 = - (1d0 - delta_dRPA) * lambda + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i, a, j, b, i_eq_j, ia, jb0, jb) & + !$OMP SHARED (nC, nO, nR, nBas, nn, ct2, e, ERI, Aph) + !$OMP DO COLLAPSE(2) + do i = nC+1, nO + do a = nO+1, nBas-nR + ia = a - nO + (i - nC - 1) * nn + + do j = nC+1, nO + i_eq_j = i == j + jb0 = (j - nC - 1) * nn - nO + + do b = nO+1, nBas-nR + jb = b + jb0 + + Aph(ia,jb) = ct2 * ERI(b,j,a,i) + if(i_eq_j) then + if(a == b) Aph(ia,jb) = Aph(ia,jb) + e(a) - e(i) + endif + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + end if + +end subroutine diff --git a/src/LR/complex_phRLR_B.f90 b/src/LR/complex_phRLR_B.f90 new file mode 100644 index 0000000..02751cb --- /dev/null +++ b/src/LR/complex_phRLR_B.f90 @@ -0,0 +1,121 @@ +subroutine complex_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) + +! Compute the coupling block of the ph channel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: lambda + complex*16,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_dRPA + + integer :: i,j,a,b,ia,jb + integer :: nn,jb0 + complex*16 :: ct1,ct2 + +! Output variables + + complex*16,intent(out) :: Bph(nS,nS) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build B matrix for singlet manifold + + if(ispin == 1) then + + nn = nBas - nR - nO + ct1 = 2d0 * lambda + ct2 = - (1d0 - delta_dRPA) * lambda + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i, a, j, b, ia, jb0, jb) & + !$OMP SHARED (nC, nO, nR, nBas, nn, ct1, ct2, ERI, Bph) + !$OMP DO COLLAPSE(2) + do i = nC+1, nO + do a = nO+1, nBas-nR + ia = a - nO + (i - nC - 1) * nn + + do j = nC+1, nO + jb0 = (j - nC - 1) * nn - nO + + do b = nO+1, nBas-nR + jb = b + jb0 + + Bph(ia,jb) = ct1 * ERI(b,i,j,a) + ct2 * ERI(b,j,i,a) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + !ia = 0 + !do i=nC+1,nO + ! do a=nO+1,nBas-nR + ! ia = ia + 1 + ! jb = 0 + ! do j=nC+1,nO + ! do b=nO+1,nBas-nR + ! jb = jb + 1 + ! Bph(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) + ! end do + ! end do + ! end do + !end do + + end if + +! Build B matrix for triplet manifold + + if(ispin == 2) then + + nn = nBas - nR - nO + ct2 = - (1d0 - delta_dRPA) * lambda + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i, a, j, b, ia, jb0, jb) & + !$OMP SHARED (nC, nO, nR, nBas, nn, ct2, ERI, Bph) + !$OMP DO COLLAPSE(2) + do i = nC+1, nO + do a = nO+1, nBas-nR + ia = a - nO + (i - nC - 1) * nn + + do j = nC+1, nO + jb0 = (j - nC - 1) * nn - nO + + do b = nO+1, nBas-nR + jb = b + jb0 + + Bph(ia,jb) = ct2 * ERI(b,j,i,a) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +! ia = 0 +! do i=nC+1,nO +! do a=nO+1,nBas-nR +! ia = ia + 1 +! jb = 0 +! do j=nC+1,nO +! do b=nO+1,nBas-nR +! jb = jb + 1 +! Bph(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) +! end do +! end do +! end do +! end do + + end if + +end subroutine