Merge branch 'master' of github.com:pfloos/QuAcK

src/GT/RG0T0pp.f90

@@ -64,6 +64,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   double precision,allocatable  :: Z(:)
   double precision,allocatable  :: eGT(:)
   double precision,allocatable  :: eGTlin(:)
+  integer,          allocatable :: supp_data_int(:)
+  double precision, allocatable :: supp_data_dbl(:)
   double precision, allocatable :: Om(:), R(:,:)
 
 
@@ -94,11 +96,11 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
 ! Dimensions of the pp-RPA linear reponse matrices
 
-  nOOs = nO*(nO + 1)/2
+  !nOOs = nO*(nO + 1)/2
-  nVVs = nV*(nV + 1)/2
+  !nVVs = nV*(nV + 1)/2
 
-  !nOOs = nO*nO
+  nOOs = nO*nO
-  !nVVs = nV*nV
+  nVVs = nV*nV
 
   nOOt = nO*(nO - 1)/2
   nVVt = nV*(nV - 1)/2
@@ -118,8 +120,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 !----------------------------------------------
 
   ispin  = 1
-  iblock = 1
+  !iblock = 1
-  !iblock = 3
+  iblock = 3
 
 ! Compute linear response
 
@@ -138,7 +140,19 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   !n_states = nOOs + 5
   !n_states_diag = n_states + 4
   !allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag))
-  !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag, "RPA")
+  !allocate(supp_data_dbl(1), supp_data_int(1))
+  !supp_data_int(1) = 0
+  !supp_data_dbl(1) = 0.d0
+  !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, &
+  !                   1.d0,                                        & ! lambda
+  !                   eHF(1),                                      &
+  !                   0.d0,                                        & ! eF
+  !                   ERI(1,1,1,1),                                &
+  !                   supp_data_int(1), 1,                         &
+  !                   supp_data_dbl(1), 1,                         &
+  !                   Om(1), R(1,1), n_states, n_states_diag, "RPA")
+  !deallocate(supp_data_dbl, supp_data_int)
+  !deallocate(Om, R)
   !stop
 
   if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s)
@@ -169,7 +183,19 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   !n_states = nOOt + 5
   !n_states_diag = n_states + 4
   !allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag))
-  !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag, "RPA")
+  !allocate(supp_data_dbl(1), supp_data_int(1))
+  !supp_data_int(1) = 0
+  !supp_data_dbl(1) = 0.d0
+  !call ppLR_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, &
+  !                   1.d0,                                        & ! lambda
+  !                   eHF(1),                                      &
+  !                   0.d0,                                        & ! eF
+  !                   ERI(1,1,1,1),                                &
+  !                   supp_data_int(1), 1,                         &
+  !                   supp_data_dbl(1), 1,                         &
+  !                   Om(1), R(1,1), n_states, n_states_diag, "RPA")
+  !deallocate(Om, R)
+  !deallocate(supp_data_dbl)
   !stop
 
   if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)

src/GW/RGW_ppBSE.f90

@@ -1,4 +1,4 @@
-subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
+subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies at the pp level
 
@@ -15,16 +15,16 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
   logical,intent(in)            :: triplet
 
   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)   :: eW(nBas)
+  double precision,intent(in)   :: eW(nOrb)
-  double precision,intent(in)   :: eGW(nBas)
+  double precision,intent(in)   :: eGW(nOrb)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
 
 ! Local variables
 
@@ -37,6 +37,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
   integer                       :: nOO
   integer                       :: nVV
 
+  integer                       :: p, q, m
+  integer                       :: i_data, supp_data_dbl_size, supp_data_int_size
+  integer                       :: n_states, n_states_diag
+
   double precision,allocatable  :: Aph(:,:)
   double precision,allocatable  :: Bph(:,:)
 
@@ -62,6 +66,10 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
   double precision,allocatable  :: KC_sta(:,:)
   double precision,allocatable  :: KD_sta(:,:)
   
+  integer,         allocatable  :: supp_data_int(:)
+  double precision,allocatable  :: supp_data_dbl(:)
+  double precision,allocatable  :: Om(:), R(:,:)
+
 ! Output variables
 
   double precision,intent(out)  :: EcBSE(nspin)
@@ -74,16 +82,16 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
   isp_W = 1
   EcRPA = 0d0
 
-  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
+  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
            Aph(nS,nS),Bph(nS,nS))
  
