added Davidson for ppLR_RG0T0pp

src/GT/RG0T0pp.f90

@@ -1,5 +1,5 @@
 subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & 
-                   linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+                   linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform one-shot calculation with a T-matrix self-energy (G0T0)
 
@@ -25,7 +25,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   double precision,intent(in)   :: eta
   logical,intent(in)            :: regularize
 
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
@@ -33,9 +33,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   integer,intent(in)            :: nS
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: ERHF
-  double precision,intent(in)   :: eHF(nBas)
+  double precision,intent(in)   :: eHF(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
 
@@ -45,6 +45,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   integer                       :: iblock
   integer                       :: nOOs,nOOt
   integer                       :: nVVs,nVVt
+  integer                       :: n_states, n_states_diag
   double precision              :: EcRPA(nspin)
   double precision              :: EcBSE(nspin)
   double precision              :: EcGM
@@ -63,6 +64,7 @@ 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(:)
+  double precision, allocatable :: Om(:), R(:,:)
 
 
 ! Output variables
@@ -92,8 +94,8 @@ 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
   nVVs = nV*nV
@@ -105,30 +107,39 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs),    & 
            Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs),    & 
-           rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & 
+           rho1s(nOrb,nOrb,nVVs),rho2s(nOrb,nOrb,nOOs), & 
            Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt),    & 
            Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt),    & 
-           rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & 
+           rho1t(nOrb,nOrb,nVVt),rho2t(nOrb,nOrb,nOOt), & 
-           Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas))
+           Sig(nOrb),Z(nOrb),eGT(nOrb),eGTlin(nOrb))
 
 !----------------------------------------------
 ! alpha-beta block
 !----------------------------------------------
 
   ispin  = 1
-! iblock = 1
+  !iblock = 1
   iblock = 3
 
 ! Compute linear response
 
   allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
-  call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
+  call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
-  call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
+  call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
-  if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
 
   call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
 
+  !n_states = nOOs + 5
+  !n_states_diag = n_states + 4
+  !allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag))
+  !call ppLR_RG0T0pp_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag)
+  !print*, 'LAPACK:'
+  !print*, Om2s
+  !print*, Om1s
+  !stop
+
   deallocate(Bpp,Cpp,Dpp)
 
   if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s)
@@ -146,12 +157,21 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
-  call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
+  call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
-  call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
+  call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
-  if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
 
   call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
 
+  n_states = nOOt + 5
+  n_states_diag = n_states + 4
+  allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag))
+  call ppLR_RG0T0pp_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag)
+  print*, 'LAPACK:'
+  print*, Om2t
+  print*, Om1t
+  stop
+
   deallocate(Bpp,Cpp,Dpp)
 
   if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
@@ -164,23 +184,23 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 ! iblock = 1
   iblock = 3
 
-  call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+  call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
 ! iblock = 2
   iblock = 4
 
-  call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+  call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
 
   if(regularize) then 
-    call GTpp_regularization(nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Om1s,rho1s,Om2s,rho2s)
+    call GTpp_regularization(nOrb,nC,nO,nV,nR,nOOs,nVVs,eHF,Om1s,rho1s,Om2s,rho2s)
-    call GTpp_regularization(nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Om1t,rho1t,Om2t,rho2t)
+    call GTpp_regularization(nOrb,nC,nO,nV,nR,nOOt,nVVt,eHF,Om1t,rho1t,Om2t,rho2t)
   end if
 
-  call RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
+  call RGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
                              Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
 
 !----------------------------------------------
@@ -201,12 +221,12 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
     write(*,*) ' *** Quasiparticle energies obtained by root search *** '
     write(*,*)
      
-   call RGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
+   call RGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & 
                       Om1t,rho1t,Om2t,rho2t,eGTlin,eHF,eGT,Z)
 
   end if
 
-! call RGTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
+! call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
 !                            Om1t,rho1t,Om2t,rho2t)
 
 !----------------------------------------------
@@ -221,9 +241,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
-  call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
+  call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
-  call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
+  call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
-  if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
 
   call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
 
@@ -235,9 +255,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
-  call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
+  call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
-  call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
+  call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
-  if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
 
   call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
 
@@ -246,13 +266,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
   EcRPA(1) = EcRPA(1) - EcRPA(2)
   EcRPA(2) = 3d0*EcRPA(2)
 
-  call print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
+  call print_RG0T0pp(nOrb,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)
 
 ! Perform BSE calculation
 
   if(dophBSE) then
 
-    call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,   & 
+    call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt,   & 
                     Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
                     ERI,dipole_int,eHF,eGT,EcBSE)
 
@@ -288,7 +308,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
       end if
 
-      call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
+      call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
                         nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,     & 
                         Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE)
 
@@ -314,7 +334,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
   if(doppBSE) then
 
