debug ppRPA

input/options

@@ -5,7 +5,7 @@
 # CC: maxSCF thresh   DIIS n_diis
       64     0.00001  F    1
 # CIS/TDHF/BSE: singlet triplet
-                T       F
+                T       T
 # GF: maxSCF thresh  DIIS n_diis renormalization
       64     0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta

src/QuAcK/linear_response_B_pp.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
+subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
 ! Compute the B matrix of the pp channel
 
@@ -8,14 +8,11 @@ subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 ! Input variables
 
   integer,intent(in)            :: ispin
-  logical,intent(in)            :: dRPA
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nOO,nVV
   double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
-  double precision              :: delta_spin
-  double precision              :: delta_dRPA
   double precision,external     :: Kronecker_delta
 
   integer                       :: a,b,i,j,ab,ij
@@ -24,33 +21,48 @@ subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
 
   double precision,intent(out)  :: B_pp(nVV,nOO)
 
-! Singlet or triplet manifold?
+! Build B matrix for the singlet manifold
  
-  delta_spin = 0d0
+  if(ispin == 1) then
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
-! Direct RPA
-
-  delta_dRPA = 0d0
-  if(dRPA) delta_dRPA = 1d0
-
-! Build A matrix
 
     ab = 0
     do a=nO+1,nBas-nR
-   do b=a+1,nBas-nR
+     do b=nO+1,a
         ab = ab + 1
         ij = 0
         do i=nC+1,nO
-       do j=i+1,nO
+         do j=nC+1,i
             ij = ij + 1
  
-          B_pp(ab,ij) = (1d0 + delta_spin)*ERI(a,b,i,j) - (1d0 - delta_dRPA)*ERI(a,b,j,i)
+            B_pp(ab,ij) = (ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
  
-        enddo
+          end do
-      enddo
+        end do
-    enddo
+      end do
-  enddo
+    end do
+
+  end if
+
+! Build the B matrix for the triplet manifold
+
+  if(ispin == 2) then
+
+    ab = 0
+    do a=nO+1,nBas-nR
+     do b=nO+1,a-1
+        ab = ab + 1
+        ij = 0
+        do i=nC+1,nO
+         do j=nC+1,i-1
+            ij = ij + 1
+ 
+            B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
 
 end subroutine linear_response_B_pp

src/QuAcK/linear_response_C_pp.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
+subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
 ! Compute the C matrix of the pp channel
 
@@ -8,14 +8,11 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 ! Input variables
 
   integer,intent(in)            :: ispin
-  logical,intent(in)            :: dRPA
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nOO,nVV
   double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
-  double precision              :: delta_spin
-  double precision              :: delta_dRPA
   double precision              :: eF
   double precision,external     :: Kronecker_delta
 
@@ -25,36 +22,54 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 
   double precision,intent(out)  :: C_pp(nVV,nVV)
 
-! Singlet or triplet manifold?
+! Define the chemical potential
-
-  delta_spin = 0d0
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
-! Direct RPA
-
-  delta_dRPA = 0d0
-  if(dRPA) delta_dRPA = 1d0
-
-! Build C matrix
  
   eF = e(nO) + e(nO+1)
 
+! Build C matrix for the singlet manifold
+
+  if(ispin == 1) then
+
     ab = 0
     do a=nO+1,nBas-nR
-   do b=a+1,nBas-nR
+     do b=nO+1,a
         ab = ab + 1
         cd = 0
         do c=nO+1,nBas-nR
-       do d=c+1,nBas-nR
+         do d=nO+1,c
             cd = cd + 1
  
             C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
-                        + (1d0 + delta_spin)*ERI(a,b,c,d) - (1d0 - delta_dRPA)*ERI(a,b,d,c)
+                          + (ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
  
-        enddo
+          end do
-      enddo
+        end do
-    enddo
+      end do
-  enddo
+    end do
+
+  end if
+
+! Build C matrix for the triplet manifold
+
+  if(ispin == 2) then
+
+    ab = 0
+    do a=nO+1,nBas-nR
+     do b=nO+1,a-1
+        ab = ab + 1
+        cd = 0
+        do c=nO+1,nBas-nR
+         do d=nO+1,c-1
+            cd = cd + 1
+ 
+            C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+                          + ERI(a,b,c,d) - ERI(a,b,d,c)
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
 
 end subroutine linear_response_C_pp

src/QuAcK/linear_response_D_pp.f90

@@ -1,4 +1,4 @@
-subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
+subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
 ! Compute the D matrix of the pp channel
 
@@ -8,15 +8,12 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 ! Input variables
 
   integer,intent(in)            :: ispin
-  logical,intent(in)            :: dRPA
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nOO,nVV
   double precision,intent(in)   :: e(nBas),ERI(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
   double precision              :: eF
-  double precision              :: delta_spin
-  double precision              :: delta_dRPA
   double precision,external     :: Kronecker_delta
 
   integer                       :: i,j,k,l,ij,kl
@@ -25,36 +22,54 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 
   double precision,intent(out)  :: D_pp(nOO,nOO)
 
-! Singlet or triplet manifold?
+! Define the chemical potential
-
-  delta_spin = 0d0
-  if(ispin == 1) delta_spin = +1d0
-  if(ispin == 2) delta_spin = -1d0
-
-! Direct RPA
-
-  delta_dRPA = 0d0
-  if(dRPA) delta_dRPA = 1d0
-
-! Build D matrix
 
   eF = e(nO) + e(nO+1)
  
+! Build the D matrix for the singlet manifold
+
+  if(ispin == 1) then
+
     ij = 0
     do i=nC+1,nO
-   do j=i+1,nO
+     do j=nC+1,i
         ij = ij + 1
         kl = 0
         do k=nC+1,nO
-       do l=k+1,nO
+         do l=nC+1,k
             kl = kl + 1
  
             D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
-                        + (1d0 + delta_spin)*ERI(k,l,i,j) - (1d0 - delta_dRPA)*ERI(k,l,j,i)
+                          + ERI(k,l,i,j) - ERI(k,l,j,i)
  
-        enddo
+          end do
-      enddo
+        end do
-    enddo
+      end do
-  enddo
+    end do
+
+  end if
+
+! Build the D matrix for the triplet manifold
+
+  if(ispin == 2) then
+
+    ij = 0
+    do i=nC+1,nO
+     do j=nC+1,i-1
+        ij = ij + 1
+        kl = 0
+        do k=nC+1,nO
+         do l=nC+1,k-1
+            kl = kl + 1
+ 
+            D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+                          + (ERI(k,l,i,j) + ERI(k,l,j,i))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
+ 
+          end do
+        end do
+      end do
+    end do
+
+  end if
 
 end subroutine linear_response_D_pp

src/QuAcK/linear_response_pp.f90

@@ -1,5 +1,4 @@
-subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & 
+subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA)
-                              Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA)
 
 ! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013)
 