-  call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+  call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
 
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  if(.not.TDA_W) call phLR_B(isp_W,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 RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
   deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
 
@@ -104,42 +112,100 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
     nOO = nO*(nO+1)/2
     nVV = nV*(nV+1)/2
 
-    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),       &
+    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV))
-             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),       &
+    allocate(Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
-             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
-             KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
 
     ! Compute BSE excitation energies
 
-    call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+    ! ---
-    call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    ! LAPACK
-    if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    ! ---
+
+    allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
+    allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
+
+    call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+    call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
     
-    call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+    call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
-    call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+    call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
-    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+    if(.not.TDA) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
 
     Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
     Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
     Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
 
     call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+    deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
 
-    ! TODO
+    print*, 'LAPACK:'
-    !call ppLR_RGW_ppBSE_davidson(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+    print*, Om2
+    print*, Om1
 
-    call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+    ! ---
+
+
+
+    ! ---
+    ! Davidson
+    ! ---
+
+    n_states = nOO + 5
+    n_states_diag = n_states + 4
+    allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag))
+
+    supp_data_int = 1
+    allocate(supp_data_int(supp_data_int_size))
+    supp_data_int(1) = nS
+
+    supp_data_dbl_size = nS + nOrb*nOrb*nS + 1
+    allocate(supp_data_dbl(supp_data_dbl_size))
+    ! scalars
+    supp_data_dbl(1) = eta
+    i_data = 1
+    ! rho_RPA
+    do m = 1, nS
+      do q = 1, nOrb
+        do p = 1, nOrb
+          i_data = i_data + 1
+          supp_data_dbl(i_data) = rho_RPA(p,q,m)
+        enddo
+      enddo
+    enddo
+    ! OmRPA
+    do m = 1, nS
+      i_data = i_data + 1
+      supp_data_dbl(i_data) = OmRPA(m)
+    enddo
+
+    call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, &
+                       1.d0,                                   & ! lambda
+                       eGW(1),                                 &
+                       0.d0,                                   & ! eF
+                       ERI(1,1,1,1),                           &
+                       supp_data_int(1), supp_data_int_size,   &
+                       supp_data_dbl(1), supp_data_dbl_size,   &
+                       Om(1), R(1,1), n_states, n_states_diag, "GW")
+
+    deallocate(Om, R)
+    deallocate(supp_data_dbl, supp_data_int)
+    stop
+
+    ! ---
+
+    call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
     !----------------------------------------------------!
     ! Compute the dynamical screening at the ppBSE level !
     !----------------------------------------------------!
 
     if(dBSE) &
-        call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
+        call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2)
   end if
 
 !-------------------
@@ -159,40 +225,97 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
     nOO = nO*(nO-1)/2
     nVV = nV*(nV-1)/2
 
-    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),       &
+    allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV))
-             Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO),       &
+    allocate(Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO))
-             Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
-             KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
 
     ! Compute BSE excitation energies
 
-    call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+    ! ---
-    call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    ! LAPACK
-    if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    ! ---
 
-    call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+    allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO))
-    call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+    allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO))
-    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+
+    call RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+    call RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+
+    call ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+    call ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+    if(.not.TDA) call ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
 
     Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
     Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
     Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
 
     call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+    deallocate(Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
+
+    !print*, 'LAPACK:'
+    !print*, Om2
+    !print*, Om1
+
+    ! ---
+    ! Davidson
+    ! ---
+
+    !n_states = nOO + 5
+    !n_states_diag = n_states + 4
+    !allocate(Om(nOO+nVV), R(nOO+nVV,n_states_diag))
+
+    !supp_data_int = 1
+    !allocate(supp_data_int(supp_data_int_size))
+    !supp_data_int(1) = nS
+
+    !supp_data_dbl_size = nS + nOrb*nOrb*nS + 1
+    !allocate(supp_data_dbl(supp_data_dbl_size))
+    !! scalars
+    !supp_data_dbl(1) = eta
+    !i_data = 1
+    !! rho_RPA
+    !do m = 1, nS
+    !  do q = 1, nOrb
+    !    do p = 1, nOrb
+    !      i_data = i_data + 1
+    !      supp_data_dbl(i_data) = rho_RPA(p,q,m)
+    !    enddo
+    !  enddo
+    !enddo
+    !! OmRPA
+    !do m = 1, nS
+    !  i_data = i_data + 1
+    !  supp_data_dbl(i_data) = OmRPA(m)
+    !enddo
+
+    !call ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, &
+    !                   1.d0,                                   & ! lambda
+    !                   eGW(1),                                 &
+    !                   0.d0,                                   & ! eF
+    !                   ERI(1,1,1,1),                           &
+    !                   supp_data_int(1), supp_data_int_size,   &
+    !                   supp_data_dbl(1), supp_data_dbl_size,   &
+    !                   Om(1), R(1,1), n_states, n_states_diag, "GW")
+
+    !deallocate(Om, R)
+    !deallocate(supp_data_dbl, supp_data_int)
+    !stop
+
+    ! ---
 
     EcBSE(ispin) = 3d0*EcBSE(ispin)
 
-    call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+    call ppLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
     !----------------------------------------------------!
     ! Compute the dynamical screening at the ppBSE level !
     !----------------------------------------------------!
 
     if(dBSE) &
-        call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
+        call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
-    deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
+    deallocate(Om1,X1,Y1,Om2,X2,Y2)
 
   end if
 

src/GW/RGW_ppBSE_static_kernel_B.f90

@@ -1,4 +1,4 @@
-subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB)
+subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB)
 
 ! Compute the VVOO block of the static screening W for the pp-BSE
 
@@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
 ! Input variables
 
   integer,intent(in)            :: ispin
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -18,9 +18,9 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
   double precision,intent(in)   :: Om(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
 
 ! Local variables
 
@@ -44,8 +44,8 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
   if(ispin == 1) then
 
     ab = 0
-    do a=nO+1,nBas-nR
+    do a=nO+1,nOrb-nR
-      do b=a,nBas-nR
+      do b=a,nOrb-nR
         ab = ab + 1
         ij = 0
         do i=nC+1,nO
@@ -75,8 +75,8 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
   if(ispin == 2) then
 
     ab = 0
-    do a=nO+1,nBas-nR
+    do a=nO+1,nOrb-nR
-     do b=a+1,nBas-nR
+     do b=a+1,nOrb-nR
         ab = ab + 1
         ij = 0
         do i=nC+1,nO

src/GW/RGW_ppBSE_static_kernel_C.f90

@@ -1,4 +1,4 @@
-subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC)
+subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nOrb,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC)
 
 ! Compute the VVVV block of the static screening W for the pp-BSE
 