-    call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,      &
+    call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, &
                     Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
                     ERI,dipole_int,eHF,eGT,EcBSE)
 

src/GW/RGW_ppBSE.f90

@@ -125,6 +125,9 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS
 
     call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
 
+    ! TODO
+    !call ppLR_RGW_ppBSE_davidson(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+
     call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
 
     !----------------------------------------------------!

src/LR/ppLR_B.f90

@@ -1,4 +1,4 @@
-subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
+subroutine ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
 
 ! Compute the B matrix of the pp channel
 
@@ -8,7 +8,7 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
 ! 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
@@ -16,7 +16,7 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
   integer,intent(in)            :: nOO
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb) 
   
 ! Local variables
 
@@ -32,8 +32,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
   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
@@ -54,8 +54,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
   if(ispin == 2 .or. ispin == 4) 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
@@ -76,8 +76,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
   if(ispin == 3) then
 
     ab = 0
-    do a=nO+1,nBas-nR
+    do a=nO+1,nOrb-nR
-     do b=nO+1,nBas-nR
+     do b=nO+1,nOrb-nR
         ab = ab + 1
         ij = 0
         do i=nC+1,nO

src/LR/ppLR_C.f90

@@ -1,4 +1,4 @@
-subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
+subroutine ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 
 ! Compute the C matrix of the pp channel
 
@@ -8,14 +8,14 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 ! 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
   integer,intent(in)            :: nR
   integer,intent(in)            :: nVV
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: e(nOrb),ERI(nOrb,nOrb,nOrb,nOrb) 
   
 ! Local variables
 
@@ -39,15 +39,15 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 
   if(ispin == 1) then
 
-    a0 = nBas - nR - nO
+    a0 = nOrb - nR - nO
 
     !$OMP PARALLEL DEFAULT(NONE)                                                   &
     !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, e_ab, tmp_ab, delta_ac, tmp_cd) &
-    !$OMP          SHARED(nO, nBas, nR, a0, eF, lambda, e, ERI, Cpp)
+    !$OMP          SHARED(nO, nOrb, nR, a0, eF, lambda, e, ERI, Cpp)
     !$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
 
         e_ab = e(a) + e(b) - eF
@@ -58,14 +58,14 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
         endif
 
         cd = 0
-        do c = nO+1, nBas-nR
+        do c = nO+1, nOrb-nR
 
           delta_ac = 0.d0
           if(a .eq. c) then
             delta_ac = 1.d0
           endif
 
-         do d = c, nBas-nR
+         do d = c, nOrb-nR
             cd = cd + 1
 
             tmp_cd = tmp_ab
@@ -88,12 +88,12 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 
 
 !    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
 ! 
 !            Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
@@ -110,16 +110,16 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
 
   if(ispin == 2 .or. ispin == 4) then
     !$OMP PARALLEL &
-    !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nBas,nR) &
+    !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nOrb,nR) &
     !$OMP PRIVATE(c,d,a,b,ab,cd) &
     !$OMP DEFAULT(NONE)
     !$OMP DO
