adding ufGW to methods

input/methods

@@ -12,8 +12,8 @@
   F   F    F    
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F     F      F    F
-# G0W0* evGW* qsGW*
+# G0W0* evGW* qsGW* ufG0W0 ufGW
-  T    F    F
+  T    F    F       T      T
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

src/MBPT/ufG0W0.f90

@@ -0,0 +1,213 @@
+subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
+
+! Unfold G0W0 equations
+
+  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
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: ERHF
+  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eHF(nBas)
+
+! Local variables
+
+  integer                       :: p
+  integer                       :: s
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+  integer                       :: klc,kcd,ija,iab
+
+  integer                       :: n2h1p,n2p1h,nH
+  double precision,external     :: Kronecker_delta
+  double precision,allocatable  :: H(:,:)
+  double precision,allocatable  :: cGW(:,:)
+  double precision,allocatable  :: eGW(:)
+  double precision,allocatable  :: Z(:)
+
+! Output variables
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'**********************************************'
+  write(*,*)'|        Unfolded G0W0 calculation           |'
+  write(*,*)'**********************************************'
+  write(*,*)
+
+! TDA for W
+
+  write(*,*) 'Tamm-Dancoff approximation for dynamic screening by default!'
+  write(*,*)
+
+! Dimension of the supermatrix
+
+  n2h1p = nO*nO*nS
+  n2p1h = nV*nV*nO
+  nH = 1 + n2h1p + n2p1h
+
+! Memory allocation
+
+  allocate(H(nH,nH),cGW(nH,nH),eGW(nH),Z(nH))
+
+! Initialization
+
+  H(:,:) = 0d0
+
+!---------------------------!
+!  Compute GW supermatrix   !
+!---------------------------!
+!                           !
+!     |   F   V2h1p V2p1h | ! 
+!     |                   | ! 
+! H = | V2h1p C2h1p     0 | ! 
+!     |                   | ! 
+!     | V2p1h   0   C2p1h | ! 
+!                           !
+!---------------------------!
+
+  !-------------!
+  ! Block C2h1p !
+  !-------------!
+
+  ija = 0
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        ija = ija + 1
+
+        klc = 0
+        do k=nC+1,nO
+          do l=nC+1,nO
+            do c=nO+1,nBas-nR
+              klc = klc + 1
+
+              H(1+ija,1+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
+
+  !-------------!
+  ! Block C2p1h !
+  !-------------!
+
+  iab = 0
+  do i=nC+1,nO
+    do a=nO+1,nBas-nR
+      do b=nO+1,nBas-nR
+        iab = iab + 1
+
+        kcd = 0
+        do k=nC+1,nO
+          do c=nO+1,nBas-nR
+            do d=nO+1,nBas-nR
+              kcd = kcd + 1
+
+              H(1+n2h1p+iab,1+n2h1p+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
+
+  do p=nC+1,nBas
+
+  !---------!
+  ! Block F !
+  !---------!
+
+    H(1,1) = eHF(p)
+
+  !-------------!
+  ! Block V2h1p !
+  !-------------!
+
+    klc = 0
+    do k=nC+1,nO
+      do l=nC+1,nO
+        do c=nO+1,nBas-nR
+          klc = klc + 1
+
+          H(1       ,1+klc) = sqrt(2d0)*ERI(p,c,k,l)
+          H(1+klc,1       ) = sqrt(2d0)*ERI(p,c,k,l)
+
+        end do
+      end do
+    end do
+
+  !-------------!
+  ! Block V2p1h !
+  !-------------!
+
+    kcd = 0
+    do k=nC+1,nO
+      do c=nO+1,nBas-nR
+        do d=nO+1,nBas-nR
+          kcd = kcd + 1
+
+          H(1              ,1+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
+          H(1+n2h1p+kcd,1              ) = sqrt(2d0)*ERI(p,k,d,c)
+
+        end do
+      end do
+    end do
+
+  !-------------------------!
+  ! Diagonalize supermatrix !
+  !-------------------------!
+
+   cGW(:,:) = H(:,:)
+   call diagonalize_matrix(nH,cGW,eGW)
+
+  !-----------------!
+  ! Compute weights !
+  !-----------------!
+
+    Z(:) = 0d0
+    do s=1,nH
+      Z(s) = Z(s) + cGW(1,s)**2
+    end do
+
+  !--------------!
+  ! Dump results !
+  !--------------!
+
+    write(*,*)'-------------------------------------------'
+    write(*,'(A35,I3)')' G0W0 energies (eV) for orbital ',p
+    write(*,*)'-------------------------------------------'
+    write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
+              '|','#','|','e_QP (eV)','|','Z','|'
+    write(*,*)'-------------------------------------------'
+ 
+    do s=1,nH
+      write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+      '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
+    enddo
+ 
+    write(*,*)'-------------------------------------------'
+    write(*,*)
+
+  end do
+
+end subroutine ufG0W0

src/MBPT/ufGW.f90

@@ -115,7 +115,7 @@ subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
         do d=nO+1,nBas-nR
           kcd = kcd + 1
 
-          H(p              ,nBas+n2h1P+kcd) = sqrt(2d0)*ERI(p,k,d,c)
+          H(p              ,nBas+n2h1p+kcd) = sqrt(2d0)*ERI(p,k,d,c)
           H(nBas+n2h1p+kcd,p              ) = sqrt(2d0)*ERI(p,k,d,c)
 
         end do
@@ -201,24 +201,21 @@ subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
 ! Dump results !
 !--------------!
 
-  write(*,*) '--------------------------------------------'
-  write(*,*) ' GW supermatrix quasiparticle energies (eV) '
-  write(*,*) '--------------------------------------------'
-  write(*,*) 
-  call matout(nH,1,HaToeV*eGW(:))
-  write(*,*) 
-  write(*,*) '-----------------------'
-  write(*,*) ' Quasiparticle weights '
-  write(*,*) '-----------------------'
-  write(*,*) 
-  call matout(nH,1,Z(:))
 
-! write(*,*) 
+
-! write(*,*) '-------------------------'
+  write(*,*)'-------------------------------------------'
-! write(*,*) ' GW supermatrix orbitals '
+  write(*,*)'        unfolded GW energies (eV)          '
-! write(*,*) '-------------------------'
+  write(*,*)'-------------------------------------------'
-! write(*,*) 
+  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
-! call matout(nH,nH,H)
+            '|','#','|','e_QP (eV)','|','Z','|'
-! write(*,*) 
+  write(*,*)'-------------------------------------------'
+
+  do s=1,nH
+    write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
+    '|',s,'|',eGW(s)*HaToeV,'|',Z(s),'|'
+  enddo
+
+  write(*,*)'-------------------------------------------'
+  write(*,*)
 
 end subroutine ufGW

src/QuAcK/QuAcK.f90

@@ -15,7 +15,7 @@ program QuAcK
   logical                       :: doRPA,doRPAx,doppRPA
   logical                       :: doADC
   logical                       :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
-  logical                       :: doG0W0,doevGW,doqsGW,doufGW
+  logical                       :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical                       :: doG0T0,doevGT,doqsGT
   logical                       :: doMCMP2,doMinMCMP2
   logical                       :: doGTGW = .false.
@@ -169,6 +169,7 @@ program QuAcK
                     doG0F2,doevGF2,doqsGF2,           & 
                     doG0F3,doevGF3,                   &
                     doG0W0,doevGW,doqsGW,             &
+                    doufG0W0,doufGW,                  &
                     doG0T0,doevGT,doqsGT,             &
                     doMCMP2)
 
@@ -1025,10 +1026,24 @@ program QuAcK
   end if
 
 !------------------------------------------------------------------------
-! Perform ufGW calculatiom
+! Perform ufG0W0 calculatiom
 !------------------------------------------------------------------------
 
-  doufGW = .true.
+  if(doufG0W0) then
+    
+    call cpu_time(start_ufGW)
+    call ufG0W0(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
+    call cpu_time(end_ufGW)
+  
+    t_ufGW = end_ufGW - start_ufGW
+    write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_ufGW,' seconds'
+    write(*,*)
+
+  end if
+
+!------------------------------------------------------------------------
+! Perform ufGW calculatiom
+!------------------------------------------------------------------------
 
   if(doufGW) then
     

src/QuAcK/read_methods.f90

@@ -7,6 +7,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
                         doG0F2,doevGF2,doqsGF2,           &
                         doG0F3,doevGF3,                   & 
                         doG0W0,doevGW,doqsGW,             & 
+                        doufG0W0,doufGW,                  & 
                         doG0T0,doevGT,doqsGT,             & 
                         doMCMP2)
 
@@ -23,7 +24,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
   logical,intent(out)           :: doCIS,doCIS_D,doCID,doCISD,doFCI
   logical,intent(out)           :: doRPA,doRPAx,doppRPA
   logical,intent(out)           :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3  
-  logical,intent(out)           :: doG0W0,doevGW,doqsGW
+  logical,intent(out)           :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
   logical,intent(out)           :: doG0T0,doevGT,doqsGT
   logical,intent(out)           :: doMCMP2
 
@@ -73,12 +74,14 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
   doevGF3 = .false.
 
   doG0W0   = .false.
-  doevGT = .false.
-  doqsGT = .false.
-
-  doG0T0 = .false.
   doevGW   = .false.
   doqsGW   = .false.
+  doufG0W0 = .false.
+  doufGW   = .false.
+
+  doG0T0 = .false.
+  doevGT = .false.
+  doqsGT = .false.
 
   doMCMP2 = .false.
 
@@ -146,10 +149,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM,           &
 ! Read GW methods
 
   read(1,*) 
-  read(1,*) answer1,answer2,answer3
+  read(1,*) answer1,answer2,answer3,answer4,answer5
   if(answer1 == 'T') doG0W0   = .true.
   if(answer2 == 'T') doevGW   = .true.
   if(answer3 == 'T') doqsGW   = .true.
+  if(answer4 == 'T') doufG0W0 = .true.
+  if(answer5 == 'T') doufGW   = .true.
 
 ! Read GT methods