@@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 ! Input variables
 
   integer,intent(in)            :: ispin
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -17,9 +17,9 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
   double precision,intent(in)   :: Om(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
 
 ! Local variables
 
@@ -44,21 +44,21 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 
   if(ispin == 1) then
 
-    a0 = nBas - nR - nO
+    a0 = nOrb - nR - nO
     lambda4 = 4.d0 * lambda
     eta2 = eta * eta
 
-    allocate(tmp_m(nBas,nBas,nS))
+    allocate(tmp_m(nOrb,nOrb,nS))
-    allocate(tmp(nBas,nBas,nBas,nBas))
+    allocate(tmp(nOrb,nOrb,nOrb,nOrb))
 
     !$OMP PARALLEL DEFAULT(NONE)         &
     !$OMP          PRIVATE(m, c, a, eps) &
-    !$OMP          SHARED(nS, nBas, eta2, Om, rho, tmp_m)
+    !$OMP          SHARED(nS, nOrb, eta2, Om, rho, tmp_m)
     !$OMP DO
     do m = 1, nS
       eps = Om(m) / (Om(m)*Om(m) + eta2)
-      do c = 1, nBas
+      do c = 1, nOrb
-        do a = 1, nBas
+        do a = 1, nOrb
           tmp_m(a,c,m) = eps * rho(a,c,m)
         enddo
       enddo
@@ -66,19 +66,19 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
     !$OMP END DO
     !$OMP END PARALLEL
 
-    call dgemm("N", "T", nBas*nBas, nBas*nBas, nS, 1.d0,       &
+    call dgemm("N", "T", nOrb*nOrb, nOrb*nOrb, nS, 1.d0,       &
-               tmp_m(1,1,1), nBas*nBas, rho(1,1,1), nBas*nBas, &
+               tmp_m(1,1,1), nOrb*nOrb, rho(1,1,1), nOrb*nOrb, &
-               0.d0, tmp(1,1,1,1), nBas*nBas)
+               0.d0, tmp(1,1,1,1), nOrb*nOrb)
 
     deallocate(tmp_m)
 
     !$OMP PARALLEL DEFAULT(NONE)                           &
     !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, tmp_ab) &
-    !$OMP          SHARED(nO, nBas, nR, nS, a0, lambda4, tmp, KC)
+    !$OMP          SHARED(nO, nOrb, nR, nS, a0, lambda4, tmp, KC)
     !$OMP DO
-    do a = nO+1, nBas-nR
+    do a = nO+1, nOrb-nR
       aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO
-      do b = a, nBas-nR
+      do b = a, nOrb-nR
         ab = aa + b
 
         tmp_ab = lambda4
@@ -87,8 +87,8 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
         endif
 
         cd = 0
-        do c = nO+1, nBas-nR
+        do c = nO+1, nOrb-nR
-          do d = c, nBas-nR
+          do d = c, nOrb-nR
             cd = cd + 1
 
             KC(ab,cd) = -tmp_ab * (tmp(a,c,b,d) + tmp(a,d,b,c))
@@ -105,12 +105,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
     deallocate(tmp)
 
 
-!    do a=nO+1,nBas-nR
+!    do a=nO+1,nOrb-nR
-!      do b=a,nBas-nR
+!      do b=a,nOrb-nR
 !        ab = ab + 1
 !        cd = 0
-!        do c=nO+1,nBas-nR
+!        do c=nO+1,nOrb-nR
-!          do d=c,nBas-nR
+!          do d=c,nOrb-nR
 !            cd = cd + 1
 !
 !              chi = 0d0
@@ -125,7 +125,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !    OpenMP implementation
 !    --- --- ---
 !
-!    a0 = nBas - nR - nO
+!    a0 = nOrb - nR - nO
 !    lambda4 = 4.d0 * lambda
 !    eta2 = eta * eta
 !
@@ -141,11 +141,11 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !
 !    !$OMP PARALLEL DEFAULT(NONE)                              &
 !    !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, m, tmp_ab) &
-!    !$OMP          SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC)
+!    !$OMP          SHARED(nO, nOrb, nR, nS, a0, lambda4, Om_tmp, rho, KC)
 !    !$OMP DO
-!    do a = nO+1, nBas-nR
+!    do a = nO+1, nOrb-nR
 !      aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO
-!      do b = a, nBas-nR
+!      do b = a, nOrb-nR
 !        ab = aa + b
 !
 !        tmp_ab = lambda4
@@ -154,8 +154,8 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !        endif
 !
 !        cd = 0
-!        do c = nO+1, nBas-nR
+!        do c = nO+1, nOrb-nR
-!          do d = c, nBas-nR
+!          do d = c, nOrb-nR
 !            cd = cd + 1
 !
 !            KC(ab,cd) = 0d0
@@ -184,12 +184,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !    --- --- ---
 !
 !    ab = 0
-!    do a=nO+1,nBas-nR
+!    do a=nO+1,nOrb-nR
-!      do b=a,nBas-nR
+!      do b=a,nOrb-nR
 !        ab = ab + 1
 !        cd = 0
-!        do c=nO+1,nBas-nR
+!        do c=nO+1,nOrb-nR
-!          do d=c,nBas-nR
+!          do d=c,nOrb-nR
 !            cd = cd + 1
 !
 !              chi = 0d0
@@ -215,21 +215,21 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 
   if(ispin == 2) then
 
-    a0 = nBas - nR - nO - 1
+    a0 = nOrb - nR - nO - 1
     lambda4 = 4.d0 * lambda
     eta2 = eta * eta
 