-    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 = (c-(nO+1))*(nBas-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c
+        cd = (c-(nO+1))*(nOrb-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c
-          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 = (a-(nO+1))*(nBas-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a
+              ab = (a-(nO+1))*(nOrb-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a
  
               Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
                    + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c))
@@ -138,12 +138,12 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
   if(ispin == 3) then
 
     ab = 0
-    do a=nO+1,nBas-nR
+    do a=nO+1,nOrb-nR
-     do b=nO+1,nBas-nR
+     do b=nO+1,nOrb-nR
         ab = ab + 1
         cd = 0
-        do c=nO+1,nBas-nR
+        do c=nO+1,nOrb-nR
-         do d=nO+1,nBas-nR
+         do d=nO+1,nOrb-nR
             cd = cd + 1
  
             Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 

src/LR/ppLR_D.f90

@@ -1,4 +1,4 @@
-subroutine ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
+subroutine ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
 
 ! Compute the D matrix of the pp channel
 
@@ -8,14 +8,14 @@ subroutine ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
 ! 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
   integer,intent(in)            :: nR
   integer,intent(in)            :: nOO
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: e(nOrb),ERI(nOrb,nOrb,nOrb,nOrb) 
   
 ! Local variables
 

src/LR/ppLR_RG0T0pp_davidson.f90

@@ -0,0 +1,982 @@
+
+! ---
+
+subroutine ppLR_RG0T0pp_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, Om, R, n_states, n_states_diag)
+
+  ! 
+  ! Extract the low n_states 
+  ! Om(i) (eigenvalues) and 
+  ! R(:,i) (right-eigenvectors) 
+  ! of the pp-RPA matrix
+  !
+  !   (+C    +B) 
+  !   (        )  
+  !   (-B.T  -D)
+  !
+
+  implicit none
+
+  logical,          intent(in)  :: TDA
+  integer,          intent(in)  :: ispin
+  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
+  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) :: 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_omax(n_states)
+  logical                       :: converged
+  character*(16384)             :: write_buffer
+  double precision              :: r1, r2, dtwo_pi
+  double precision              :: lambda_tmp
+  double precision              :: to_print(2,n_states)
+  double precision              :: mem
+  double precision, allocatable :: H_diag(:)
+  double precision, allocatable :: W(:,:)
+  double precision, allocatable :: U(:,:)
+  double precision, allocatable :: h(:,:), h_vec(:,:), h_val(:)
+  double precision, allocatable :: residual_norm(:)
+  double precision, allocatable :: overlap(:)
+  double precision, allocatable :: S_check(:,:)
+
+  double precision, external    :: u_dot_u
+
+  dtwo_pi = 6.283185307179586d0
+
+  N = nOO + nVV
+  itermax = 8
+  M = n_states_diag * itermax
+
+  if(M .ge. N) then
+    print*, 'N = ', N
+    print*, 'M = ', M
+    print*, ' use Lapack or decrease n_states and/or itermax '
+    stop
+  endif
+
+  write(6,'(A)') ''
+  write(6,'(A)') 'Davidson Diagonalization'
+  write(6,'(A)') '------------------------'
+  write(6,'(A)') ''
+
+  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
+
+
+
+
+  allocate(H_diag(N))
+  allocate(U(N,M))
+  allocate(W(N,M))
+  allocate(h(M,M), h_vec(M,M), h_val(M))
+  allocate(overlap(n_states_diag))
+  allocate(residual_norm(n_states_diag))
+
+  mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) / 1d6
+  write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem
+
+  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)') 'tmp mem (MB) = ', mem
+
+
+  call ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag)
+
+  !print*, "H_diag:"
+  !do ab = 1, N
+  !  print*, ab, H_diag(ab)
+  !enddo
+
+  ! initialize guess
+  R = 0.d0
+  do k = 1, n_states
+    R(k,k) = 1.d0
+  enddo
+  do k = n_states+1, n_states_diag
+    do i = 1, N
+      call random_number(r1)
+      call random_number(r2)
+      r1 = dsqrt(-2.d0*dlog(r1))
+      r2 = dtwo_pi*r2
+      R(i,k) = r1*dcos(r2)
+    enddo
+    R(k,k) = R(k,k) + 10.d0
+  enddo
+
+  do k = 1, n_states_diag
+    call normalize(R(1,k), N)
+  enddo
+
+  !print*, 'guess'
+  !do k = 1, N
+  !  write(*,'(100(F15.7,2X))') (R(k,i), i = 1, n_states_diag)
+  !enddo
+
+  ! working vectors
+  do k = 1, n_states_diag
+    do i = 1, N
+      U(i,k) = R(i,k)
+    enddo
+  enddo
+
+  !print*, 'working vectors'
+  !do k = 1, N
+  !  write(*,'(100(F15.7,2X))') (U(k,i), i = 1, n_states_diag)
+  !enddo
+
+  write(6,'(A)') ''
+  write_buffer = '====='
+  do i = 1, n_states
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*n_states)
+  write_buffer = 'Iter'
+  do i = 1, n_states
+    write_buffer = trim(write_buffer)//'       Energy         Residual '
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*n_states)
+  write_buffer = '====='
+  do i = 1, n_states
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*n_states)
+
+
+  converged = .False.
+  itertot = 0
+
+  do while (.not.converged)
+
+    itertot = itertot + 1
+    if(itertot == itermax) then
+      print*, 'exit before convergence !'
+      print*, 'itertot == itermax', itertot
+      exit
+    endif
+
+    do iter = 1, itermax-1
+
+      shift1 = n_states_diag * (iter - 1)
+      shift2 = shift1 + n_states_diag
+      !print*, iter, shift1, shift2
+
+      if((iter > 1) .or. (itertot == 1)) then
+
+        call ortho_qr(U(1,1), size(U, 1), N, shift2)
+        !call ortho_qr(U(1,1), size(U, 1), N, shift2)
+
+        !print*, 'working vectors after qr'
+        !do k = 1, N
+        !  write(*,'(100(F15.7,2X))') (U(k,i), i = 1, n_states_diag)
+        !enddo
+        !allocate(S_check(shift2,shift2))
+        !call dgemm("T", "N", shift2, shift2, N, 1.d0, U(1,1), size(U, 1), U(1,1), size(U, 1), 0.d0, S_check(1,1), size(S_check, 1))
+        !do k = 1, shift2
+        !  write(*,'(100(F15.7,2X))') (S_check(k,i), i = 1, shift2)
+        !enddo
+        !deallocate(S_check)
+
+        call ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI(1,1,1,1), U(1,shift1+1), W(1,shift1+1))
+
+      else
+
+        ! computed below
+        continue
+      endif
+
+      ! h = U.T H U
+      call dgemm('T', 'N', shift2, shift2, N, 1.d0,      &
+                 U(1,1), size(U, 1), W(1,1), size(W, 1), &
+                 0.d0, h(1,1), size(h, 1))
+
+      ! h h_vec = h_val h_vec
+      call diag_nonsym_right(shift2, h(1,1), size(h, 1), h_vec(1,1), size(h_vec, 1), h_val(1), size(h_val, 1))
+      !print*, 'h_val', h_val(1:shift2)
+
+      ! U1 = U0 h_vec
+      call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0,       &
+                 U(1,1), size(U, 1), h_vec(1,1), size(h_vec, 1), &
+                 0.d0, U(1,shift2+1), size(U, 1))
+
+      do k = 1, n_states_diag
+        call normalize(U(1,shift2+k), N)
+      enddo
+
+      do l = 1, n_states
+        !do k = 1, n_states_diag
+        !  overlap(k) = 0.d0
+        !  do i = 1, N
+        !    overlap(k) = overlap(k) + U(i,shift2+k) * R(i,l)
+        !  enddo
+        !  overlap(k) = dabs(overlap(k))
+        !  !print *, ' overlap =', k, overlap(k)
+        !enddo
+        !lambda_tmp = 0.d0
+        !do k = 1, n_states_diag
+        !  if(overlap(k) .gt. lambda_tmp) then
+        !    i_omax(l) = k
+        !    lambda_tmp = overlap(k)
+        !  endif
+        !enddo
+        !if(lambda_tmp .lt. 0.7d0) then
+        !  print *, ' small overlap ...', l, i_omax(l)
+        !  print *, ' max overlap =', lambda_tmp
+        !  !stop
+        !endif
+        !if(i_omax(l) .ne. l) then
+        !  print *, ' !!! WARNING !!!'
+        !  print *, ' index of state', l, i_omax(l)
+        !endif
+      enddo
+
+      ! W1 = W0 h_vec
+      call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0,       &
+                 W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), &
+                 0.d0, W(1,shift2+1), size(W, 1))
+
+      ! check if W1 = U1 h_val
+      !$OMP PARALLEL DEFAULT(NONE) &
+      !$OMP          PRIVATE(i, k) &
+      !$OMP          SHARED(n_states, n_states_diag, N, shift2, U, h_val, W, H_diag, residual_norm, to_print)
+      !$OMP DO 
+      do k = 1, n_states_diag
+        do i = 1, N
+          U(i,shift2+k) = (h_val(k) * U(i,shift2+k) - W(i,shift2+k)) / max(H_diag(i) - h_val(k), 1.d-2)
+        enddo
+        if(k <= n_states) then
+          residual_norm(k) = u_dot_u(U(1,shift2+k), N)
+          to_print(1,k) = h_val(k)
+          to_print(2,k) = residual_norm(k)
+        endif
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+      !print*, " to_print", to_print
+
+      if((itertot > 1) .and. (iter == 1)) then
+        continue
+      else
+        write(*,'(1X, I3, 1X, 100(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states)
+        !write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, F16.10, 1X, F16.10))') iter-1, to_print(1:2,1:n_states)
+      endif
+
+      !print*, 'iter = ', iter
+      if(iter > 1) then
+        converged = dabs(maxval(residual_norm(1:n_states))) < 1d-15
+      endif
+
+      do k = 1, n_states
+        if(residual_norm(k) > 1.d8) then
+          print *, 'Davidson failed'
+          stop -1
+        endif
+      enddo
+
+      if(converged) exit
+
+    enddo ! loop over iter
+
+
+    ! Re-contract U and update W
+    ! --------------------------------
+
+    call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0,       &
+               W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), &
+               0.d0, R, size(R, 1))
+    do k = 1, n_states_diag
+      do i = 1, N
+        W(i,k) = R(i,k)
+      enddo
+    enddo
+
+    call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0,       &
+               U(1,1), size(U, 1), h_vec(1,1), size(h_vec, 1), &
+               0.d0, R(1,1), size(R, 1))
+
+    do k = 1, n_states_diag
+      do i = 1, N
+        U(i,k) = R(i,k)
+      enddo
+    enddo
+
+    call ortho_qr(U(1,1), size(U, 1), N, n_states_diag)
+    !call ortho_qr(U(1,1), size(U, 1), N, n_states_diag)
+
+    do j = 1, n_states_diag
+      k = 1
+      do while((k < N) .and. (U(k,j) == 0.d0))
+        k = k+1
+      enddo
+      if(U(k,j) * R(k,j) < 0.d0) then
+        do i = 1, N
+          W(i,j) = -W(i,j)
+        enddo
+      endif
+    enddo
+
+  enddo ! loop over while
+
+  ! ---
+
+  write_buffer = '====='
+  do i = 1, n_states
+    write_buffer = trim(write_buffer)//' ================  ==========='
+  enddo
+  write(6,'(A)') trim(write_buffer)
+  write(6,'(A)') ''
+
+
+  print*, " Davidson eigenvalues"
+  do k = 1, n_states
+    Om(k) = h_val(k)
+    print*, k, Om(k)
+  enddo
+
+  deallocate(H_diag)
+  deallocate(U)
+  deallocate(W)
+  deallocate(h)
+  deallocate(h_vec)
+  deallocate(h_val)
+  deallocate(overlap)
+  deallocate(residual_norm)
+
+  return
+end
+
+! ---
+
+subroutine ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W)
+
+  implicit none
+
+  integer,          intent(in)  :: ispin
+  integer,          intent(in)  :: n_states_diag
+  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR
+  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)
+
+  integer                       :: i, j, ij, k, l, kl
+  integer                       :: a, b, c, d, ab, cd
+  integer                       :: state
+  double precision              :: mat_tmp
+  double precision              :: diff_loc, diff_tot
+  double precision, allocatable :: M_ref(:,:), W_ref(:,:)
+  double precision, allocatable :: Cpp_ref(:,:), Dpp_ref(:,:), Bpp_ref(:,:)
+
+  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)))
+
+              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)))
+  
+              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
+  
+          cd = 0
+          do c = nO+1, nOrb-nR
+            do d = c, nOrb-nR
+              cd = cd + 1
+    
+              mat_tmp = lambda * (ERI(c,d,i,j) + ERI(c,d,j,i)) / dsqrt( (1.d0 + Kronecker_delta(c, d)) &
+                                                                      * (1.d0 + Kronecker_delta(i, j)))
+
+              W(ij,state) = W(ij,state) + mat_tmp * U(cd,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)))
+
+              W(ij,state) = W(ij,state) - mat_tmp * U(kl,state)
+            enddo
+          enddo
+  
+        enddo ! state
+      enddo ! j
+    enddo ! i
+
+  elseif((ispin .eq. 2) .or. (ispin .eq. 4)) 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))
+
+              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))
+  
+              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
+  
+          cd = 0
+          do c = nO+1, nOrb-nR
+            do d = c+1, nOrb-nR
+              cd = cd + 1
+    
+              mat_tmp = lambda * (ERI(c,d,i,j) - ERI(c,d,j,i))
+
+              W(ij,state) = W(ij,state) + mat_tmp * U(cd,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))
+
+              W(ij,state) = W(ij,state) - mat_tmp * U(kl,state)
+            enddo
+          enddo
+  
+        enddo ! state
+      enddo ! j
+    enddo ! i
+
+  elseif(ispin .eq. 3) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = nO+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 = nO+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)
+
+              W(ab,state) = W(ab,state) + mat_tmp * U(cd,state)
+            enddo
+          enddo
+
+          ij = nVV
+          do i = nC+1, nO
+            do j = nC+1, nO
+              ij = ij + 1
+  
+              mat_tmp = lambda * ERI(a,b,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 = nC+1, nO
+        ij = ij + 1
+  
+        do state = 1, n_states_diag
+  
+          W(ij,state) = 0.d0
+  
+          cd = 0
+          do c = nO+1, nOrb-nR
+            do d = nO+1, nOrb-nR
+              cd = cd + 1
+    
+              mat_tmp = lambda * ERI(c,d,i,j)
+
+              W(ij,state) = W(ij,state) + mat_tmp * U(cd,state)
+            enddo
+          enddo
+  
+          kl = nVV
+          do k = nC+1, nO
+            do l = nC+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)
+
+              W(ij,state) = W(ij,state) - mat_tmp * U(kl,state)
+            enddo
+          enddo
+  
+        enddo ! state
+      enddo ! j
+    enddo ! i
+
+  else
+
+    print*, ' ispin is not supported'
+    print*, ' ispin = ', ispin
+    stop
+
+  endif
+
+
+!  print*, ' debug ppLR_HR_calc:'
+!  print*, ispin, nOO, nVV
+!  allocate(M_ref(nOO+nVV,nOO+nVV))
+!  allocate(Bpp_ref(nVV,nOO), Cpp_ref(nVV,nVV), Dpp_ref(nOO,nOO))
+!  allocate(W_ref(nOO+nVV,n_states_diag))
+!
+!  call ppLR_C(ispin, nOrb, nC, nO, nOrb-nO, nR, nVV, 1d0, e, ERI, Cpp_ref)
+!  call ppLR_D(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, 1d0, e, ERI, Dpp_ref)
+!  call ppLR_B(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, nVV, 1d0, ERI, Bpp_ref)
+!  M_ref = 0.d0
+!  M_ref(    1:nVV    ,    1:nVV)     = + Cpp_ref(1:nVV,1:nVV)
+!  M_ref(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp_ref(1:nOO,1:nOO)
+!  M_ref(    1:nVV    ,nVV+1:nOO+nVV) = - Bpp_ref(1:nVV,1:nOO)
+!  M_ref(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp_ref(1:nVV,1:nOO))
+!
+!  call dgemm('N', 'N', nOO+nVV, n_states_diag, nOO+nVV, 1.d0, &
+!             M_ref(1,1), size(M_ref, 1), U(1,1), size(U, 1),  &
+!             0.d0, W_ref(1,1), size(W_ref, 1))
+!
+!  diff_tot = 0.d0
+!  do state = 1, n_states_diag
+!    do ab = 1, nOO
+!      diff_loc = dabs(W(ab,state) - W_ref(ab,state))
+!      if(diff_loc .gt. 1d-12) then
+!        print*, ' important diff on:', ab, state
+!        print*, W(ab,state), W_ref(ab,state)
+!        stop
+!      endif
+!      diff_tot = diff_tot + diff_loc
+!    enddo
+!    do ij = nVV+1, nVV+nOO
+!      diff_loc = dabs(W(ij,state) - W_ref(ij,state))
+!      if(diff_loc .gt. 1d-12) then
+!        print*, ' important diff on:', ij, state
+!        print*, W(ij,state), W_ref(ij,state)
+!        stop
+!      endif
+!      diff_tot = diff_tot + diff_loc
+!    enddo
+!  enddo
+!  print*, 'diff_tot = ', diff_tot
+!
+!  deallocate(M_ref)
+!  deallocate(Bpp_ref, Cpp_ref, Dpp_ref)
+!  deallocate(W_ref)
+
+  return
+end
+
+! ---
+
+subroutine ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag)
+
+  implicit none
+
+  integer,          intent(in)  :: ispin
+  integer,          intent(in)  :: nOO, nVV, nOrb, nC, nO, nR
+  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)
+
+  integer                       :: i, j, ij, k, l, kl
+  integer                       :: a, b, c, d, ab, cd
+  double precision              :: diff_loc, diff_tot
+  double precision, allocatable :: M_ref(:,:)
+  double precision, allocatable :: Cpp_ref(:,:), Dpp_ref(:,:), Bpp_ref(:,:)
+
+  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)))
+          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)))
+          enddo
+        enddo
+      enddo ! j
+    enddo ! i
+
+  elseif((ispin .eq. 2) .or. (ispin .eq. 4)) 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))
+          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))
+          enddo
+        enddo
+      enddo ! j
+    enddo ! i
+
+  elseif(ispin .eq. 3) then
+
+    ab = 0
+    do a = nO+1, nOrb-nR
+      do b = nO+1, nOrb-nR
+        ab = ab + 1
+        cd = 0
+        do c = nO+1, nOrb-nR
+          do d = nO+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)
+          enddo
+        enddo
+      enddo ! b
+    enddo ! a
+
+    ij = nVV
+    do i = nC+1, nO
+      do j = nC+1, nO
+        ij = ij + 1 
+        kl = 0
+        do k = nC+1, nO
+          do l = nC+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)
+          enddo
+        enddo
+      enddo ! j
+    enddo ! i
+
+  else
+
+    print*, ' ispin is not supported'
+    print*, ' ispin = ', ispin
+    stop
+
+  endif
+
+
+!  print*, ' debug ppLR_H_diag:'
+!  print*, ispin, nOO, nVV
+!  allocate(M_ref(nOO+nVV,nOO+nVV))
+!  allocate(Bpp_ref(nVV,nOO), Cpp_ref(nVV,nVV), Dpp_ref(nOO,nOO))
+!
+!  call ppLR_C(ispin, nOrb, nC, nO, nOrb-nO, nR, nVV, 1d0, e, ERI, Cpp_ref)
+!  call ppLR_D(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, 1d0, e, ERI, Dpp_ref)
+!  call ppLR_B(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, nVV, 1d0, ERI, Bpp_ref)
+!  M_ref = 0.d0
+!  M_ref(    1:nVV    ,    1:nVV)     = + Cpp_ref(1:nVV,1:nVV)
+!  M_ref(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp_ref(1:nOO,1:nOO)
+!  M_ref(    1:nVV    ,nVV+1:nOO+nVV) = - Bpp_ref(1:nVV,1:nOO)
+!  M_ref(nVV+1:nOO+nVV,    1:nVV)     = + transpose(Bpp_ref(1:nVV,1:nOO))
+!
+!  diff_tot = 0.d0
+!  do ab = 1, nOO
+!    diff_loc = dabs(H_diag(ab) - M_ref(ab,ab))
+!    if(diff_loc .gt. 1d-12) then
+!      print*, ' important diff on:', ab
+!      print*, H_diag(ab), M_ref(ab,ab)
+!      stop
+!    endif
+!    diff_tot = diff_tot + diff_loc
+!  enddo
+!  do ij = nVV+1, nVV+nOO
+!    diff_loc = dabs(H_diag(ij) - M_ref(ij,ij))
+!    if(diff_loc .gt. 1d-12) then
+!      print*, ' important diff on:', ij
+!      print*, H_diag(ij), M_ref(ij,ij)
+!      stop
+!    endif
+!    diff_tot = diff_tot + diff_loc
+!  enddo
+!  print*, 'diff_tot = ', diff_tot
+!
+!  deallocate(M_ref)
+!  deallocate(Bpp_ref, Cpp_ref, Dpp_ref)
+
+  return
+end
+
+! ---
+
+subroutine diag_nonsym_right(n, A, A_ldim, V, V_ldim, energy, E_ldim)
+
+  implicit none
+
+  integer,          intent(in)  :: n, A_ldim, V_ldim, E_ldim
+  double precision, intent(in)  :: A(A_ldim,n)
+  double precision, intent(out) :: energy(E_ldim), V(V_ldim,n)
+
+  character*1                   :: JOBVL, JOBVR, BALANC, SENSE
+  integer                       :: i, j
+  integer                       :: ILO, IHI, lda, ldvl, ldvr, LWORK, INFO
+  double precision              :: ABNRM
+  integer,          allocatable :: iorder(:), IWORK(:)
+  double precision, allocatable :: WORK(:), SCALE_array(:), RCONDE(:), RCONDV(:)
+  double precision, allocatable :: Atmp(:,:), WR(:), WI(:), VL(:,:), VR(:,:), Vtmp(:)
+  double precision, allocatable :: energy_loc(:), V_loc(:,:)
+
+  !print*, " n = ", n
+  allocate(Atmp(n,n))
+  allocate(WI(n))
+  allocate(WR(n))
+  allocate(VR(n,n))
+  allocate(VL(1,1))
+
+  do i = 1, n
+    do j = 1, n
+      Atmp(j,i) = A(j,i)
+    enddo
+  enddo
+
+  JOBVL  = "N" ! computes the left  eigenvectors
+  JOBVR  = "V" ! computes the right eigenvectors
+  BALANC = "B" ! Diagonal scaling and Permutation for optimization
+  SENSE  = "V" ! Determines which reciprocal condition numbers are computed
+  lda  = n
+  ldvr = n
+  ldvl = 1
+
+  allocate(WORK(1), SCALE_array(n), RCONDE(n), RCONDV(n), IWORK(2*n-2))
+
+  LWORK = -1 ! to ask for the optimal size of WORK
+  call dgeevx(BALANC, JOBVL, JOBVR, SENSE,                  & ! CHARACTERS
+              n, Atmp, lda,                                 & ! MATRIX TO DIAGONALIZE
+              WR, WI,                                       & ! REAL AND IMAGINARY PART OF EIGENVALUES
+              VL, ldvl, VR, ldvr,                           & ! LEFT AND RIGHT EIGENVECTORS
+              ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV, & ! OUTPUTS OF OPTIMIZATION
+              WORK, LWORK, IWORK, INFO)
+
+  if(INFO .ne. 0) then
+    print*, 'dgeevx failed !!', INFO
+    stop
+  endif
+
+  LWORK = max(int(work(1)), 1) ! this is the optimal size of WORK
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+  call dgeevx(BALANC, JOBVL, JOBVR, SENSE,                  &
+              n, Atmp, lda,                                 &
+              WR, WI,                                       &
+              VL, ldvl, VR, ldvr,                           &
+              ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV, &
+              WORK, LWORK, IWORK, INFO)
+
+  if(INFO .ne. 0) then
+    print*, 'dgeevx failed !!', INFO
+    stop
+  endif
+
+  deallocate(WORK, SCALE_array, RCONDE, RCONDV, IWORK)
+  deallocate(VL, Atmp)
+
+
+  allocate(energy_loc(n), V_loc(n,n))
+  energy_loc = 0.d0
+  V_loc = 0.d0
+
+  i = 1
+  do while(i .le. n)
+
+    !print*, i, WR(i), WI(i)
+
+    if(dabs(WI(i)) .gt. 1d-7) then
+
+      print*, ' Found an imaginary component to eigenvalue'
+      print*, ' Re(i) + Im(i)', i, WR(i), WI(i)
+
+      energy_loc(i) = WR(i)
+      do j = 1, n
+        V_loc(j,i) = WR(i) * VR(j,i) - WI(i) * VR(j,i+1)
+      enddo
+      energy_loc(i+1) = WI(i)
+      do j = 1, n
+        V_loc(j,i+1) = WR(i) * VR(j,i+1) + WI(i) * VR(j,i)
+      enddo
+      i = i + 2
+
+    else
+
+      energy_loc(i) = WR(i)
+      do j = 1, n
+        V_loc(j,i) = VR(j,i)
+      enddo
+      i = i + 1
+
+    endif
+
+  enddo
+
+  deallocate(WR, WI, VR)
+
+
+  ! ordering
+!  do j = 1, n
+!    write(444, '(100(1X, F16.10))') (V_loc(j,i), i=1,5)
+!  enddo
+  allocate(iorder(n))
+  do i = 1, n
+    iorder(i) = i
+  enddo
+  call quick_sort(energy_loc, iorder, n)
+  do i = 1, n
+    energy(i) = energy_loc(i)
+    do j = 1, n
+      V(j,i) = V_loc(j,iorder(i))
+    enddo
+  enddo
+  deallocate(iorder)
+!  do j = 1, n
+!    write(445, '(100(1X, F16.10))') (V_loc(j,i), i=1,5)
+!  enddo
+  deallocate(V_loc, energy_loc)
+
+end
+
+! ---
+

src/utils/non_sym_diag.f90

@@ -519,7 +519,7 @@ subroutine impose_biorthog_svd(n, m, L, R)
   threshold               = 1.d-6
   num_linear_dependencies = 0
   do i = 1, m
-    if(abs(D(i)) <= threshold) then
+    if(dabs(D(i)) <= threshold) then
       D(i) = 0.d0
       num_linear_dependencies = num_linear_dependencies + 1
     else

src/utils/utils.f90

@@ -642,3 +642,113 @@ subroutine wall_time(t)
   t = dble(c)/dble(rate)
 end subroutine 
 
+! ---
+! Compute <u|u>
+double precision function u_dot_u(u, sze)
+
+  implicit none
+
+  integer,          intent(in) :: sze
+  double precision, intent(in) :: u(sze)
+
+  double precision, external   :: ddot
+
+  !DIR$ FORCEINLINE
+  u_dot_u = ddot(sze,u,1,u,1)
+
+end
+
+! ---
+
+! Orthogonalization using Q.R factorization
+subroutine ortho_qr(A, LDA, m, n)
+
+  ! A : matrix to orthogonalize
+  ! LDA : leftmost dimension of A
+  ! m : Number of rows of A
+  ! n : Number of columns of A
+  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+
+  integer,          intent(in)    :: m, n, LDA
+  double precision, intent(inout) :: A(LDA,n)
+
+  integer                        :: LWORK, INFO
+  double precision, allocatable  :: TAU(:), WORK(:)
+
+  allocate(TAU(min(m, n)))
+
+  allocate(WORK(1))
+
+  LWORK = -1
+
+  call dgeqrf(m, n, A(1,1), LDA, TAU, WORK, LWORK, INFO)
+
+  if(INFO .lt. 0) then
+    print*, ' dgeqrf failed !! ', INFO
+    stop
+  endif
+
+  LWORK = max(n, int(WORK(1)))
+
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+
+  call dgeqrf(m, n, A(1,1), LDA, TAU, WORK, LWORK, INFO)
+
+  if(INFO .lt. 0) then
+    print*, ' dgeqrf failed !! ', INFO
+    stop
+  endif
+
+  LWORK = -1
+  call dorgqr(m, n, n, A(1,1), LDA, TAU, WORK, LWORK, INFO)
+
+  if(INFO .lt. 0) then
+    print*, ' dorgqr failed !! ', INFO
+    stop
+  endif
+
+  LWORK = max(n, int(WORK(1)))
+ 
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+
+  call dorgqr(m, n, n, A(1,1), LDA, TAU, WORK, LWORK, INFO)
+
+  if(INFO .lt. 0) then
+    print*, ' dorgqr failed !! ', INFO
+    stop
+  endif
+
+  deallocate(WORK, TAU)
+
+end
+
+! ---
+
+! Normalizes vector u
+subroutine normalize(u, sze)
+
+  implicit none
+
+  integer,          intent(in)    :: sze
+  double precision, intent(inout) :: u(sze)
+
+  integer                         :: i
+  double precision                :: d
+  double precision, external      :: dnrm2
+
+  !DIR$ FORCEINLINE
+  d = dnrm2(sze, u, 1)
+  if (d /= 0.d0) then
+    d = 1.d0/d 
+  endif
+  if (d /= 1.d0) then
+    !DIR$ FORCEINLINE
+    call dscal(sze, d, u, 1)
+  endif
+
+end
+
+! ---
+