@@ -8,8 +7,6 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,
 
 ! Input variables
 
-  logical,intent(in)            :: dRPA
-  logical,intent(in)            :: TDA
   logical,intent(in)            :: BSE
   integer,intent(in)            :: ispin,nBas,nC,nO,nV,nR
   integer,intent(in)            :: nOO
@@ -43,9 +40,9 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,
 
 ! Build B, C and D matrices for the pp channel
 
-  call linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B)
+  call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B)
-  call linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C)
+  call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C)
-  call linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D)
+  call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D)
 
 !------------------------------------------------------------------------
 ! Solve the p-p eigenproblem

src/QuAcK/ppRPA.f90

@@ -21,8 +21,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
 
 ! Local variables
 
-  logical                       :: dRPA
-  logical                       :: TDA
   logical                       :: BSE
   integer                       :: ispin
   integer                       :: nOO
@@ -44,45 +42,40 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
   write(*,*)'****************************************'
   write(*,*)
 
-! Useful quantities
-
- nOO = nO*(nO-1)/2
- nVV = nV*(nV-1)/2
-
 ! Initialization
 
   Ec_ppRPA(:) = 0d0
 
-! Switch on exchange for TDHF
-
-  dRPA = .false.
- 
-! Switch off Tamm-Dancoff approximation for TDHF
-
-  TDA = .false.
- 
 ! Switch off Bethe-Salpeter equation for TDHF
 
   BSE = .false. 
 
-! Memory allocation
-
-  allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
-           Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
-
 ! Singlet manifold
 
   if(singlet_manifold) then 
 
     ispin = 1
 
-    call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & 
+   ! Useful quantities
+
+    nOO = nO*(nO+1)/2
+    nVV = nV*(nV+1)/2
+
+   ! Memory allocation
+
+    allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
+             Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
+
+    call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,    & 
                             Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & 
                             Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & 
                             Ec_ppRPA(ispin))
+
     call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
     call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
 
+    deallocate(Omega1,X1,Y1,Omega2,X2,Y2)
+
   endif
 
 ! Triplet manifold 
@@ -91,13 +84,27 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
 
     ispin = 2
 
-    call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
+    ! Useful quantities
+
+    nOO = nO*(nO-1)/2
+    nVV = nV*(nV-1)/2
+
+  ! Memory allocation
+
+    allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & 
+             Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
+
+
+    call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,    &
                             Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & 
                             Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & 
                             Ec_ppRPA(ispin))
+
     call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))
     call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin))
 
+    deallocate(Omega1,X1,Y1,Omega2,X2,Y2)
+
   endif
 
   write(*,*)