-    allocate(tmp_m(nBas,nBas,nS))
+    allocate(tmp_m(nOrb,nOrb,nS))
-    allocate(tmp(nBas,nBas,nBas,nBas))
+    allocate(tmp(nOrb,nOrb,nOrb,nOrb))
 
     !$OMP PARALLEL DEFAULT(NONE)         &
     !$OMP          PRIVATE(m, c, a, eps) &
-    !$OMP          SHARED(nS, nBas, eta2, Om, rho, tmp_m)
+    !$OMP          SHARED(nS, nOrb, eta2, Om, rho, tmp_m)
     !$OMP DO
     do m = 1, nS
       eps = Om(m) / (Om(m)*Om(m) + eta2)
-      do c = 1, nBas
+      do c = 1, nOrb
-        do a = 1, nBas
+        do a = 1, nOrb
           tmp_m(a,c,m) = eps * rho(a,c,m)
         enddo
       enddo
@@ -237,24 +237,24 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
     !$OMP END DO
     !$OMP END PARALLEL
 
-    call dgemm("N", "T", nBas*nBas, nBas*nBas, nS, 1.d0,       &
+    call dgemm("N", "T", nOrb*nOrb, nOrb*nOrb, nS, 1.d0,       &
-               tmp_m(1,1,1), nBas*nBas, rho(1,1,1), nBas*nBas, &
+               tmp_m(1,1,1), nOrb*nOrb, rho(1,1,1), nOrb*nOrb, &
-               0.d0, tmp(1,1,1,1), nBas*nBas)
+               0.d0, tmp(1,1,1,1), nOrb*nOrb)
 
     deallocate(tmp_m)
 
     !$OMP PARALLEL DEFAULT(NONE)                   &
     !$OMP          PRIVATE(a, b, aa, ab, c, d, cd) &
-    !$OMP          SHARED(nO, nBas, nR, nS, a0, lambda4, tmp, KC)
+    !$OMP          SHARED(nO, nOrb, nR, nS, a0, lambda4, tmp, KC)
     !$OMP DO
-    do a = nO+1, nBas-nR
+    do a = nO+1, nOrb-nR
       aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1
-      do b = a+1, nBas-nR
+      do b = a+1, nOrb-nR
         ab = aa + b
 
         cd = 0
-        do c = nO+1, nBas-nR
+        do c = nO+1, nOrb-nR
-          do d = c+1, nBas-nR
+          do d = c+1, nOrb-nR
             cd = cd + 1
 
             KC(ab,cd) = lambda4 * (-tmp(a,c,b,d) + tmp(a,d,b,c))
@@ -272,7 +272,7 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !    OpenMP implementation
 !    --- --- ---
 !
-!    a0 = nBas - nR - nO - 1
+!    a0 = nOrb - nR - nO - 1
 !    lambda4 = 4.d0 * lambda
 !    eta2 = eta * eta
 !
@@ -288,16 +288,16 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !
 !    !$OMP PARALLEL DEFAULT(NONE)                      &
 !    !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, m) &
-!    !$OMP          SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC)
+!    !$OMP          SHARED(nO, nOrb, nR, nS, a0, lambda4, Om_tmp, rho, KC)
 !    !$OMP DO
-!    do a = nO+1, nBas-nR
+!    do a = nO+1, nOrb-nR
 !      aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1
-!      do b = a+1, nBas-nR
+!      do b = a+1, nOrb-nR
 !        ab = aa + b
 !
 !        cd = 0
-!        do c = nO+1, nBas-nR
+!        do c = nO+1, nOrb-nR
-!          do d = c+1, nBas-nR
+!          do d = c+1, nOrb-nR
 !            cd = cd + 1
 !
 !            KC(ab,cd) = 0d0
@@ -323,12 +323,12 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
 !    --- --- ---
 !
 !    ab = 0
-!    do a=nO+1,nBas-nR
+!    do a=nO+1,nOrb-nR
-!      do b=a+1,nBas-nR
+!      do b=a+1,nOrb-nR
 !        ab = ab + 1
 !        cd = 0
-!        do c=nO+1,nBas-nR
+!        do c=nO+1,nOrb-nR
-!          do d=c+1,nBas-nR
+!          do d=c+1,nOrb-nR
 !            cd = cd + 1
 !
 !            chi = 0d0

src/GW/RGW_ppBSE_static_kernel_D.f90

@@ -1,4 +1,4 @@
-subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD)
+subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nOrb,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD)
 
 ! Compute the OOOO block of the static screening W for the pp-BSE
 
@@ -8,7 +8,7 @@ subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ER
 ! Input variables
 
   integer,intent(in)            :: ispin
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -17,9 +17,9 @@ subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ER
   integer,intent(in)            :: nOO
   double precision,intent(in)   :: eta
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
   double precision,intent(in)   :: Om(nS)
-  double precision,intent(in)   :: rho(nBas,nBas,nS)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
 
 ! Local variables
 

src/LR/ppLR.f90

@@ -44,7 +44,6 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA)
   double precision, external    :: trace_matrix
 
 
-
   N = nOO + nVV
 
   allocate(M(N,N), Z(N,N), Om(N))

src/LR/ppLR_GW_davidson.f90

