From 94f5fd416f04a3cdf94722cdb266f5a64516da36 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 16 Sep 2024 22:30:23 +0200 Subject: [PATCH] CC-based GW --- src/GW/CCG0W0.f90 | 353 ----------------------------- src/GW/RGW.f90 | 18 +- src/GW/ccRG0W0.f90 | 302 ++++++++++++++++++++++++ src/GW/ccRG0W0_mat.f90 | 278 +++++++++++++++++++++++ src/GW/{CCGW.f90 => ccRGW.f90} | 2 +- src/GW/ccRGW_mat.f90 | 295 ++++++++++++++++++++++++ src/GW/{ufG0W0.f90 => ufRG0W0.f90} | 10 +- 7 files changed, 891 insertions(+), 367 deletions(-) delete mode 100644 src/GW/CCG0W0.f90 create mode 100644 src/GW/ccRG0W0.f90 create mode 100644 src/GW/ccRG0W0_mat.f90 rename src/GW/{CCGW.f90 => ccRGW.f90} (99%) create mode 100644 src/GW/ccRGW_mat.f90 rename src/GW/{ufG0W0.f90 => ufRG0W0.f90} (98%) diff --git a/src/GW/CCG0W0.f90 b/src/GW/CCG0W0.f90 deleted file mode 100644 index a1d1a4d..0000000 --- a/src/GW/CCG0W0.f90 +++ /dev/null @@ -1,353 +0,0 @@ -subroutine CCG0W0(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) - -! CC-based GW module - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: maxSCF - double precision,intent(in) :: thresh - - integer,intent(in) :: nBas - integer,intent(in) :: nOrb - integer,intent(in) :: nC - integer,intent(in) :: nO - integer,intent(in) :: nV - integer,intent(in) :: nR - 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) - -! Local variables - - integer :: p,q - integer :: i,j,k,l - integer :: a,b,c,d - - integer :: nSCF - double precision :: Conv - - double precision,allocatable :: OVVO(:,:,:,:) - double precision,allocatable :: VOOV(:,:,:,:) - - double precision,allocatable :: delta_2h1p(:,:,:,:) - double precision,allocatable :: delta_2p1h(:,:,:,:) - - double precision,allocatable :: V_2h1p(:,:,:,:) - double precision,allocatable :: V_2p1h(:,:,:,:) - - double precision,allocatable :: r_2h1p(:,:,:,:) - double precision,allocatable :: r_2p1h(:,:,:,:) - - double precision,allocatable :: t_2h1p(:,:,:,:) - double precision,allocatable :: t_2p1h(:,:,:,:) - - double precision,allocatable :: x_2h1p(:,:) - double precision,allocatable :: x_2p1h(:,:) - - double precision,allocatable :: eGW(:) - double precision,allocatable :: SigGW(:,:) - double precision,allocatable :: cGW(:,:) - double precision,allocatable :: Z(:) - - integer,allocatable :: order(:) - -! Hello world - - write(*,*) - write(*,*)'*****************************' - write(*,*)'* CC-based G0W0 Calculation *' - write(*,*)'*****************************' - write(*,*) - -! Create integral batches - - allocate(OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) - - OVVO(:,:,:,:) = ERI( 1:nO ,nO+1:nOrb,nO+1:nOrb, 1:nO ) - VOOV(:,:,:,:) = ERI(nO+1:nOrb , 1:nO , 1:nO ,nO+1:nOrb) - -! Form energy denominator and guess amplitudes - - allocate(delta_2h1p(nO,nO,nV,nOrb),delta_2p1h(nO,nV,nV,nOrb)) - allocate(V_2h1p(nOrb,nO,nO,nV),V_2p1h(nOrb,nO,nV,nV)) - allocate(t_2h1p(nO,nO,nV,nOrb),t_2p1h(nO,nV,nV,nOrb)) - allocate(x_2h1p(nOrb,nOrb),x_2p1h(nOrb,nOrb)) - - do k=nC+1,nO - do l=nC+1,nO - do c=1,nV-nR - do p=nC+1,nOrb-nR - - V_2h1p(p,k,l,c) = sqrt(2d0)*ERI(p,nO+c,k,l) - - end do - end do - end do - end do - - do k=nC+1,nO - do c=1,nV-nR - do d=1,nV-nR - do p=nC+1,nOrb-nR - - V_2p1h(p,k,c,d) = sqrt(2d0)*ERI(p,k,nO+d,nO+c) - - end do - end do - end do - end do - -! Initialization - - allocate(r_2h1p(nO,nO,nV,nOrb),r_2p1h(nO,nV,nV,nOrb)) - allocate(eGW(nOrb),SigGW(nOrb,nOrb),cGW(nOrb,nOrb),Z(nOrb)) - allocate(order(nOrb)) - - Conv = 1d0 - nSCF = 0 - eGW(:) = eHF(:) - - t_2h1p(:,:,:,:) = 0d0 - t_2p1h(:,:,:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - write(*,*) - write(*,*)'----------------------------------------------' - write(*,*)'| CC-based G0W0 calculation |' - write(*,*)'----------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','HOMO','|','LUMO','|','Conv','|' - write(*,*)'----------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - - ! Increment - - nSCF = nSCF + 1 - - ! Compute energy differences - - do i=nC+1,nO - do j=nC+1,nO - do a=1,nV-nR - do p=nC+1,nOrb-nR - - delta_2h1p(i,j,a,p) = eGW(i) + eGW(j) - eGW(nO+a) - eHF(p) - - end do - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - do p=nC+1,nOrb-nR - - delta_2p1h(i,a,b,p) = eGW(nO+a) + eGW(nO+b) - eGW(i) - eHF(p) - - end do - end do - end do - end do - - ! Compute intermediates - - x_2h1p(:,:) = 0d0 - - do p=nC+1,nOrb-nR - do q=nC+1,nOrb-nR - - do k=nC+1,nO - do l=nC+1,nO - do c=1,nV-nR - - x_2h1p(p,q) = x_2h1p(p,q) + V_2h1p(q,k,l,c)*t_2h1p(k,l,c,p) - - end do - end do - end do - - end do - end do - - x_2p1h(:,:) = 0d0 - - do p=nC+1,nOrb-nR - do q=nC+1,nOrb-nR - - do k=nC+1,nO - do c=1,nV-nR - do d=1,nV-nR - - x_2p1h(p,q) = x_2p1h(p,q) + V_2p1h(q,k,c,d)*t_2p1h(k,c,d,p) - - end do - end do - end do - - end do - end do - - ! Compute residual for 2h1p sector - - do i=nC+1,nO - do j=nC+1,nO - do a=1,nV-nR - - do p=nC+1,nOrb-nR - - r_2h1p(i,j,a,p) = V_2h1p(p,i,j,a) + delta_2h1p(i,j,a,p)*t_2h1p(i,j,a,p) - - do k=nC+1,nO - do c=1,nV-nR - - r_2h1p(i,j,a,p) = r_2h1p(i,j,a,p) - 2d0*OVVO(j,c,a,k)*t_2h1p(i,k,c,p) - - end do - end do - - do q=nC+1,nOrb-nR - - r_2h1p(i,j,a,p) = r_2h1p(i,j,a,p) - t_2h1p(i,j,a,q)*x_2h1p(p,q) - t_2h1p(i,j,a,q)*x_2p1h(p,q) - - end do - - end do - - end do - end do - end do - - ! Compute residual for 2p1h sector - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - do p=nC+1,nOrb-nR - - r_2p1h(i,a,b,p) = V_2p1h(p,i,a,b) + delta_2p1h(i,a,b,p)*t_2p1h(i,a,b,p) - - do k=nC+1,nO - do c=1,nV-nR - - r_2p1h(i,a,b,p) = r_2p1h(i,a,b,p) + 2d0*VOOV(a,k,i,c)*t_2p1h(k,c,b,p) - - end do - end do - - do q=nC+1,nOrb-nR - - r_2p1h(i,a,b,p) = r_2p1h(i,a,b,p) - t_2p1h(i,a,b,q)*x_2h1p(p,q) - t_2p1h(i,a,b,q)*x_2p1h(p,q) - - end do - - end do - - end do - end do - end do - - ! Check convergence - - Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h))) - - ! Update amplitudes - - t_2h1p(:,:,:,:) = t_2h1p(:,:,:,:) - r_2h1p(:,:,:,:)/delta_2h1p(:,:,:,:) - t_2p1h(:,:,:,:) = t_2p1h(:,:,:,:) - r_2p1h(:,:,:,:)/delta_2p1h(:,:,:,:) - - ! Compute self-energy - - SigGW(:,:) = 0d0 - - do p=nC+1,nOrb-nR - - SigGW(p,p) = SigGW(p,p) + eHF(p) - - do q=nC+1,nOrb-nR - - do i=nC+1,nO - do j=nC+1,nO - do a=1,nV-nR - - SigGW(p,q) = SigGW(p,q) + V_2h1p(p,i,j,a)*t_2h1p(i,j,a,q) - - end do - end do - end do - - do i=nC+1,nO - do a=1,nV-nR - do b=1,nV-nR - - SigGW(p,q) = SigGW(p,q) + V_2p1h(p,i,a,b)*t_2p1h(i,a,b,q) - - end do - end do - end do - - end do - end do - - ! Diagonalize non-Hermitian matrix - - call diagonalize_general_matrix(nOrb,SigGW,eGW,cGW) - - do p=1,nOrb - order(p) = p - end do - - call quick_sort(eGW,order,nOrb) - call set_order(cGW,order,nOrb,nOrb) - - ! Renormalization factor - - Z(:) = 1d0 - - ! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',eGW(nO)*HaToeV,'|',eGW(nO+1)*HaToeV,'|',Conv,'|' - - end do - write(*,*)'----------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - end if - - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' CCGW 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)') & - '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' - write(*,*)'-------------------------------------------------------------------------------' - - do p=1,nOrb - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' - end do - write(*,*)'-------------------------------------------------------------------------------' - -end subroutine diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index fc8eacb..852c5d4 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -68,7 +68,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii double precision :: start_GW ,end_GW ,t_GW - logical :: doCCGW + logical :: doccG0W0,doccGW !------------------------------------------------------------------------ ! Perform G0W0 calculation @@ -131,7 +131,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii call wall_time(start_GW) ! TODO - call ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -161,12 +161,13 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii ! Perform CC-based G0W0 calculation !------------------------------------------------------------------------ - doCCGW = .false. + doccG0W0 = .false. - if(doCCGW) then + if(doccG0W0) then call wall_time(start_GW) - call CCG0W0(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) +! call ccRG0W0(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call ccRG0W0_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -180,12 +181,13 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii ! Perform CC-based GW calculation !------------------------------------------------------------------------ - doCCGW = .false. + doccGW = .true. - if(doCCGW) then + if(doccGW) then call wall_time(start_GW) - call CCGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) +! call ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call ccRGW_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/ccRG0W0.f90 b/src/GW/ccRG0W0.f90 new file mode 100644 index 0000000..6370d5e --- /dev/null +++ b/src/GW/ccRG0W0.f90 @@ -0,0 +1,302 @@ +subroutine ccRG0W0(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! CC-based GW module + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + 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) + +! Local variables + + integer :: p,q + integer :: i,j,k,l + integer :: a,b,c,d + + integer :: nSCF + double precision :: Conv + + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VOOV(:,:,:,:) + + double precision,allocatable :: delta_2h1p(:,:,:) + double precision,allocatable :: delta_2p1h(:,:,:) + + double precision,allocatable :: V_2h1p(:,:,:) + double precision,allocatable :: V_2p1h(:,:,:) + + double precision,allocatable :: r_2h1p(:,:,:) + double precision,allocatable :: r_2p1h(:,:,:) + + double precision,allocatable :: t_2h1p(:,:,:) + double precision,allocatable :: t_2p1h(:,:,:) + + double precision :: x_2h1p + double precision :: x_2p1h + + double precision,allocatable :: eGW(:) + double precision,allocatable :: Z(:) + + +! Hello world + + write(*,*) + write(*,*)'*****************************' + write(*,*)'* CC-based G0W0 Calculation *' + write(*,*)'*****************************' + write(*,*) + +! Create integral batches + + allocate(OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) + + OVVO(:,:,:,:) = ERI( 1:nO ,nO+1:nOrb,nO+1:nOrb, 1:nO ) + VOOV(:,:,:,:) = ERI(nO+1:nOrb , 1:nO , 1:nO ,nO+1:nOrb) + +! Form energy denominator and guess amplitudes + + allocate(delta_2h1p(nO,nO,nV),delta_2p1h(nO,nV,nV)) + allocate(V_2h1p(nO,nO,nV),V_2p1h(nO,nV,nV)) + allocate(t_2h1p(nO,nO,nV),t_2p1h(nO,nV,nV)) + allocate(r_2h1p(nO,nO,nV),r_2p1h(nO,nV,nV)) + allocate(eGW(nOrb),Z(nOrb)) + +! Initialization + + eGW(:) = eHF(:) + +!-------------------------! +! Main loop over orbitals ! +!-------------------------! + + do p=nO,nO + + ! Initialization + + Conv = 1d0 + nSCF = 0 + + t_2h1p(:,:,:) = 0d0 + t_2p1h(:,:,:) = 0d0 + + do k=nC+1,nO + do l=nC+1,nO + do c=1,nV-nR + + V_2h1p(k,l,c) = sqrt(2d0)*ERI(p,nO+c,k,l) + + end do + end do + end do + + do k=nC+1,nO + do c=1,nV-nR + do d=1,nV-nR + + V_2p1h(k,c,d) = sqrt(2d0)*ERI(p,k,nO+d,nO+c) + + end do + end do + end do + + !----------------------! + ! Loop over amplitudes ! + !----------------------! + + write(*,*) + write(*,*)'----------------------------------------------' + write(*,*)'| CC-based G0W0 calculation |' + write(*,*)'----------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','HF','|','G0W0','|','Conv','|' + write(*,*)'----------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Compute energy differences + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + delta_2h1p(i,j,a) = eHF(i) + eHF(j) - eHF(nO+a) - eGW(p) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + delta_2p1h(i,a,b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - eGW(p) + + end do + end do + end do + + ! Compute intermediates + + x_2h1p = 0d0 + + do k=nC+1,nO + do l=nC+1,nO + do c=1,nV-nR + + x_2h1p = x_2h1p + V_2h1p(k,l,c)*t_2h1p(k,l,c) + + end do + end do + end do + + x_2p1h = 0d0 + + do k=nC+1,nO + do c=1,nV-nR + do d=1,nV-nR + + x_2p1h = x_2p1h + V_2p1h(k,c,d)*t_2p1h(k,c,d) + + end do + end do + end do + + ! Compute residual for 2h1p sector + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + r_2h1p(i,j,a) = V_2h1p(i,j,a) + delta_2h1p(i,j,a)*t_2h1p(i,j,a) + + do k=nC+1,nO + do c=1,nV-nR + + r_2h1p(i,j,a) = r_2h1p(i,j,a) - 2d0*OVVO(j,c,a,k)*t_2h1p(i,k,c) + + end do + end do + + r_2h1p(i,j,a) = r_2h1p(i,j,a) - t_2h1p(i,j,a)*x_2h1p - t_2h1p(i,j,a)*x_2p1h + + end do + end do + end do + + ! Compute residual for 2p1h sector + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + r_2p1h(i,a,b) = V_2p1h(i,a,b) + delta_2p1h(i,a,b)*t_2p1h(i,a,b) + + do k=nC+1,nO + do c=1,nV-nR + + r_2p1h(i,a,b) = r_2p1h(i,a,b) + 2d0*VOOV(a,k,i,c)*t_2p1h(k,c,b) + + end do + end do + + r_2p1h(i,a,b) = r_2p1h(i,a,b) - t_2p1h(i,a,b)*x_2h1p - t_2p1h(i,a,b)*x_2p1h + + end do + end do + end do + + ! Check convergence + + Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h))) + + ! Update amplitudes + + t_2h1p(:,:,:) = t_2h1p(:,:,:) - r_2h1p(:,:,:)/delta_2h1p(:,:,:) + t_2p1h(:,:,:) = t_2p1h(:,:,:) - r_2p1h(:,:,:)/delta_2p1h(:,:,:) + + ! Compute self-energy + + eGW(p) = eHF(p) + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + + eGW(p) = eGW(p) + V_2h1p(i,j,a)*t_2h1p(i,j,a) + + end do + end do + end do + + do i=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + eGW(p) = eGW(p) + V_2p1h(i,a,b)*t_2p1h(i,a,b) + + end do + end do + end do + + ! Renormalization factor + + Z(:) = 1d0 + + ! Dump results + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|' + + end do + + write(*,*)'----------------------------------------------' + !------------------------------------------------------------------------ + ! End of SCF loop + !------------------------------------------------------------------------ + + ! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' CC-G0W0 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)') & + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + 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,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' + write(*,*)'-------------------------------------------------------------------------------' + + end do + +end subroutine diff --git a/src/GW/ccRG0W0_mat.f90 b/src/GW/ccRG0W0_mat.f90 new file mode 100644 index 0000000..d75b68d --- /dev/null +++ b/src/GW/ccRG0W0_mat.f90 @@ -0,0 +1,278 @@ +subroutine ccRG0W0_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! CC-based GW module + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + 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) + +! Local variables + + integer :: p,q + integer :: i,j,k,l + integer :: a,b,c,d + + integer :: nSCF + double precision :: Conv + + double precision,allocatable :: delta_2h1p(:) + double precision,allocatable :: delta_2p1h(:) + + double precision,allocatable :: C_2h1p(:,:) + double precision,allocatable :: C_2p1h(:,:) + + double precision,allocatable :: V_2h1p(:) + double precision,allocatable :: V_2p1h(:) + + double precision,allocatable :: r_2h1p(:) + double precision,allocatable :: r_2p1h(:) + + double precision,allocatable :: t_2h1p(:) + double precision,allocatable :: t_2p1h(:) + + double precision,allocatable :: eGW(:) + double precision,allocatable :: Z(:) + + double precision,external :: Kronecker_delta + + integer :: n1h,n1p,n1h1p,n2h1p,n2p1h + integer :: ija,klc,iab,kcd + + +! Hello world + + write(*,*) + write(*,*)'*****************************' + write(*,*)'* CC-based G0W0 Calculation *' + write(*,*)'*****************************' + write(*,*) + +! Form energy denominator and guess amplitudes + + n1h = nO + n1p = nV + n1h1p = n1h*n1p + n2h1p = n1h*n1h1p + n2p1h = n1p*n1h1p + + allocate(delta_2h1p(n2h1p),delta_2p1h(n2p1h)) + allocate(C_2h1p(n2h1p,n2h1p),C_2p1h(n2p1h,n2p1h)) + allocate(V_2h1p(n2h1p),V_2p1h(n2p1h)) + allocate(t_2h1p(n2h1p),t_2p1h(n2p1h)) + allocate(r_2h1p(n2h1p),r_2p1h(n2p1h)) + allocate(eGW(nOrb),Z(nOrb)) + + ! Compute C2h1p and C2p1h + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + ija = ija + 1 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nOrb-nR + klc = klc + 1 + + C_2h1p(ija,klc) = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) & + - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k) + + end do + end do + end do + + end do + end do + end do + + iab = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + iab = iab + 1 + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR + kcd = kcd + 1 + + C_2p1h(iab,kcd) = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d) + + end do + end do + end do + + end do + end do + end do + +!-------------------------! +! Main loop over orbitals ! +!-------------------------! + + eGW(:) = eHF(:) + + do p=nO,nO + + ! Initialization + + Conv = 1d0 + nSCF = 0 + + t_2h1p(:) = 0d0 + t_2p1h(:) = 0d0 + + ! Compute energy differences + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + ija = ija + 1 + + delta_2h1p(ija) = eHF(i) + eHF(j) - eHF(a) - eHF(p) + + end do + end do + end do + + iab = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + iab = iab + 1 + + delta_2p1h(iab) = eHF(a) + eHF(b) - eHF(i) - eHF(p) + + end do + end do + end do + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nOrb-nR + klc = klc + 1 + + V_2h1p(klc) = sqrt(2d0)*ERI(p,c,k,l) + + end do + end do + end do + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR + kcd = kcd + 1 + + V_2p1h(kcd) = sqrt(2d0)*ERI(p,k,d,c) + + end do + end do + end do + + !----------------------! + ! Loop over amplitudes ! + !----------------------! + + write(*,*) + write(*,*)'----------------------------------------------' + write(*,*)'| CC-based G0W0 calculation |' + write(*,*)'----------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','HF','|','G0W0','|','Conv','|' + write(*,*)'----------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Compute residual for 2h1p sector + + r_2h1p = V_2h1p + matmul(C_2h1p,t_2h1p) - t_2h1p*eHF(p) & + - dot_product(t_2h1p,V_2h1p)*t_2h1p - t_2h1p*dot_product(V_2p1h,t_2p1h) + + ! Compute residual for 2p1h sector + + r_2p1h = V_2p1h + matmul(C_2p1h,t_2p1h) - t_2p1h*eHF(p) & + - t_2p1h*dot_product(V_2h1p,t_2h1p) - dot_product(t_2p1h,V_2p1h)*t_2p1h + + ! Check convergence + + Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h))) + + ! Update amplitudes + + t_2h1p(:) = t_2h1p(:) - r_2h1p(:)/delta_2h1p(:) + t_2p1h(:) = t_2p1h(:) - r_2p1h(:)/delta_2p1h(:) + + ! Compute quasiparticle energies + + eGW(p) = eHF(p) + dot_product(V_2h1p,t_2h1p) + dot_product(V_2p1h,t_2p1h) + + ! Renormalization factor + + Z(:) = 1d0 + + ! Dump results + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|' + + end do + + write(*,*)'----------------------------------------------' + !------------------------------------------------------------------------ + ! End of SCF loop + !------------------------------------------------------------------------ + + ! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' CC-G0W0 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)') & + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + 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,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' + write(*,*)'-------------------------------------------------------------------------------' + + end do + +end subroutine diff --git a/src/GW/CCGW.f90 b/src/GW/ccRGW.f90 similarity index 99% rename from src/GW/CCGW.f90 rename to src/GW/ccRGW.f90 index 8b36a4c..c55dbd7 100644 --- a/src/GW/CCGW.f90 +++ b/src/GW/ccRGW.f90 @@ -1,4 +1,4 @@ -subroutine CCGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) +subroutine ccRGW(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! CC-based GW module diff --git a/src/GW/ccRGW_mat.f90 b/src/GW/ccRGW_mat.f90 new file mode 100644 index 0000000..e8077cf --- /dev/null +++ b/src/GW/ccRGW_mat.f90 @@ -0,0 +1,295 @@ +subroutine ccRGW_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + +! CC-based GW module + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + 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) + +! Local variables + + integer :: p,q + integer :: i,j,k,l + integer :: a,b,c,d + + integer :: nSCF + double precision :: Conv + + double precision,allocatable :: delta_2h1p(:,:) + double precision,allocatable :: delta_2p1h(:,:) + + double precision,allocatable :: C_2h1p(:,:) + double precision,allocatable :: C_2p1h(:,:) + + double precision,allocatable :: V_2h1p(:,:) + double precision,allocatable :: V_2p1h(:,:) + + double precision,allocatable :: r_2h1p(:,:) + double precision,allocatable :: r_2p1h(:,:) + + double precision,allocatable :: t_2h1p(:,:) + double precision,allocatable :: t_2p1h(:,:) + + double precision,allocatable :: SigGW(:,:) + double precision,allocatable :: eGW(:) + double precision,allocatable :: F(:,:) + double precision,allocatable :: Z(:) + + double precision,external :: Kronecker_delta + + integer :: n1h,n1p,n1h1p,n2h1p,n2p1h + integer :: ija,klc,iab,kcd + + +! Hello world + + write(*,*) + write(*,*)'*****************************' + write(*,*)'* CC-based G0W0 Calculation *' + write(*,*)'*****************************' + write(*,*) + +! Form energy denominator and guess amplitudes + + n1h = nO + n1p = nV + n1h1p = n1h*n1p + n2h1p = n1h*n1h1p + n2p1h = n1p*n1h1p + + allocate(delta_2h1p(n2h1p,nOrb),delta_2p1h(n2p1h,nOrb)) + allocate(C_2h1p(n2h1p,n2h1p),C_2p1h(n2p1h,n2p1h)) + allocate(V_2h1p(nOrb,n2h1p),V_2p1h(nOrb,n2p1h)) + allocate(t_2h1p(n2h1p,nOrb),t_2p1h(n2p1h,nOrb)) + allocate(r_2h1p(n2h1p,nOrb),r_2p1h(n2p1h,nOrb)) + allocate(F(nOrb,nOrb),eGW(nOrb),SigGW(nOrb,nOrb),Z(nOrb)) + + F(:,:) = 0d0 + do p=nC+1,nOrb-nR + F(p,p) = eHF(p) + end do + + ! Compute C2h1p and C2p1h + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + ija = ija + 1 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nOrb-nR + klc = klc + 1 + + C_2h1p(ija,klc) = ((eHF(i) + eHF(j) - eHF(a))*Kronecker_delta(j,l)*Kronecker_delta(a,c) & + - 2d0*ERI(j,c,a,l))*Kronecker_delta(i,k) + + end do + end do + end do + + end do + end do + end do + + iab = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + iab = iab + 1 + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR + kcd = kcd + 1 + + C_2p1h(iab,kcd) = ((eHF(a) + eHF(b) - eHF(i))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + + 2d0*ERI(a,k,i,c))*Kronecker_delta(b,d) + + end do + end do + end do + + end do + end do + end do + + ! Initialization + + Conv = 1d0 + nSCF = 0 + + t_2h1p(:,:) = 0d0 + t_2p1h(:,:) = 0d0 + + ! Compute energy differences + + ija = 0 + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nOrb-nR + ija = ija + 1 + + do p=nC+1,nOrb-nR + + delta_2h1p(ija,p) = eHF(i) + eHF(j) - eHF(a) - eGW(p) + + end do + + end do + end do + end do + + iab = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR + iab = iab + 1 + + do p=nC+1,nOrb-nR + + delta_2p1h(iab,p) = eHF(a) + eHF(b) - eHF(i) - eGW(p) + + end do + + end do + end do + end do + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nOrb-nR + klc = klc + 1 + + do p=nC+1,nOrb-nR + + V_2h1p(p,klc) = sqrt(2d0)*ERI(p,c,k,l) + + end do + + end do + end do + end do + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR + kcd = kcd + 1 + + do p=nC+1,nOrb-nR + + V_2p1h(p,kcd) = sqrt(2d0)*ERI(p,k,d,c) + + end do + + end do + end do + end do + + !----------------------! + ! Loop over amplitudes ! + !----------------------! + + write(*,*) + write(*,*)'----------------------------------------------' + write(*,*)'| CC-based G0W0 calculation |' + write(*,*)'----------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','HF','|','G0W0','|','Conv','|' + write(*,*)'----------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Compute residual for 2h1p sector + + r_2h1p = transpose(V_2h1p) + matmul(C_2h1p,t_2h1p) - matmul(t_2h1p,F) & + - matmul(matmul(t_2h1p,V_2h1p),t_2h1p) - matmul(matmul(t_2h1p,V_2p1h),t_2p1h) + + ! Compute residual for 2p1h sector + + r_2p1h = transpose(V_2p1h) + matmul(C_2p1h,t_2p1h) - matmul(t_2p1h,F) & + - matmul(matmul(t_2p1h,V_2h1p),t_2h1p) - matmul(matmul(t_2p1h,V_2p1h),t_2p1h) + + ! Check convergence + + Conv = max(maxval(abs(r_2h1p)),maxval(abs(r_2p1h))) + + ! Update amplitudes + + t_2h1p(:,:) = t_2h1p(:,:) - r_2h1p(:,:)/delta_2h1p(:,:) + t_2p1h(:,:) = t_2p1h(:,:) - r_2p1h(:,:)/delta_2p1h(:,:) + + ! Compute quasiparticle energies + + SigGW(:,:) = F(:,:) + matmul(V_2h1p,t_2h1p) + matmul(V_2p1h,t_2p1h) + + call diagonalize_matrix(nOrb,SigGW,eGW) + + ! Renormalization factor + + Z(:) = 0d0 + + ! Dump results + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',eGW(nO)*HaToeV,'|',eGW(nO+1)*HaToeV,'|',Conv,'|' + + end do + + write(*,*)'----------------------------------------------' + !------------------------------------------------------------------------ + ! End of SCF loop + !------------------------------------------------------------------------ + + ! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' CC-G0W0 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)') & + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=nC+1,nOrb-nR + 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,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' + end do + write(*,*)'-------------------------------------------------------------------------------' + +end subroutine diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufRG0W0.f90 similarity index 98% rename from src/GW/ufG0W0.f90 rename to src/GW/ufRG0W0.f90 index 6e9b2d0..6548151 100644 --- a/src/GW/ufG0W0.f90 +++ b/src/GW/ufRG0W0.f90 @@ -1,6 +1,6 @@ -subroutine ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufRG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) -! Unfold G0W0 equations +! Upfolded G0W0 equations implicit none include 'parameters.h' @@ -47,7 +47,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) double precision,allocatable :: XmY(:,:) double precision,allocatable :: rho(:,:,:) - logical :: verbose = .true. + logical :: verbose = .false. double precision,parameter :: cutoff1 = 0.01d0 double precision,parameter :: cutoff2 = 0.01d0 double precision :: eF @@ -372,8 +372,8 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) write(*,*)'-------------------------------------------' do s=1,nH - if(eGW(s) < eF .and. eGW(s) > eF - window) then -! if(Z(s) > cutoff1) then +! if(eGW(s) < eF .and. eGW(s) > eF - window) then + if(Z(s) > cutoff1) then write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|' end if