diff --git a/input/methods.default b/input/methods.default index 5ed2560..134dcbd 100644 --- a/input/methods.default +++ b/input/methods.default @@ -18,8 +18,6 @@ F F F F # G0T0eh evGTeh qsGTeh F F F -# cG0W0 cG0F2 - F F # Parquet F # Rtest Utest Gtest diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 07acb38..13e488b 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -1,4 +1,4 @@ -subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2, & +subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3, & renorm,maxSCF, & thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, & eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & @@ -14,7 +14,6 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2, logical,intent(in) :: dotest logical,intent(in) :: doG0F2 - logical,intent(in) :: docG0F2 logical,intent(in) :: doevGF2 logical,intent(in) :: doqsGF2 logical,intent(in) :: doufG0F02 @@ -171,21 +170,4 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2, end if -!------------------------------------------------------------------------ -! Compute complex G0F2 electronic binding energies -!------------------------------------------------------------------------ - - if(docG0F2) then - - call wall_time(start_GF) - call cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS, & - ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF) - call wall_time(end_GF) - - t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' - write(*,*) - - end if end subroutine diff --git a/src/GF/cRG0F2.f90 b/src/GF/cRG0F2.f90 deleted file mode 100644 index 75a8ba7..0000000 --- a/src/GF/cRG0F2.f90 +++ /dev/null @@ -1,152 +0,0 @@ -subroutine cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & - nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,CAP,dipole_int,eHF) - -! Perform a one-shot second-order Green function calculation - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: dotest - - logical,intent(in) :: dophBSE - logical,intent(in) :: doppBSE - 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) :: regularize - - integer,intent(in) :: nBas - integer,intent(in) :: nOrb - integer,intent(in) :: nO - integer,intent(in) :: nC - 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) :: eHF(nOrb) - double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) - double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) - double precision,intent(in) :: CAP(nOrb,nOrb) - -! Local variables - - integer :: p - double precision :: Ec - double precision :: EcBSE(nspin) - double precision,allocatable :: Re_SigC(:) - double precision,allocatable :: Im_SigC(:) - double precision,allocatable :: Re_Z(:) - double precision,allocatable :: Im_Z(:) - double precision,allocatable :: Re_eGFlin(:) - double precision, allocatable :: Im_eGFlin(:) - double precision,allocatable :: Re_eGF(:) - double precision,allocatable :: Im_eGF(:) - double precision, allocatable :: e_CAP(:) - - ! Hello world - write(*,*) - write(*,*)'*******************************' - write(*,*)'* Restricted G0F2 Calculation *' - write(*,*)'*******************************' - write(*,*) - -! Memory allocation - - allocate(Re_SigC(nOrb),Im_SigC(nOrb), Re_Z(nOrb),Im_Z(nOrb),& - Re_eGFlin(nOrb),Im_eGFlin(nOrb), Re_eGF(nOrb),Im_eGF(nOrb),e_CAP(nOrb)) - do p = 1, nOrb - e_CAP(p) = CAP(p,p) - end do - -! Frequency-dependent second-order contribution - - if(regularize) then - write(*,*) "Regularisation not implemented (yet)" - !call RGF2_reg_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z) - - else - - call cRGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_CAP) - - end if - - Re_eGFlin(:) = eHF(:) + Re_Z(:)*Re_SigC(:) - Im_Z(:)*Im_SigC(:) - Im_eGFlin(:) = e_CAP(:) + Re_Z(:)*Im_SigC(:) + Im_Z(:)*Re_SigC(:) - - if(linearize) then - - write(*,*) '*** Quasiparticle energies obtained by linearization ***' - - Re_eGF(:) = Re_eGFlin(:) - Im_eGF(:) = Im_eGFlin(:) - - else - - write(*,*) ' *** Quasiparticle energies obtained by root search *** ' - write(*,*) - - call cRGF2_QP_graph(eta,nOrb,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,eHF,e_cap,Re_eGF,Im_eGF,Re_Z,Im_Z) - - end if - - ! Print results - -! call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) - call print_cRG0F2(nOrb,nO,eHF,e_CAP,Re_SigC,Im_SigC,Re_eGF,Im_eGF,Re_Z,Im_Z,ENuc,ERHF,Ec) - -! Perform BSE@GF2 calculation -! -! if(dophBSE) then -! -! call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) -! -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (singlet) =',EcBSE(1) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy (triplet) =',EcBSE(2) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 correlation energy =',sum(EcBSE) -! write(*,'(2X,A50,F20.10)') 'Tr@phBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) -! -! end if -! -!! Perform ppBSE@GF2 calculation -! -! if(doppBSE) then -! -! call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) -! -! EcBSE(2) = 3d0*EcBSE(2) -! -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy (triplet) =',EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 correlation energy =',sum(EcBSE),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0F2 total energy =',ENuc + ERHF + sum(EcBSE),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) -! -! end if -! -!! Testing zone -! -! if(dotest) then -! -! call dump_test_value('R','G0F2 correlation energy',Ec) -! call dump_test_value('R','G0F2 HOMO energy',eGF(nO)) -! call dump_test_value('R','G0F2 LUMO energy',eGF(nO+1)) -! -! end if -! -! deallocate(SigC, Z, eGFlin, eGF) -! -end subroutine diff --git a/src/GF/cRGF2_Im_SigC.f90 b/src/GF/cRGF2_Im_SigC.f90 deleted file mode 100644 index 6583314..0000000 --- a/src/GF/cRGF2_Im_SigC.f90 +++ /dev/null @@ -1,59 +0,0 @@ -double precision function cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: i,j,a,b - double precision :: eps - double precision :: eta_tilde - double precision :: num - -! Initialize - - cRGF2_Im_SigC = 0d0 - - ! Occupied part of the correlation self-energy - - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - - eps = Re_w + eHF(a) - eHF(i) - eHF(j) - eta_tilde = eta - Im_w + e_cap(i) -( e_cap(a) - e_cap(j)) - num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - cRGF2_Im_SigC = cRGF2_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2) - - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = Re_w + eHF(i) - eHF(a) - eHF(b) - num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) - cRGF2_Im_SigC = cRGF2_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2) - end do - end do - end do - -end function diff --git a/src/GF/cRGF2_Im_dSigC.f90 b/src/GF/cRGF2_Im_dSigC.f90 deleted file mode 100644 index 9bdf7e7..0000000 --- a/src/GF/cRGF2_Im_dSigC.f90 +++ /dev/null @@ -1,56 +0,0 @@ -double precision function cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: i,j,a,b - double precision :: eps - double precision :: eta_tilde - double precision :: num - - ! Initialize - - cRGF2_Im_dSigC = 0d0 - - ! Occupied part of the correlation self-energy - - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - eps = Re_w + eHF(a) - eHF(i) - eHF(j) - eta_tilde = eta - Im_w + e_cap(i) -( e_cap(a) - e_cap(j)) - num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - cRGF2_Im_dSigC = cRGF2_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - eps = Re_w + eHF(i) - eHF(a) - eHF(b) - num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) - cRGF2_Im_dSigC = cRGF2_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 - end do - end do - end do - -end function diff --git a/src/GF/cRGF2_QP_graph.f90 b/src/GF/cRGF2_QP_graph.f90 deleted file mode 100644 index 2f2a6ff..0000000 --- a/src/GF/cRGF2_QP_graph.f90 +++ /dev/null @@ -1,97 +0,0 @@ -subroutine cRGF2_QP_graph(eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI,Re_eGFlin,Im_eGFlin,Re_eOld,Im_eold,Re_eGF,Im_eGF,Re_Z,Im_Z) - -! Compute the graphical solution of the GF2 QP equation - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: Re_eGFlin(nBas) - double precision,intent(in) :: Im_eGFlin(nBas) - double precision,intent(in) :: Re_eOld(nBas) - double precision,intent(in) :: Im_eOld(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: p - integer :: nIt - integer,parameter :: maxIt = 64 - double precision,parameter :: thresh = 1d-6 - double precision,external :: cRGF2_Re_SigC,cRGF2_Im_SigC,cRGF2_Re_dSigC,cRGF2_Im_dSigC - double precision :: Re_SigC,Im_SigC,Re_dSigC,Im_dSigC - double precision :: Re_f,Im_f,Re_df,Im_df - double precision :: Re_w,Im_w - -! Output variables - - double precision,intent(out) :: Re_eGF(nBas),Im_eGF(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,A15,1X,A15,1X,A10)') 'Orb.','It.','Re(e_GFlin) (eV)','Re(e_GF) (eV)','Re(Z)' - write(*,'(A5,1X,A3,1X,A15,1X,A15,1X,A10)') 'Orb.','It.','Im(e_GFlin) (eV)','Im(e_GF) (eV)','Im(Z)' - write(*,*)'-----------------------------------------------------' - - do p=nC+1,nBas-nR - - Re_w = Re_eGFlin(p) - Im_w = Im_eGFlin(p) - nIt = 0 - Re_f = 1d0 - Im_f = 0d0 - - do while (abs(cmplx(Re_f,Im_f,kind=8)) > thresh .and. nIt < maxIt) - - nIt = nIt + 1 - - - Re_SigC = cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) - Im_SigC = cRGF2_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) - Re_dSigC = cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) - Im_dSigC = cRGF2_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,Re_eOld,Im_eOld,ERI) - - Re_f = Re_w - eHF(p) - Re_SigC - Im_f = Im_w - e_cap(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_eGF(p) = Re_eGFlin(p) - Im_eGF(p) = Im_eGFlin(p) - write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Re_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p),'Cvg Failed!' - write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6,1X,A12)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p),'Cvg Failed!' - - else - - Re_eGF(p) = Re_w - Im_eGF(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_eGFlin(p)*HaToeV,Re_eGF(p)*HaToeV,Re_Z(p) - write(*,'(I5,1X,I3,1X,F15.9,1X,F15.9,1X,F10.6)') p,nIt,Im_eGFlin(p)*HaToeV,Im_eGF(p)*HaToeV,Im_Z(p) - - write(*,*)'-----------------------------------------------------' - end if - - end do - -end subroutine diff --git a/src/GF/cRGF2_Re_SigC.f90 b/src/GF/cRGF2_Re_SigC.f90 deleted file mode 100644 index eafec42..0000000 --- a/src/GF/cRGF2_Re_SigC.f90 +++ /dev/null @@ -1,58 +0,0 @@ -double precision function cRGF2_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: i,j,a,b - double precision :: eps - double precision :: eta_tilde - double precision :: num - -! Initialize - - cRGF2_Re_SigC = 0d0 - - ! Occupied part of the correlation self-energy - - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - - eps = Re_w + eHF(a) - eHF(i) - eHF(j) - eta_tilde = eta - Im_w + e_cap(i) -(e_cap(a) - e_cap(j)) - num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2) - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = Re_w + eHF(i) - eHF(a) - eHF(b) - num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) - cRGF2_Re_SigC = cRGF2_Re_SigC + num*eps/(eps**2 + eta_tilde**2) - end do - end do - end do - -end function diff --git a/src/GF/cRGF2_Re_dSigC.f90 b/src/GF/cRGF2_Re_dSigC.f90 deleted file mode 100644 index 4f75e9f..0000000 --- a/src/GF/cRGF2_Re_dSigC.f90 +++ /dev/null @@ -1,59 +0,0 @@ -double precision function cRGF2_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,eHF,e_cap,ERI) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: i,j,a,b - double precision :: eps - double precision :: eta_tilde - double precision :: num -! Initialize - - cRGF2_Re_dSigC = 0d0 - - ! Occupied part of the correlation self-energy - - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - - eps = Re_w + eHF(a) - eHF(i) - eHF(j) - eta_tilde = eta - Im_w + e_cap(i) - (e_cap(a) - e_cap(j)) - cRGF2_Re_dSigC = cRGF2_Re_dSigC -& - (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = Re_w + eHF(i) - eHF(a) - eHF(b) - eta_tilde = eta + Im_w - e_cap(a) - e_cap(b) + e_cap(i) - cRGF2_Re_dSigC = cRGF2_Re_dSigC -& - (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - - end do - end do - end do - -end function diff --git a/src/GF/cRGF2_self_energy_diag.f90 b/src/GF/cRGF2_self_energy_diag.f90 deleted file mode 100644 index 0c701dc..0000000 --- a/src/GF/cRGF2_self_energy_diag.f90 +++ /dev/null @@ -1,87 +0,0 @@ -subroutine cRGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Re_SigC,Im_SigC,Re_Z,Im_Z,e_cap) - -! Compute diagonal part of the GF2 self-energy and its renormalization factor - - implicit none - include 'parameters.h' - -! Input variables - - double precision,intent(in) :: eta - integer,intent(in) :: nBas - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: i,j,a,b - integer :: p - double precision :: eps - double precision :: eta_tilde - double precision :: num - double precision,allocatable :: Re_DS(:) - double precision,allocatable :: Im_DS(:) - -! Output variables - - double precision,intent(out) :: Re_SigC(nBas) - double precision,intent(out) :: Im_SigC(nBas) - double precision,intent(out) :: Re_Z(nBas) - double precision,intent(out) :: Im_Z(nBas) - -! Initialize - allocate(Re_DS(nBas),Im_DS(nBas)) - Re_SigC(:) = 0d0 - Im_SigC(:) = 0d0 - Re_DS(:) = 0d0 - Im_DS(:) = 0d0 - - -! Compute GF2 self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - - eps = e(p) + e(a) - e(i) - e(j) - eta_tilde = eta - e_cap(p) + e_cap(i) - (e_cap(a) - e_cap(j)) - num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - - Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2) - Im_SigC(p) = Im_SigC(p) + num*eta_tilde/(eps**2 + eta_tilde**2) - Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - Im_DS(p) = Im_DS(p) - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 - - end do - end do - end do - end do - - do p=nC+1,nBas-nR - do i=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR - - eps = e(p) + e(i) - e(a) - e(b) - eta_tilde = eta + e_cap(p) - e_cap(a) - e_cap(b) + e_cap(i) - num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - - Re_SigC(p) = Re_SigC(p) + num*eps/(eps**2 + eta_tilde**2) - Im_SigC(p) = Im_SigC(p) - num*eta_tilde/(eps**2 + eta_tilde**2) - Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 - - end do - end do - end do - end do - 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/RGW.f90 b/src/GW/RGW.f90 index c94d04f..1f7baa6 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -1,4 +1,4 @@ -subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh,max_diis,doACFDT, & +subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,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,eGW) @@ -17,7 +17,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh logical,intent(in) :: doqsGW logical,intent(in) :: doufG0W0 logical,intent(in) :: doufGW - logical,intent(in) :: docG0W0 integer,intent(in) :: maxSCF integer,intent(in) :: max_diis @@ -93,22 +92,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,docG0W0,maxSCF,thresh end if -!------------------------------------------------------------------------ -! Perform cG0W0 calculation -!------------------------------------------------------------------------ - - if(docG0W0) then - call wall_time(start_GW) - call 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 G0W0 = ',t_GW,' seconds' - write(*,*) - - end if - !------------------------------------------------------------------------ ! Perform evGW calculation !------------------------------------------------------------------------ diff --git a/src/GW/cRG0W0.f90 b/src/GW/cRG0W0.f90 deleted file mode 100644 index e3ff028..0000000 --- a/src/GW/cRG0W0.f90 +++ /dev/null @@ -1,242 +0,0 @@ -subroutine 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 G0W0 calculation - - 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 phRLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eHF,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) - - if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) - -!--------------------------! -! Compute spectral weights ! -!--------------------------! - - call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) - -!------------------------! -! Compute GW self-energy ! -!------------------------! - call 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/cRGW_Im_SigC.f90 b/src/GW/cRGW_Im_SigC.f90 deleted file mode 100644 index d038c23..0000000 --- a/src/GW/cRGW_Im_SigC.f90 +++ /dev/null @@ -1,56 +0,0 @@ -double precision function cRGW_Im_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - 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) :: Re_e(nBas) - double precision,intent(in) :: Im_e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: i,a,m - double precision :: num,eps,eta_tilde - -! Initialize - - cRGW_Im_SigC = 0d0 - -! Occupied part of the correlation self-energy - - do i=nC+1,nO - do m=1,nS - eps = Re_w - Re_e(i) + Om(m) - eta_tilde = eta - Im_w + Im_e(i) - num = 2d0*rho(p,i,m)**2 - cRGW_Im_SigC = cRGW_Im_SigC + num*eta_tilde/(eps**2 + eta_tilde**2) - end do - end do - -! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - do m=1,nS - eps = Re_w - Re_e(a) - Om(m) - eta_tilde = eta + Im_w - Im_e(a) - num = 2d0*rho(p,a,m)**2 - cRGW_Im_SigC =cRGW_Im_SigC - num*eta_tilde/(eps**2 + eta_tilde**2) - end do - end do - -end function diff --git a/src/GW/cRGW_Im_dSigC.f90 b/src/GW/cRGW_Im_dSigC.f90 deleted file mode 100644 index 13c8fce..0000000 --- a/src/GW/cRGW_Im_dSigC.f90 +++ /dev/null @@ -1,56 +0,0 @@ -double precision function cRGW_Im_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho) - -! Compute the derivative of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - 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) :: Re_e(nBas) - double precision,intent(in) :: Im_e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: i,a,m - double precision :: num,eps,eta_tilde - -! Initialize - - cRGW_Im_dSigC = 0d0 - -! Occupied part of the correlation self-energy - - do i=nC+1,nO - do m=1,nS - eps = Re_w - Re_e(i) + Om(m) - eta_tilde = eta - Im_w + Im_e(i) - num = 2d0*rho(p,i,m)**2 - cRGW_Im_dSigC = cRGW_Im_dSigC - 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 - end do - end do - -! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - do m=1,nS - eps = Re_w - Re_e(a) - Om(m) - eta_tilde = eta + Im_w - Im_e(a) - num = 2d0*rho(p,a,m)**2 - cRGW_Im_dSigC = cRGW_Im_dSigC + 2d0*num*eps*eta_tilde/(eps**2 + eta_tilde**2)**2 - end do - end do - -end function diff --git a/src/GW/cRGW_QP_graph.f90 b/src/GW/cRGW_QP_graph.f90 deleted file mode 100644 index 8e55b3e..0000000 --- a/src/GW/cRGW_QP_graph.f90 +++ /dev/null @@ -1,106 +0,0 @@ -subroutine cRGW_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,nS,eHF,e_cap,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) :: eHF(nBas) - double precision,intent(in) :: e_cap(nBas) - double precision,intent(in) :: Om(nS) - double precision,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 - eHF(p) - Re_SigC - Im_f = Im_w - e_cap(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/cRGW_Re_SigC.f90 b/src/GW/cRGW_Re_SigC.f90 deleted file mode 100644 index ca85e39..0000000 --- a/src/GW/cRGW_Re_SigC.f90 +++ /dev/null @@ -1,57 +0,0 @@ -double precision function cRGW_Re_SigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - 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) :: Re_e(nBas) - double precision,intent(in) :: Im_e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: i,a,m - double precision :: num,eps - double precision :: eta_tilde - -! Initialize - - cRGW_Re_SigC = 0d0 - -! Occupied part of the correlation self-energy - - do i=nC+1,nO - do m=1,nS - eps = Re_w - Re_e(i) + Om(m) - eta_tilde = eta - Im_w + Im_e(i) - num = 2d0*rho(p,i,m)**2 - cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2) - end do - end do - -! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - do m=1,nS - eps = Re_w - Re_e(a) - Om(m) - eta_tilde = eta + Im_w - Im_e(a) - num = 2d0*rho(p,a,m)**2 - cRGW_Re_SigC = cRGW_Re_SigC + num*eps/(eps**2 + eta_tilde**2) - end do - end do - -end function diff --git a/src/GW/cRGW_Re_dSigC.f90 b/src/GW/cRGW_Re_dSigC.f90 deleted file mode 100644 index 33ba960..0000000 --- a/src/GW/cRGW_Re_dSigC.f90 +++ /dev/null @@ -1,56 +0,0 @@ -double precision function cRGW_Re_dSigC(p,Re_w,Im_w,eta,nBas,nC,nO,nV,nR,nS,Re_e,Im_e,Om,rho) - -! Compute the derivative of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: p - double precision,intent(in) :: Re_w - double precision,intent(in) :: Im_w - double precision,intent(in) :: eta - 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) :: Re_e(nBas) - double precision,intent(in) :: Im_e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: i,a,m - double precision :: num,eps,eta_tilde - -! Initialize - - cRGW_Re_dSigC = 0d0 - -! Occupied part of the correlation self-energy - - do i=nC+1,nO - do m=1,nS - eps = Re_w - Re_e(i) + Om(m) - eta_tilde = eta - Im_w + Im_e(i) - num = 2d0*rho(p,i,m)**2 - cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - end do - end do - -! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - do m=1,nS - eps = Re_w - Re_e(a) - Om(m) - eta_tilde = eta + Im_w - Im_e(a) - num = 2d0*rho(p,a,m)**2 - cRGW_Re_dSigC = cRGW_Re_dSigC - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - end do - end do - -end function diff --git a/src/GW/cRGW_self_energy_diag.f90 b/src/GW/cRGW_self_energy_diag.f90 deleted file mode 100644 index 4f940ff..0000000 --- a/src/GW/cRGW_self_energy_diag.f90 +++ /dev/null @@ -1,102 +0,0 @@ -subroutine cRGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,Re_Sig,Im_Sig,Re_Z,Im_Z,e_cap) - -! 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) :: e(nBas) - double precision,intent(in) :: Om(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: e_cap(nBas) - -! Local variables - - integer :: i,a,p,m - double precision :: num,eps - double precision :: eta_tilde - double precision,allocatable :: Re_DS(:) - double precision,allocatable :: Im_DS(:) - -! 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) - double precision,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 = e(p) - e(i) + Om(m) - eta_tilde = eta - e_cap(p) + e_cap(i) - num = 2d0*rho(p,i,m)**2 - Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2) - Im_Sig(p) = Im_Sig(p) + num*eta_tilde/(eps**2 + eta_tilde**2) - Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - Im_DS(p) = Im_DS(p) - 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 - - 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 = e(p) - e(a) - Om(m) - eta_tilde = eta + e_cap(p) - e_cap(a) - num = 2d0*rho(p,a,m)**2 - Re_Sig(p) = Re_Sig(p) + num*eps/(eps**2 + eta_tilde**2) - Im_Sig(p) = Im_Sig(p) - num*eta_tilde/(eps**2 + eta_tilde**2) - Re_DS(p) = Re_DS(p) - num*(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2 - Im_DS(p) = Im_DS(p) + 2*num*eta_tilde*eps/(eps**2 + eta_tilde**2)**2 - end do - end do - end do - -! Galitskii-Migdal correlation energy -! MAYBE MODIFY THIS TERM - EcGM = 0d0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do m=1,nS - - eps = e(a) - e(i) + Om(m) - num = 4d0*rho(a,i,m)**2 - EcGM = EcGM - num*eps/(eps**2 + eta**2) - - 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/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 5781bcb..1f80ad4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -137,12 +137,14 @@ program QuAcK doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, & doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0,docG0F2, & doParquet, & doRtest,doUtest,doGtest) - doCAP = docG0W0 .or. docG0F2 .or. docRHF ! Add different cases if they need CAP - docG0W0 = docG0W0 .or. (doG0W0 .and. docRHF) - docG0F2 = docG0F2 .or. (doG0F2 .and. docRHF) +! Determine complex function calls + + doCAP = docRHF + docG0W0 = doG0W0 .and. docRHF + docG0F2 = doG0F2 .and. docRHF + !--------------------------! ! Read options for methods ! !--------------------------! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 519e08c..4db910a 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -386,7 +386,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 if(doGF .and. .not. docRHF) then call wall_time(start_GF) - call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,renorm_GF,maxSCF_GF, & + call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,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,complex_ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -424,7 +424,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 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, & + call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,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,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF,eGW) diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 873abb0..889df88 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -10,7 +10,6 @@ subroutine read_methods(working_dir, & doG0W0,doevGW,doqsGW,doufG0W0,doufGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp, & doG0T0eh,doevGTeh,doqsGTeh, & - docG0W0,docG0F2, & doParquet, & doRtest,doUtest,doGtest) @@ -34,7 +33,6 @@ subroutine read_methods(working_dir, & logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(out) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(out) :: doG0T0eh,doevGTeh,doqsGTeh - logical,intent(out) :: docG0W0,docG0F2 logical,intent(out) :: doParquet logical,intent(out) :: doRtest,doUtest,doGtest @@ -207,17 +205,6 @@ subroutine read_methods(working_dir, & if(ans2 == 'T') doevGTeh = .true. if(ans3 == 'T') doqsGTeh = .true. - - ! Read Complex methods - - docG0W0 = .false. - docG0F2 = .false. - - read(1,*) - read(1,*) ans1,ans2 - if(ans1 == 'T') docG0W0 = .true. - if(ans2 == 'T') docG0F2 = .true. - ! Read coupled channels methods doParquet = .false.