@@ -0,0 +1,402 @@
+
+
+! ---
+
+subroutine ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF, n_states_diag, &
+                           ERI, eta, rho, Om, U, W)
+
+
+  implicit none
+
+  integer,          intent(in)  :: ispin
+  integer,          intent(in)  :: n_states_diag
+  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR, nS
+  double precision, intent(in)  :: lambda, eF, eta
+  double precision, intent(in)  :: e(nOrb)
+  double precision, intent(in)  :: ERI(nOrb,nOrb,nOrb,nOrb)
+  double precision, intent(in)  :: rho(nOrb,nOrb,nS), Om(nS)
+  double precision, intent(in)  :: U(nOO+nVV,n_states_diag)
+  double precision, intent(out) :: W(nOO+nVV,n_states_diag)
+
+  integer                       :: i, j, ij, k, l, kl
+  integer                       :: a, b, c, d, ab, cd
+  integer                       :: m
+  integer                       :: state
+  double precision              :: mat_tmp, chi, eps
+
+  double precision, external    :: Kronecker_delta
+
+  if(ispin .eq. 1) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = a, nOrb-nR
+        ab = ab + 1
+  
+        do state = 1, n_states_diag
+  
+          W(ab,state) = 0.d0
+
+          cd = 0
+          do c = nO+1, nOrb-nR
+            do d = c, nOrb-nR
+              cd = cd + 1
+    
+              mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) &
+                      + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                                      * (1.d0 + Kronecker_delta(c, d)))
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps &
+                          - rho(a,d,m) * rho(b,c,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                             * (1.d0 + Kronecker_delta(c, d)))
+
+              W(ab,state) = W(ab,state) + mat_tmp * U(cd,state)
+            enddo
+          enddo
+
+          ij = nVV
+          do i = nC+1, nO
+            do j = i, nO
+              ij = ij + 1
+  
+              mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                                      * (1.d0 + Kronecker_delta(i, j)))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps &
+                          - rho(i,b,m) * rho(a,j,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                             * (1.d0 + Kronecker_delta(i, j)))
+  
+              W(ab,state) = W(ab,state) - mat_tmp * U(ij,state)
+            enddo
+          enddo
+
+        enddo ! state
+      enddo ! b
+    enddo ! a
+
+    ! ---
+  
+    ij = nVV
+    do i = nC+1, nO
+      do j = i, nO
+        ij = ij + 1
+  
+        do state = 1, n_states_diag
+  
+          W(ij,state) = 0.d0
+  
+          ab = 0
+          do a = nO+1, nOrb-nR
+            do b = a, nOrb-nR
+              ab = ab + 1
+    
+              mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                                      * (1.d0 + Kronecker_delta(i, j)))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps &
+                          - rho(i,b,m) * rho(a,j,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                             * (1.d0 + Kronecker_delta(i, j)))
+
+              W(ij,state) = W(ij,state) + mat_tmp * U(ab,state)
+            enddo
+          enddo
+
+          kl = nVV
+          do k = nC+1, nO
+            do l = k, nO
+              kl = kl + 1
+    
+              mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) &
+                      + lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt( (1.d0 + Kronecker_delta(i, j)) &
+                                                                      * (1.d0 + Kronecker_delta(k, l)))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps &
+                          - rho(i,l,m) * rho(j,k,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(i, j)) &
+                                                             * (1.d0 + Kronecker_delta(k, l)))
+
+              W(ij,state) = W(ij,state) - mat_tmp * U(kl,state)
+            enddo
+          enddo
+  
+        enddo ! state
+      enddo ! j
+    enddo ! i
+
+
+  elseif(ispin .eq. 2) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = a+1, nOrb-nR
+        ab = ab + 1
+  
+        do state = 1, n_states_diag
+  
+          W(ab,state) = 0.d0
+
+          cd = 0
+          do c = nO+1, nOrb-nR
+            do d = c+1, nOrb-nR
+              cd = cd + 1
+    
+              mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) &
+                      + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps &
+                          + rho(a,d,m) * rho(b,c,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi
+
+              W(ab,state) = W(ab,state) + mat_tmp * U(cd,state)
+            enddo
+          enddo
+
+          ij = nVV
+          do i = nC+1, nO
+            do j = i+1, nO
+              ij = ij + 1
+  
+              mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps &
+                          + rho(i,b,m) * rho(a,j,m) * Om(m) / eps
+              end do
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi
+  
+              W(ab,state) = W(ab,state) - mat_tmp * U(ij,state)
+            enddo
+          enddo
+
+        enddo ! state
+      enddo ! b
+    enddo ! a
+
+    ! ---
+  
+    ij = nVV
+    do i = nC+1, nO
+      do j = i+1, nO
+        ij = ij + 1
+  
+        do state = 1, n_states_diag
+  
+          W(ij,state) = 0.d0
+  
+          ab = 0
+          do a = nO+1, nOrb-nR
+            do b = a+1, nOrb-nR
+              ab = ab + 1
+    
+              mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,a,m) * rho(j,b,m) * Om(m) / eps &
+                          + rho(i,b,m) * rho(a,j,m) * Om(m) / eps
+              end do
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi
+
+              W(ij,state) = W(ij,state) + mat_tmp * U(ab,state)
+            enddo
+          enddo
+
+          kl = nVV
+          do k = nC+1, nO
+            do l = k+1, nO
+              kl = kl + 1
+    
+              mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) &
+                      + lambda * (ERI(i,j,k,l) - ERI(i,j,l,k))
+
+              chi = 0.d0
+              do m = 1, nS
+                eps = Om(m)**2 + eta**2
+                chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps &
+                          + rho(i,l,m) * rho(j,k,m) * Om(m) / eps
+              enddo
+              mat_tmp = mat_tmp + 4.d0 * lambda * chi
+
+              W(ij,state) = W(ij,state) - mat_tmp * U(kl,state)
+            enddo
+          enddo
+
+        enddo ! state
+      enddo ! j
+    enddo ! i
+
+  else
+
+    print*, ' Error in ppLR_GW_HR_calc'
+    print*, ' ispin is not supported'
+    print*, ' ispin = ', ispin
+    stop
+
+  endif
+
+  return
+end
+
+! ---
+
+subroutine ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e, eF, ERI, eta, rho, Om, H_diag)
+
+  implicit none
+
+  integer,          intent(in)  :: ispin
+  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR, nS
+  double precision, intent(in)  :: lambda, eF, eta
+  double precision, intent(in)  :: e(nOrb)
+  double precision, intent(in)  :: ERI(nOrb,nOrb,nOrb,nOrb)
+  double precision, intent(in)  :: rho(nOrb,nOrb,nS), Om(nS)
+  double precision, intent(out) :: H_diag(nOO+nVV)
+
+  integer                       :: i, j, ij, k, l, kl
+  integer                       :: a, b, c, d, ab, cd
+  integer                       :: m
+  double precision              :: chi, eps
+
+  double precision, external    :: Kronecker_delta
+
+
+  if(ispin .eq. 1) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = a, nOrb-nR
+        ab = ab + 1
+        cd = 0
+        do c = nO+1, nOrb-nR
+          do d = c, nOrb-nR
+            cd = cd + 1
+            if(a .ne. c) cycle
+            if(b .ne. d) cycle
+            H_diag(ab) = e(a) + e(b) - eF &
+                       + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                                       * (1.d0 + Kronecker_delta(c, d)))
+            chi = 0.d0
+            do m = 1, nS
+              eps = Om(m)**2 + eta**2
+              chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps &
+                        - rho(a,d,m) * rho(b,c,m) * Om(m) / eps
+            end do
+            H_diag(ab) = H_diag(ab) + 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
+                                                                 * (1.d0 + Kronecker_delta(c, d)))
+          enddo
+        enddo
+      enddo ! b
+    enddo ! a
+
+    ij = nVV
+    do i = nC+1, nO
+      do j = i, nO
+        ij = ij + 1
+        kl = 0
+        do k = nC+1, nO
+          do l = k, nO
+            kl = kl + 1
+            if(i .ne. k) cycle
+            if(j .ne. l) cycle
+            H_diag(ij) = e(i) + e(j) - eF &
+                       - lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt( (1.d0 + Kronecker_delta(i, j)) &
+                                                                       * (1.d0 + Kronecker_delta(k, l)))
+            chi = 0.d0
+            do m = 1, nS
+              eps = Om(m)**2 + eta**2
+              chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps &
+                        - rho(i,l,m) * rho(j,k,m) * Om(m) / eps
+            enddo
+            H_diag(ij) = H_diag(ij) - 4.d0 * lambda * chi / dsqrt( (1.d0 + Kronecker_delta(i, j)) &
+                                                                 * (1.d0 + Kronecker_delta(k, l)))
+          enddo
+        enddo
+      enddo ! j
+    enddo ! i
+
+  elseif(ispin .eq. 2) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = a+1, nOrb-nR
+        ab = ab + 1
+        cd = 0
+        do c = nO+1, nOrb-nR
+          do d = c+1, nOrb-nR
+            cd = cd + 1
+            if(a .ne. c) cycle
+            if(b .ne. d) cycle
+            H_diag(ab) = e(a) + e(b) - eF + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c))
+            chi = 0.d0
+            do m = 1, nS
+              eps = Om(m)**2 + eta**2
+              chi = chi - rho(a,c,m) * rho(b,d,m) * Om(m) / eps &
+                        + rho(a,d,m) * rho(b,c,m) * Om(m) / eps
+            enddo
+            H_diag(ab) = H_diag(ab) + 4.d0 * lambda * chi
+          enddo
+        enddo
+      enddo ! b
+    enddo ! a
+  
+    ij = nVV
+    do i = nC+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
+            if(i .ne. k) cycle
+            if(j .ne. l) cycle
+            H_diag(ij) = e(i) + e(j) - eF - lambda * (ERI(i,j,k,l) - ERI(i,j,l,k))
+            chi = 0.d0
+            do m = 1, nS
+              eps = Om(m)**2 + eta**2
+              chi = chi - rho(i,k,m) * rho(j,l,m) * Om(m) / eps &
+                        + rho(i,l,m) * rho(j,k,m) * Om(m) / eps
+            end do
+            H_diag(ij) = H_diag(ij) - 4.d0 * lambda * chi
+          enddo
+        enddo
+      enddo ! j
+    enddo ! i
+
+  else
+
+    print*, ' Error in ppLR_GW_H_diag'
+    print*, ' ispin is not supported'
+    print*, ' ispin = ', ispin
+    stop
+
+  endif
+
+  return
+end
+
+! ---
+

