From 9b663eeca6f1a42334847c74c285601db7206bbd Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Wed, 2 Apr 2025 10:45:46 +0200 Subject: [PATCH] First go for fully complex RG0W0 --- src/GW/complex_RGW.f90 | 85 +++++++++ src/GW/complex_RGW_QP_graph.f90 | 107 +++++++++++ src/GW/complex_RGW_excitation_density.f90 | 46 +++++ src/GW/complex_RGW_self_energy_diag.f90 | 97 ++++++++++ src/GW/complex_cRG0W0.f90 | 212 +++++++++++---------- src/GW/complex_cRGW_excitation_density.f90 | 80 -------- src/GW/print_complex_cRG0W0.f90 | 66 +++++++ src/LR/complex_phRLR.f90 | 26 +-- src/QuAcK/RQuAcK.f90 | 16 +- src/utils/complex_orthogonalization.f90 | 47 +++++ src/utils/complex_sort_eigenvalues_RPA.f90 | 36 ++++ src/utils/wrap_lapack.f90 | 41 ++++ 12 files changed, 660 insertions(+), 199 deletions(-) create mode 100644 src/GW/complex_RGW.f90 create mode 100644 src/GW/complex_RGW_QP_graph.f90 create mode 100644 src/GW/complex_RGW_excitation_density.f90 create mode 100644 src/GW/complex_RGW_self_energy_diag.f90 delete mode 100644 src/GW/complex_cRGW_excitation_density.f90 create mode 100644 src/GW/print_complex_cRG0W0.f90 create mode 100644 src/utils/complex_sort_eigenvalues_RPA.f90 diff --git a/src/GW/complex_RGW.f90 b/src/GW/complex_RGW.f90 new file mode 100644 index 0000000..4e18d8a --- /dev/null +++ b/src/GW/complex_RGW.f90 @@ -0,0 +1,85 @@ +subroutine complex_RGW(dotest,docG0W0,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) + +! Restricted GW module + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + logical,intent(in) :: docG0W0 + + 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 + logical,intent(in) :: linearize + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + complex*16,intent(in) :: ERI_MO(nOrb,nOrb,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) + +! Local variables + + double precision :: start_GW ,end_GW ,t_GW + + logical :: doccG0W0,doccGW + +!------------------------------------------------------------------------ +! Perform cG0W0 calculation +!------------------------------------------------------------------------ + + if(docG0W0) then + call wall_time(start_GW) + call 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_MO,CAP_MO,dipole_int_MO,eHF) + call wall_time(end_GW) + + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for complex G0W0 = ',t_GW,' seconds' + write(*,*) + + end if +end subroutine diff --git a/src/GW/complex_RGW_QP_graph.f90 b/src/GW/complex_RGW_QP_graph.f90 new file mode 100644 index 0000000..3916927 --- /dev/null +++ b/src/GW/complex_RGW_QP_graph.f90 @@ -0,0 +1,107 @@ +subroutine complex_RGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,& + Om,rho,Re_eGWlin,Im_eGWlin, & + Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_Z) + +! Compute the graphical solution of the QP equation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + logical,intent(in) :: doSRG + double precision,intent(in) :: eta + double precision,intent(in) :: flow + double precision,intent(in) :: Re_eHF(nBas) + double precision,intent(in) :: Im_eHF(nBas) + complex*16,intent(in) :: Om(nS) + complex*16,intent(in) :: rho(nBas,nBas,nS) + + double precision,intent(in) :: Re_eGWlin(nBas) + double precision,intent(in) :: Im_eGWlin(nBas) + double precision,external :: cRGW_Re_SigC,cRGW_Re_dSigC + double precision,external :: cRGW_Im_SigC,cRGW_Im_dSigC + double precision,intent(in) :: Re_eOld(nBas) + double precision,intent(in) :: Im_eOld(nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision :: Re_SigC,Re_dSigC + double precision :: Im_SigC,Im_dSigC + double precision :: Re_f,Im_f,Re_df,Im_df + double precision :: Re_w + double precision :: Im_w + +! Output variables + + double precision,intent(out) :: Re_eGW(nBas),Im_eGW(nBas) + double precision,intent(out) :: Re_Z(nBas),Im_Z(nBas) + +! Run Newton's algorithm to find the root + + write(*,*)'-----------------------------------------------------' + write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Re(e_GWlin) (eV)','Re(e_GW (eV))','Re(Z)' + write(*,'(A5,1X,A3,1X,A16,1X,A16,1X,A10)') 'Orb.','It.','Im(e_GWlin) (eV)','Im(e_GW (eV))','Im(Z)' + write(*,*)'-----------------------------------------------------' + + do p=nC+1,nBas-nR + + Re_w = Re_eGWlin(p) + Im_w = Im_eGWlin(p) + nIt = 0 + Re_f = 1d0 + Im_f = 1d0 + + do while (sqrt(Re_f**2+Im_f**2) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + +! Re_SigC = cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho) +! Im_SigC = cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho) +! Re_dSigC = cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho) +! Im_dSigC = cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_eOld,Im_eold,Om,rho) +! + + Re_f = Re_w - Re_eHF(p) - Re_SigC + Im_f = Im_w - Im_eHF(p) - Im_SigC + Re_df = (1d0 - Re_dSigC)/((1d0 - Re_dSigC)**2 + Im_dSigC**2) + Im_df = Im_dSigC/((1d0 - Re_dSigC)**2 + Im_dSigC**2) + Re_w = Re_w - Re_df*Re_f + Im_df*Im_f + Im_w = Im_w - Re_f*Im_df - Re_df*Im_f + + end do + + if(nIt == maxIt) then + + Re_eGW(p) = Re_eGWlin(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p),'Cvg Failed!' + + else + + Re_eGW(p) = Re_w + Im_eGW(p) = Im_w + Re_Z(p) = Re_df + Im_Z(p) = Im_df + + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Re_eGWlin(p)*HaToeV,Re_eGW(p)*HaToeV,Re_Z(p) + write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGWlin(p)*HaToeV,Im_eGW(p)*HaToeV,Im_Z(p) + + end if + + write(*,*)'-----------------------------------------------------' + end do + write(*,*) + +end subroutine diff --git a/src/GW/complex_RGW_excitation_density.f90 b/src/GW/complex_RGW_excitation_density.f90 new file mode 100644 index 0000000..919f3d4 --- /dev/null +++ b/src/GW/complex_RGW_excitation_density.f90 @@ -0,0 +1,46 @@ +subroutine complex_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 + +! Output variables + + complex*16,intent(out) :: rho(nOrb,nOrb,nS) + + 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 +end subroutine diff --git a/src/GW/complex_RGW_self_energy_diag.f90 b/src/GW/complex_RGW_self_energy_diag.f90 new file mode 100644 index 0000000..1b8d61c --- /dev/null +++ b/src/GW/complex_RGW_self_energy_diag.f90 @@ -0,0 +1,97 @@ +subroutine complex_RGW_self_energy_diag(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 + 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)**2,& + -eta_tilde/(eps**2 + eta_tilde**2)**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_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90 index 17fcaab..c5c8e40 100644 --- a/src/GW/complex_cRG0W0.f90 +++ b/src/GW/complex_cRG0W0.f90 @@ -35,11 +35,11 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, 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) + complex*16,intent(in) :: ERHF + complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + complex*16,intent(in) :: CAP(nOrb,nOrb) + complex*16,intent(in) :: dipole_int(nOrb,nOrb,ncart) + complex*16,intent(in) :: eHF(nOrb) ! Local variables @@ -49,26 +49,29 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, 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(:,:) + complex*16 :: EcRPA + complex*16 :: EcBSE(nspin) + complex*16 :: EcGM + complex*16,allocatable :: Aph(:,:) + complex*16,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(:,:,:) + complex*16,allocatable :: Om(:) + double precision,allocatable :: Re_Om(:) + double precision,allocatable :: Im_Om(:) + complex*16,allocatable :: XpY(:,:) + complex*16,allocatable :: XmY(:,:) + complex*16,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(:) + double precision, allocatable :: Re_eHF(:) + double precision, allocatable :: Im_eHF(:) ! Hello world @@ -103,10 +106,9 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, ! 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 + Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),Re_Om(nS),Im_Om(nS)) + Re_eHF(:) = real(eHF(:)) + Im_eHF(:) = aimag(eHF(:)) !-------------------! ! Compute screening ! !-------------------! @@ -115,6 +117,9 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, 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) + Re_Om(:) = real(Om(:)) + Im_Om(:) = aimag(Om(:)) + call complex_vecout(nS,Om) !if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) @@ -123,22 +128,19 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, !--------------------------! call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) - + write(*,*) rho(1,1,1) !------------------------! ! 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) - + call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,rho,EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) + !-----------------------------------! ! 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(:) + Re_eGWlin(:) = Re_eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) + Im_eGWlin(:) = Im_eHF(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) if(linearize) then @@ -152,94 +154,96 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, 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) + write(*,*) "ROOT SEARCH NOT IMPLEMENTED YET" +! call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS, & +! Re_eHF,Im_eHF,Om,rho,Re_eGWlin,Im_eGWlin,Re_eHF,Im_eHF, & +! 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 +! if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) ! -! 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) +!! Cumulant expansion ! -! 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(*,*) +!! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eHF,eGW,Z) ! -! ! Compute the BSE correlation energy via the adiabatic connection fluctuation dissipation theorem +!! Compute the RPA correlation energy ! -! if(doACFDT) then +! 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 RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eHF,Re_eGW,EcBSE) +! call phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! -! 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 +!!--------------! +!! Dump results ! +!!--------------! ! + call print_complex_cRG0W0(nOrb,nO,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM) !!---------------------------! -!! Perform ppBSE calculation ! +!! Perform phBSE 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 -! +!! +!! 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 deleted file mode 100644 index e2d5d95..0000000 --- a/src/GW/complex_cRGW_excitation_density.f90 +++ /dev/null @@ -1,80 +0,0 @@ -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/GW/print_complex_cRG0W0.f90 b/src/GW/print_complex_cRG0W0.f90 new file mode 100644 index 0000000..0f930c2 --- /dev/null +++ b/src/GW/print_complex_cRG0W0.f90 @@ -0,0 +1,66 @@ +subroutine print_complex_cRG0W0(nBas,nO,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM) + +! Print one-electron energies and other stuff for G0W0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas,nO + double precision,intent(in) :: ENuc + complex*16,intent(in) :: ERHF + complex*16,intent(in) :: EcRPA + complex*16,intent(in) :: EcGM + double precision,intent(in) :: Re_eHF(nBas) + double precision,intent(in) :: Im_eHF(nBas) + double precision,intent(in) :: Re_SigC(nBas) + double precision,intent(in) :: Im_SigC(nBas) + double precision,intent(in) :: Re_Z(nBas) + double precision,intent(in) :: Im_Z(nBas) + double precision,intent(in) :: Re_eGW(nBas) + double precision,intent(in) :: Im_eGW(nBas) + + integer :: p,index_homo,index_lumo + double precision :: Re_eHOMO,Re_eLUMO,Im_eHOMO,Im_eLUMO,Re_Gap,Im_Gap + +! HOMO and LUMO + + index_homo = maxloc(Re_eGW(1:nO),1) + Re_eHOMO = Re_eGW(index_homo) + Im_eHOMO = Im_eGW(index_homo) + index_lumo = minloc(Re_eGW(nO+1:nBas),1) + nO + Re_eLUMO = Re_eGW(index_lumo) + Im_eLUMO = Im_eGW(index_lumo) + Re_Gap = Re_eLUMO-Re_eHOMO + Im_Gap = Im_eLUMO-Im_eHOMO + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' cG0W0@RHF calculation ' + 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,nBas + 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,'|',Re_eHF(p)*HaToeV,'|',Re_SigC(p)*HaToeV,'|',Re_Z(p),'|',Re_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,'|',Im_eHF(p)*HaToeV,'|',Im_SigC(p)*HaToeV,'|',Im_Z(p),'|',Im_eGW(p)*HaToeV,'|' + write(*,*)'-------------------------------------------------------------------------------' + end do + write(*,*) + write(*,*)'---------------------------------------------------------------------------------------------------' + write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO energy = ',Re_eHOMO*HaToeV,' + i*',Im_eHOMO*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF LUMO energy = ',Re_eLUMO*HaToeV,' + i*',Im_eLUMO*HaToeV,' eV' + write(*,'(2X,A60,F15.6,A5,F15.6,A3)') 'cG0W0@RHF HOMO-LUMO gap = ',Re_Gap*HaToeV,' + i*',Im_Gap*HaToeV,' eV' + write(*,*)'---------------------------------------------------------------------------------------------------' +! write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF total energy = ',ENuc + ERHF + EcRPA,' au' +! write(*,'(2X,A60,F15.6,A3)') 'phRPA@cG0W0@RHF correlation energy = ',EcRPA,' au' +! write(*,'(2X,A60,F15.6,A3)') ' GM@cG0W0@RHF total energy = ',ENuc + ERHF + EcGM,' au' +! write(*,'(2X,A60,F15.6,A3)') ' GM@cG0W0@RHF correlation energy = ',EcGM,' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +end subroutine diff --git a/src/LR/complex_phRLR.f90 b/src/LR/complex_phRLR.f90 index 2c2e85b..3e726f0 100644 --- a/src/LR/complex_phRLR.f90 +++ b/src/LR/complex_phRLR.f90 @@ -1,4 +1,4 @@ -subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) +subroutine complex_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute linear response @@ -14,7 +14,7 @@ subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Local variables - complex*16 :: complex_trace_matrix + complex*16,external :: complex_trace_matrix complex*16 :: t1, t2 complex*16,allocatable :: RPA_matrix(:,:) complex*16,allocatable :: Z(:,:) @@ -42,20 +42,22 @@ subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) else allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS)) - - call complex_diagonalize_matrix(2*nS,RPA_matrix,Om) + RPA_matrix(1:nS,1:nS) = Aph(:,:) + RPA_matrix(1:nS,nS+1:2*nS) = Bph(:,:) + RPA_matrix(nS+1:2*nS,1:nS) = -Bph(:,:) + RPA_matrix(nS+1:2*nS,nS+1:2*nS) = -Aph(:,:) + call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus) + call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix) + call complex_normalize_RPA(nS,RPA_matrix) 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) & + if(maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS))) > 1e-12) & 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) - + XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(1:nS,nS+1:2*nS) + XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(1:nS,nS+1:2*nS) + call complex_matout(nS,nS,XpY) + deallocate(RPA_matrix,OmOmminus) end if ! Compute the RPA correlation energy diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 393f9c5..85a8333 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -363,7 +363,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 call wall_time(start_GF) call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, & thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, & - eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & + eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,complex_ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GF) @@ -396,7 +396,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 !-----------! doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. docG0W0 - if(doGW) then + if(doGW .and. .not. docRHF) then call wall_time(start_GW) call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW, & @@ -416,8 +416,18 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 ! complex GW module ! !-------------------! -! IMPLEMENT LATER TO TREAT EVERYTHING COMPLEX, i.e. start from complex HF orbitals + if(doGW .and. docRHF) then + call wall_time(start_GW) + call complex_RGW(dotest,docG0W0,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) + call wall_time(end_GW) + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds' + write(*,*) + end if !-----------------! ! T-matrix module ! !-----------------! diff --git a/src/utils/complex_orthogonalization.f90 b/src/utils/complex_orthogonalization.f90 index 7112063..3154f77 100644 --- a/src/utils/complex_orthogonalization.f90 +++ b/src/utils/complex_orthogonalization.f90 @@ -19,6 +19,53 @@ subroutine complex_orthogonalize_matrix(N, vectors) vectors = matmul(vectors,transpose(Linv)) deallocate(L,Linv) end subroutine +subroutine complex_orthonormalize(N,vectors,A) + + ! Orthonormalize vectors matrix, such that vectors^T A vectors = Identity + ! A and vectors are assumed quadratic NxN matrices + + ! Input variables + implicit none + integer, intent(in) :: N + complex*16, intent(inout) :: vectors(N, N) + complex*16, intent(inout) :: A(N, N) + + ! Local variables + integer :: i, j + complex*16,allocatable :: L(:,:),Linv(:,:) + complex*16 :: proj + complex*16 :: norm + + ! Copy the input matrix to a temporary matrix + allocate(L(N,N),Linv(N,N)) + L = matmul(matmul(transpose(vectors),A),vectors) + call complex_cholesky_decomp(N,L) + call complex_inverse_matrix(N,L,Linv) + vectors = matmul(vectors,transpose(Linv)) + deallocate(L,Linv) +end subroutine +subroutine complex_normalize_RPA(nS,XYYX) + + ! Orthonormalize vectors matrix, such that RPA^T (1 0; 0 -1) RPA = Identity + + ! Input variables + implicit none + integer, intent(in) :: nS + complex*16, intent(inout) :: XYYX(2*nS, 2*nS) + + ! Local variables + integer :: i + complex*16,allocatable :: A(:,:) + + allocate(A(2*nS,2*nS)) + A(:,:) = (0d0,0d0) + do i=1,nS + A(i,i) = 1 + A(i+nS,i+nS) = -1 + end do + call complex_orthonormalize(2*nS,XYYX,A) + deallocate(A) +end subroutine subroutine complex_gram_schmidt(N, vectors) ! Input variables diff --git a/src/utils/complex_sort_eigenvalues_RPA.f90 b/src/utils/complex_sort_eigenvalues_RPA.f90 new file mode 100644 index 0000000..09320c5 --- /dev/null +++ b/src/utils/complex_sort_eigenvalues_RPA.f90 @@ -0,0 +1,36 @@ +subroutine complex_sort_eigenvalues_RPA(n, evals, evecs) + +! Sorts the eigenvecs and eigenvals like in the rpa problem, i.e. first the positive ones in ascending order than the negative ones +! in descending order + + implicit none + integer, intent(in) :: n + complex*16, intent(inout) :: evals(n) + complex*16, intent(inout) :: evecs(n, n) + integer :: i, j, k + complex*16 :: temp_val + complex*16 :: temp_vec(n) + + ! bubble sort (or any sorting algorithm can be applied) + do i = 1, n-1 + do j = i+1, n + ! sorting condition: first ascending for positives, then by absolute value for negatives + if ((real(evals(i)) > 0 .and. real(evals(j)) > 0 .and. real(evals(i)) > real(evals(j))) .or. & + (real(evals(i)) < 0 .and. real(evals(j)) < 0 .and. abs(real(evals(i))) > abs(real(evals(j)))) .or. & + (real(evals(i)) < 0 .and. real(evals(j)) > 0)) then + + ! swap eigenvalues + temp_val = evals(i) + evals(i) = evals(j) + evals(j) = temp_val + + ! swap corresponding eigenvectors + temp_vec = evecs(:, i) + evecs(:, i) = evecs(:, j) + evecs(:, j) = temp_vec + end if + end do + end do + +end subroutine + diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 2e0fba6..400841f 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -170,6 +170,47 @@ subroutine complex_diagonalize_matrix(N,A,e) end subroutine +subroutine complex_diagonalize_matrix_without_sort(N,A,e) + +! Diagonalize a general complex matrix + + implicit none + +! Input variables + + integer,intent(in) :: N + complex*16,intent(inout) :: A(N,N) + complex*16,intent(out) :: e(N) + +! Local variables + + integer :: lwork,info + double precision,allocatable :: rwork(:) + complex*16,allocatable :: work(:) + complex*16,allocatable :: VL(:,:) + complex*16,allocatable :: VR(:,:) + +! Memory allocation + allocate(work(1),rwork(2*N),VL(1,1),VR(N,N)) + lwork = -1 + call zgeev('N','V',N,A,N,e,VL,1,VR,N,work,lwork,rwork,info) + lwork = max(1,int(real(work(1)))) + + deallocate(work) + allocate(work(lwork)) + + call zgeev('N','V',N,A,N,e,VL,N,VR,N,work,lwork,rwork,info) + call complex_sort_eigenvalues(N,e,VR) + + + deallocate(work) + A = VR + + if(info /= 0) then + print*,'Problem in diagonalize_matrix (zgeev)!!' + end if + +end subroutine subroutine svd(N,A,U,D,Vt)