saving work in ppBSE@GW

src/GW/GGW_ppBSE.f90

@@ -65,22 +65,18 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
 ! Compute RPA screening !
 !-----------------------!
 
-  EcRPA = 0d0
-
   allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
            Aph(nS,nS),Bph(nS,nS))
  
                  call phGLR_A(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
   if(.not.TDA_W) call phGLR_B(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
-  call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  call phGLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
 
-  call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call GGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
 
-  EcBSE = 0d0
-
   nOO = nO*(nO-1)/2
   nVV = nV*(nV-1)/2
 
@@ -103,7 +99,7 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
   Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
   Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
 
-  call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
+  call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE)
 
   call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 

src/GW/GGW_ppBSE_dynamic_kernel_C.f90

@@ -53,30 +53,34 @@ subroutine GGW_ppBSE_dynamic_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,eGW,Om,
 
           do m=1,nS
 
-            dem = OmBSE - eGW(c) - Om(m) - eGW(b)
+            dem = OmBSE - (eGW(a) + eGW(c) + Om(m))
+!           num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
+            num = - rho(b,c,m)*rho(a,d,m)
+
+            KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
+            ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+            dem = OmBSE - (eGW(b) + eGW(d) + Om(m))
+!           num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
+            num = - rho(b,c,m)*rho(a,d,m)
+
+            KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
+            ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
+
+            dem = OmBSE - (eGW(b) + eGW(c) + Om(m))
             num = rho(a,c,m)*rho(b,d,m)
+!           num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
 
             KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
             ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
 
-            dem = OmBSE - eGW(c) - Om(m) - eGW(a)
+            dem = OmBSE - (eGW(a) + eGW(d) + Om(m))
-            num = rho(b,c,m)*rho(a,d,m)
+            num = rho(a,c,m)*rho(b,d,m)
-
+!           num = 0.5d0*(rho(a,c,m)*rho(b,d,m) - rho(b,c,m)*rho(a,d,m))
-            KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
-            ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
-
-            dem = OmBSE - eGW(d) - Om(m) - eGW(a)
-            num = rho(a,c,m)*rho(b,d,m) 
 
             KC_dyn(ab,cd) = KC_dyn(ab,cd) + num*dem/(dem**2 + eta**2)
             ZC_dyn(ab,cd) = ZC_dyn(ab,cd) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
 
-            dem = OmBSE - eGW(d) - Om(m) - eGW(b)
-            num = rho(b,c,m)*rho(a,d,m)
-
-            KC_dyn(ab,cd) = KC_dyn(ab,cd) - num*dem/(dem**2 + eta**2)
-            ZC_dyn(ab,cd) = ZC_dyn(ab,cd) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2
-
           end do
 
         end do

src/GW/GGW_ppBSE_dynamic_perturbation.f90

@@ -71,8 +71,6 @@ subroutine GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,e
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
 
-  print*,nVV,maxVV
-
   do ab=1,min(nVV,maxVV)
 
     if(dTDA) then 

src/GW/GGW_ppBSE_static_kernel_C.f90

@@ -37,14 +37,11 @@ subroutine GGW_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,r
 
   double precision,intent(out)  :: KC(nVV,nVV)
 
-!---------------!
-! SpinOrb block !
-!---------------!
-
   ab = 0
   do a=nO+1,nBas-nR
     do b=a+1,nBas-nR
       ab = ab + 1
+
       cd = 0
       do c=nO+1,nBas-nR
         do d=c+1,nBas-nR

src/GW/GGW_ppBSE_upfolded.f90

@@ -24,7 +24,8 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
   integer                       :: i,j,k,l
   integer                       :: a,b,c,d
   integer                       :: m,ij,kl,ijm
-  integer,parameter             :: maxH = 1000
+  integer                       :: ab,cd,abm
+  integer,parameter             :: maxH = 100
   double precision              :: tmp,tmp1,tmp2,tmp3,tmp4
 
   integer                       :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,n2h2p,nH
@@ -64,7 +65,8 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
   n3h1p = n2h*n1h1p
   n3p1h = n2p*n1h1p
 
-  nH = n2h + n3h1p 
+  nH = n2p + n3p1h 
+! nH = n2h + n3h1p 
 
 ! Memory allocation
 
@@ -91,21 +93,21 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
 !----------------------------------------!
 
   !---------!
-  ! Block D !
+  ! Block C !
   !---------!
 
-  ij = 0
+  ab = 0
-  do i=nC+1,nO
+  do a=nO+1,nOrb-nR
-    do j=i+1,nO
+    do b=a+1,nOrb-nR
-    ij = ij + 1
+    ab = ab + 1
  
-    kl = 0
+    cd = 0
-    do k=nC+1,nO
+    do c=nO+1,nOrb-nR
-      do l=k+1,nO
+      do d=c+1,nOrb-nR
-      kl = kl + 1
+      cd = cd + 1
  
-          H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+          H(ab,cd) = (eGW(a) + eGW(b))*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
-                     + (ERI(i,j,k,l) - ERI(i,j,l,k))
+                     + (ERI(a,b,c,d) - ERI(a,b,d,c))
  
         end do
       end do
@@ -113,29 +115,52 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
     end do
   end do
 
-  !----------------!
+  !---------!
-  ! Blocks M1 & M2 !
+  ! Block D !
-  !----------------!
+  !---------!
 
-  ijm = 0
+! ij = 0
-  do i=nC+1,nO
+! do i=nC+1,nO
-    do j=i+1,nO
+!   do j=i+1,nO
+!   ij = ij + 1
+!
+!   kl = 0
+!   do k=nC+1,nO
+!     do l=k+1,nO
+!     kl = kl + 1
+!
+!         H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+!                    + (ERI(i,j,k,l) - ERI(i,j,l,k))
+!
+!       end do
+!     end do
+!
+!   end do
+! end do
+
+  !-----------------!
+  ! Coupling Blocks !
+  !-----------------!
+
+  abm = 0
+  do a=nO+1,nOrb-nR
+    do b=a+1,nOrb-nR
       do m=1,nS
-        ijm = ijm + 1
+        abm = abm + 1
 
-        kl = 0
+        cd = 0
-        do k=nC+1,nO
+        do c=nO+1,nOrb-nR
-          do l=k+1,nO
+          do d=c+1,nOrb-nR
-            kl = kl + 1
+            cd = cd + 1
 
-            tmp1 = Kronecker_delta(j,l)*rho(i,k,m)
+            tmp1 = Kronecker_delta(b,d)*rho(a,c,m)
-            tmp2 = Kronecker_delta(j,k)*rho(i,l,m)
+            tmp2 = Kronecker_delta(b,c)*rho(a,d,m)
-            tmp3 = Kronecker_delta(i,l)*rho(j,k,m)
+            tmp3 = Kronecker_delta(a,d)*rho(b,c,m)
-            tmp4 = Kronecker_delta(i,k)*rho(j,l,m)
+            tmp4 = Kronecker_delta(a,c)*rho(b,d,m)
+
+            H(n2p+0*n3p1h+abm,cd             ) = tmp1 + tmp2
+            H(cd             ,n2p+0*n3p1h+abm) = tmp3 + tmp4
 
-            H(n2h+0*n3h1p+ijm,kl             ) = tmp1 - tmp2
-            H(kl             ,n2h+0*n3h1p+ijm) = tmp3 - tmp4
-           
 !           H(n2h+1*n3h1p+ijm,kl             ) = +tmp4
 !           H(kl             ,n2h+1*n3h1p+ijm) = +tmp2
 !          
@@ -152,27 +177,80 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
     end do
   end do
 
+! ijm = 0
+! do i=nC+1,nO
+!   do j=i+1,nO
+!     do m=1,nS
+!       ijm = ijm + 1
+
+!       kl = 0
+!       do k=nC+1,nO
+!         do l=k+1,nO
+!           kl = kl + 1
+
+!           tmp1 = Kronecker_delta(j,l)*rho(i,k,m)
+!           tmp2 = Kronecker_delta(j,k)*rho(i,l,m)
+!           tmp3 = Kronecker_delta(i,l)*rho(j,k,m)
+!           tmp4 = Kronecker_delta(i,k)*rho(j,l,m)
+
+!           H(n2h+0*n3h1p+ijm,kl             ) = tmp1 - tmp2
+!           H(kl             ,n2h+0*n3h1p+ijm) = tmp3 - tmp4
+!          
+!           H(n2h+1*n3h1p+ijm,kl             ) = +tmp4
+!           H(kl             ,n2h+1*n3h1p+ijm) = +tmp2
+!          
+!           H(n2h+2*n3h1p+ijm,kl             ) = +tmp1
+!           H(kl             ,n2h+2*n3h1p+ijm) = +tmp4
+
+!           H(n2h+3*n3h1p+ijm,kl             ) = +tmp3
+!           H(kl             ,n2h+3*n3h1p+ijm) = +tmp1
+
+!         end do
+!       end do
+
+!     end do
+!   end do
+! end do
+
   !------------!
-  ! Block 3h1p !
+  ! Block 3p1h !
   !------------!
 
-  ijm = 0
+  abm = 0
-  do i=nC+1,nO
+  do a=nO+1,nOrb-nR
-    do j=i+1,nO
+    do b=a+1,nOrb-nR
       do m=1,nS
-        ijm = ijm + 1
+        abm = abm + 1
 
-        tmp = - eGW(i) - eGW(j) + Om(m)
+        tmp = eGW(a) + eGW(b) + Om(m)
 
-        H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp
+        H(n2p+0*n3p1h+abm,n2p+0*n3p1h+abm) = tmp
-!       H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp
-!       H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp
-!       H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp
 
       end do
     end do
   end do
 
+  !------------!
+  ! Block 3h1p !
+  !------------!
+
+! ijm = 0
+! do i=nC+1,nO
+!   do j=i+1,nO
+!     do m=1,nS
+!       ijm = ijm + 1
+
+!       tmp = - eGW(i) - eGW(j) + Om(m)
+
+!       H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp
+!       H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp
+!       H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp
+!       H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp
+
+!     end do
+!   end do
+! end do
+
 !-------------------------!
 ! Diagonalize supermatrix !
 !-------------------------!
@@ -211,7 +289,7 @@ subroutine GGW_ppBSE_upfolded(nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
   do s=1,min(nH,maxH)
     if(Z(s) > 1d-7) &
       write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-      '|',s,'|',-OmBSE(s)*HaToeV,'|',Z(s),'|'
+      '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|'
   end do
 
   write(*,*)'-------------------------------------------'