src/LR/ppLR_RPA_davidson.f90

@@ -76,15 +76,15 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_
   
           W(ij,state) = 0.d0
   
-          cd = 0
+          ab = 0
-          do c = nO+1, nOrb-nR
+          do a = nO+1, nOrb-nR
-            do d = c, nOrb-nR
+            do b = a, nOrb-nR
-              cd = cd + 1
+              ab = ab + 1
     
-              mat_tmp = lambda * (ERI(c,d,i,j) + ERI(c,d,j,i)) / dsqrt( (1.d0 + Kronecker_delta(c, d)) &
+              mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) &
                                                                       * (1.d0 + Kronecker_delta(i, j)))
 
-              W(ij,state) = W(ij,state) + mat_tmp * U(cd,state)
+              W(ij,state) = W(ij,state) + mat_tmp * U(ab,state)
             enddo
           enddo
   
@@ -154,14 +154,14 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_
   
           W(ij,state) = 0.d0
   
-          cd = 0
+          ab = 0
-          do c = nO+1, nOrb-nR
+          do a = nO+1, nOrb-nR
-            do d = c+1, nOrb-nR
+            do b = a+1, nOrb-nR
-              cd = cd + 1
+              ab = ab + 1
     
-              mat_tmp = lambda * (ERI(c,d,i,j) - ERI(c,d,j,i))
+              mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i))
 
-              W(ij,state) = W(ij,state) + mat_tmp * U(cd,state)
+              W(ij,state) = W(ij,state) + mat_tmp * U(ab,state)
             enddo
           enddo
 
@@ -230,14 +230,14 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_
   
           W(ij,state) = 0.d0
   
-          cd = 0
+          ab = 0
-          do c = nO+1, nOrb-nR
+          do a = nO+1, nOrb-nR
-            do d = nO+1, nOrb-nR
+            do b = nO+1, nOrb-nR
-              cd = cd + 1
+              ab = ab + 1
     
-              mat_tmp = lambda * ERI(c,d,i,j)
+              mat_tmp = lambda * ERI(a,b,i,j)
 
-              W(ij,state) = W(ij,state) + mat_tmp * U(cd,state)
+              W(ij,state) = W(ij,state) + mat_tmp * U(ab,state)
             enddo
           enddo
   
@@ -259,6 +259,7 @@ subroutine ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_
 
   else
 
+    print*, ' Error in ppLR_RPA_HR_calc'
     print*, ' ispin is not supported'
     print*, ' ispin = ', ispin
     stop
@@ -445,6 +446,7 @@ subroutine ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI
 
   else
 
+    print*, ' Error in ppLR_RPA_H_diag'
     print*, ' ispin is not supported'
     print*, ' ispin = ', ispin
     stop

src/LR/ppLR_davidson.f90

@@ -1,7 +1,10 @@
 
 ! ---
 
-subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, Om, R, n_states, n_states_diag, kernel)
+subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, &
+                         supp_data_int, supp_data_int_size,                          &
+                         supp_data_dbl, supp_data_dbl_size,                          &
+                         Om, R, n_states, n_states_diag, kernel)
 
   ! 
   ! Extract the low n_states 
@@ -22,17 +25,22 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   integer,          intent(in)  :: nC, nO, nR, nOrb, nOO, nVV
   integer,          intent(in)  :: n_states      ! nb of physical states
   integer,          intent(in)  :: n_states_diag ! nb of states used to get n_states
+  integer,          intent(in)  :: supp_data_int_size
+  integer,          intent(in)  :: supp_data_dbl_size
   character(len=*), intent(in)  :: kernel
   double precision, intent(in)  :: lambda, eF
   double precision, intent(in)  :: e(nOrb)
   double precision, intent(in)  :: ERI(nOrb,nOrb,nOrb,nOrb)
+  integer,          intent(in)  :: supp_data_int(supp_data_int_size)
+  double precision, intent(in)  :: supp_data_dbl(supp_data_dbl_size)
   double precision, intent(out) :: Om(n_states)
   double precision, intent(out) :: R(nOO+nVV,n_states_diag)
   
   integer                       :: N, M
   integer                       :: iter, itermax, itertot
   integer                       :: shift1, shift2
-  integer                       :: i, j, ij, k, l, kl, a, b, ab, c, d, cd
+  integer                       :: i, j, k, l, ab
+  integer                       :: p, q, mm, i_data, nS
   integer                       :: i_omax(n_states)
   logical                       :: converged
   character*(16384)             :: write_buffer
@@ -40,6 +48,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   double precision              :: lambda_tmp
   double precision              :: to_print(2,n_states)
   double precision              :: mem
+  double precision              :: eta
+  character(len=len(kernel))    :: kernel_name
   double precision, allocatable :: H_diag(:)
   double precision, allocatable :: W(:,:)
   double precision, allocatable :: U(:,:)
@@ -47,6 +57,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   double precision, allocatable :: residual_norm(:)
   double precision, allocatable :: overlap(:)
   double precision, allocatable :: S_check(:,:)
+  double precision, allocatable :: rho_tmp(:,:,:), Om_tmp(:)
                                 
   double precision, external    :: u_dot_u
 
@@ -56,6 +67,8 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   itermax = 8
   M = n_states_diag * itermax
 
+  call lower_case(trim(kernel), kernel_name)
+
   if(M .ge. N) then
     print*, 'N = ', N
     print*, 'M = ', M
@@ -71,7 +84,7 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   write(*,'(A40, I12)') 'Number of states = ', n_states
   write(*,'(A40, I12)') 'Number of states in diagonalization = ', n_states_diag
   write(*,'(A40, I12)') 'Number of basis functions = ', N
-
+  write(*,'(A40, A12)') 'Kernel: ', kernel_name
 
 
 
@@ -82,14 +95,55 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   allocate(overlap(n_states_diag))
   allocate(residual_norm(n_states_diag))
 
-  mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) / 1d6
+  mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) &
-  write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem
+      + 8.d0 * dble(supp_data_dbl_size) + 4.d0 * dble(supp_data_int_size)
 
-  mem = 8.d0 * dble(N + N*M + N*M + M*M + M*M + M + n_states_diag + n_states_diag) / 1d6
+  write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem / (1024.d0*1024.d0)
-  write(*,'(A40, F12.4)') 'tmp mem (MB) = ', mem
+
+  mem = 8.d0 * dble(N + N*M + N*M + M*M + M*M + M + n_states_diag + n_states_diag)
+
+  write(*,'(A40, F12.4)') 'tmp mem (MB) = ', mem / (1024.d0*1024.d0)
 
 
-  call ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag, kernel)
+  if(kernel_name .eq. "rpa") then
+
+    call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, &
+                         ERI(1,1,1,1), H_diag(1))
+
+  elseif(kernel_name .eq. "gw") then
+
+    nS = supp_data_int(1)
+
+    allocate(rho_tmp(nOrb,nOrb,nS))
+    allocate(Om_tmp(nS))
+
+    eta = supp_data_dbl(1)
+    i_data = 1
+    do mm = 1, nS
+      do q = 1, nOrb
+        do p = 1, nOrb
+          i_data = i_data + 1
+          rho_tmp(p,q,mm) = supp_data_dbl(i_data)
+        enddo
+      enddo
+    enddo
+    do mm = 1, nS
+      i_data = i_data + 1
+      Om_tmp(mm) = supp_data_dbl(i_data)
+    enddo
+
+    call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, &
+                        ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1), H_diag(1))
+
+  !! TODO
+  !elseif(kernel_name .eq. "gf2") then
+
+  else
+
+    print*, ' kernel not supported', kernel
+    stop
+
+  endif
 
   !print*, "H_diag:"
   !do ab = 1, N
@@ -187,8 +241,22 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
         !enddo
         !deallocate(S_check)
 
-        call ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, &
+        if(kernel_name .eq. "rpa") then
-                          ERI(1,1,1,1), U(1,shift1+1), W(1,shift1+1), kernel)
+      
+          call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e(1), eF, n_states_diag, &
+                                ERI(1,1,1,1),                                                       &
+                                U(1,shift1+1), W(1,shift1+1))
+      
+        elseif(kernel_name .eq. "gw") then
+      
+          call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, nS, lambda, e(1), eF, n_states_diag, &
+                               ERI(1,1,1,1), eta, rho_tmp(1,1,1), Om_tmp(1),                           &
+                               U(1,shift1+1), W(1,shift1+1))
+      
+        !! TODO
+        !elseif(kernel_name .eq. "gf2") then
+      
+        endif
 
       else
 
@@ -353,92 +421,9 @@ subroutine ppLR_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF,
   deallocate(overlap)
   deallocate(residual_norm)
 
-  return
+  if(kernel_name .eq. "gw") then
-end
+    deallocate(rho_tmp)
-
+    deallocate(Om_tmp)
-! ---
-
-subroutine ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W, kernel)
-
-  implicit none
-
-  integer,          intent(in)  :: ispin
-  integer,          intent(in)  :: n_states_diag
-  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR
-  character(len=*), intent(in)  :: kernel
-  double precision, intent(in)  :: lambda, eF
-  double precision, intent(in)  :: e(nOrb)
-  double precision, intent(in)  :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision, intent(in)  :: U(nOO+nVV,n_states_diag)
-  double precision, intent(out) :: W(nOO+nVV,n_states_diag)
-
-  character(len=len(kernel))    :: kernel_name
-
-  call lower_case(trim(kernel), kernel_name)
-
-  if(kernel_name .eq. "rpa") then
-
-    call ppLR_RPA_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W)
-
-  !! TODO
-  !elseif(kernel_name .eq. "gw") then
-
-  !  call ppLR_GW_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W)
-
-  !! TODO
-  !elseif(kernel_name .eq. "gf2") then
-
-  !  call ppLR_GF2_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W)
-
-  else 
-
-    print*, ' Error in routine ppLR_HR_calc'
-    print*, ' kernel not supported', kernel
-    stop
-
-  endif
-
-  return
-end
-
-! ---
-
-subroutine ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag, kernel)
-
-  implicit none
-
-  integer,          intent(in)  :: ispin
-  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR
-  character(len=*), intent(in)  :: kernel
-  double precision, intent(in)  :: lambda, eF
-  double precision, intent(in)  :: e(nOrb)
-  double precision, intent(in)  :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision, intent(out) :: H_diag(nOO+nVV)
-
-  character(len=len(kernel))    :: kernel_name
-
-  call lower_case(trim(kernel), kernel_name)
-
-  if(kernel_name .eq. "rpa") then
-
-    call ppLR_RPA_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag)
-
-  !! TODO
-  !elseif(kernel_name .eq. "gw") then
-
-  !  call ppLR_GW_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag)
-
-  !! TODO
-  !elseif(kernel_name .eq. "gf2") then
-
-  !  call ppLR_GF2_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag)
-
-  else
-
-    print*, ' Error in routine ppLR_H_diag'
-    print*, ' kernel not supported', kernel
-    stop
-
   endif
 
   return

src/utils/non_sym_diag.f90

@@ -465,7 +465,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
           accu_nd = accu_nd + dabs(S(j,k))
         enddo
       enddo
-      if(accu_nd .gt. 1d-12) then
+      if(accu_nd .gt. 1d-7) then
         print*, ' accu_nd =', accu_nd
         print*, ' your strategy for degenerates orbitals failed !'
         print*, m, 'deg on', i