CID and CISD

input/methods

@@ -7,7 +7,7 @@
 # drCCD rCCD crCCD lCCD
   F     F    F     F
 # CIS* CIS(D) CID CISD FCI
-  F    F      F   F    F
+  F    F      T   T    F
 # RPA* RPAx* crRPA ppRPA 
   F    F     F      F
 # G0F2* evGF2* qsGF2* G0F3 evGF3
@@ -15,7 +15,7 @@
 # G0W0* evGW* qsGW* ufG0W0 ufGW
   F    F    F       F      F
 # G0T0 evGT qsGT 
-  T    F    F
+  F    F    F
 # MCMP2
   F
 # * unrestricted version available

input/options

@@ -15,6 +15,6 @@
 # ACFDT: AC Kx  XBS
          F  F   F
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        T    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
  
 H 0. 0. 0.
-H 0. 0. 1.0
+H 0. 0. 2.000000

src/CI/CID.f90

@@ -0,0 +1,218 @@
+subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0)
+
+! Perform configuration interaction with doubles
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: triplet_manifold
+  integer,intent(in)            :: nBasin
+  integer,intent(in)            :: nCin
+  integer,intent(in)            :: nOin
+  integer,intent(in)            :: nVin
+  integer,intent(in)            :: nRin
+  double precision,intent(in)   :: Fin(nBasin,nBasin)
+  double precision,intent(in)   :: ERIin(nBasin,nBasin,nBasin,nBasin)
+  double precision,intent(in)   :: E0
+
+! Local variables
+
+  integer                       :: nBas
+  integer                       :: nC
+  integer                       :: nO
+  integer                       :: nV
+  integer                       :: nR
+
+  double precision,allocatable  :: F(:,:)
+  double precision,allocatable  :: sERI(:,:,:,:)
+  double precision,allocatable  :: ERI(:,:,:,:)
+
+  logical                       :: dump_trans = .false.
+  integer                       :: i,j,k,l
+  integer                       :: a,b,c,d
+  integer                       :: ia,kc,iajb,kcld
+  integer                       :: ishift,jshift
+  integer                       :: ispin
+  integer                       :: nD
+  integer                       :: nH
+  integer                       :: maxH
+  double precision,external     :: Kronecker_delta
+  double precision,allocatable  :: H(:,:)
+  double precision,allocatable  :: ECID(:)
+
+  double precision              :: tmp
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'******************************************************'
+  write(*,*)'| Configuration Interaction with Singles and Doubles |'
+  write(*,*)'******************************************************'
+  write(*,*)
+
+! Spatial to spin orbitals
+
+  nBas = 2*nBasin
+  nC   = 2*nCin
+  nO   = 2*nOin
+  nV   = 2*nVin
+  nR   = 2*nRin
+
+  allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas))
+
+  call spatial_to_spin_fock(nBasin,Fin,nBas,F)
+  call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
+
+! Antysymmetrize ERIs
+
+  allocate(ERI(nBas,nBas,nBas,nBas))
+
+  call antisymmetrize_ERI(2,nBas,sERI,ERI)
+
+  deallocate(sERI)
+
+! Compute CID matrix
+
+  nD  = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2
+  nH  = 1 + nD
+
+  write(*,*) 'nD = ',nD 
+  write(*,*) 'nH = ',nH
+  write(*,*)
+
+  maxH = min(nH,21)
+
+  ! Memory allocation
+
+  allocate(H(nH,nH),ECID(nH))
+ 
+  ! 00 block
+
+  ishift = 0 
+  jshift = 0
+
+  H(ishift+1,jshift+1) = E0
+ 
+  print*,'00 block done...'
+
+  ! 0D blocks
+
+  ishift = 0
+  jshift = 1
+
+  iajb = 0
+  do i=nC+1,nO
+    do a=1,nV-nR
+      do j=i+1,nO
+        do b=a+1,nV-nR
+
+          iajb = iajb + 1
+          tmp = ERI(i,j,nO+a,nO+b)
+
+          H(ishift+1,jshift+iajb) = tmp
+          H(jshift+iajb,ishift+1) = tmp
+
+        end do
+      end do
+    end do
+  end do
+  
+  print*,'0D blocks done...'
+
+  ! DD block
+
+  ishift = 1 
+  jshift = 1
+
+  iajb = 0
+  do i=nC+1,nO
+    do a=1,nV-nR
+      do j=i+1,nO
+        do b=a+1,nV-nR
+
+          iajb = iajb + 1
+
+          kcld = 0
+          do k=nC+1,nO
+            do c=1,nV-nR
+              do l=k+1,nO
+                do d=c+1,nV-nR
+       
+                  kcld = kcld + 1
+                  tmp = &
+                        E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+                      + F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
+                      - F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
+                      - F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+                      + F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
+                      - F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+                      + F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+                      + F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
+                      - F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+                      + F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      - F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      - F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      + F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      - F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+                      - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+                      + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+                      + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & 
+                      - ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
+                      - ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+                      + ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
+                      - ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
+                      - ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+                      - ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+                      + ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
+                      + ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
+                      - ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+                      + ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) &
+                      - ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) &
+                      - ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) &
+                      + ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) &
+                      - ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l)   
+
+                  H(ishift+iajb,jshift+kcld) = tmp
+       
+                end do
+              end do
+            end do
+          end do
+
+        end do
+      end do
+    end do
+  end do
+
+  print*,'DD block done...'
+
+  write(*,*)
+  write(*,*) 'Diagonalizing CID matrix...'
+  write(*,*)
+
+  call diagonalize_matrix(nH,H,ECID)
+
+  print*,'CID energies (au)'
+  call matout(maxH,1,ECID)
+  write(*,*)
+
+  print*,'CID excitation energies (eV)'
+  call matout(maxH-1,1,(ECID(2:maxH)-ECID(1))*HaToeV)
+  write(*,*)
+
+  if(dump_trans) then
+    print*,'Singlet CID transition vectors'
+    call matout(nH,nH,H)
+    write(*,*)
+  endif
+
+end subroutine CID

src/CI/CISD.f90

@@ -1,4 +1,4 @@
-subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
+subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0)
 
 ! Perform configuration interaction with singles and doubles
 
@@ -9,23 +9,42 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
 
   logical,intent(in)            :: singlet_manifold
   logical,intent(in)            :: triplet_manifold
-  integer,intent(in)            :: nBas,nC,nO,nV,nR
+  integer,intent(in)            :: nBasin
-  double precision,intent(in)   :: eHF(nBas)
+  integer,intent(in)            :: nCin
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  integer,intent(in)            :: nOin
+  integer,intent(in)            :: nVin
+  integer,intent(in)            :: nRin
+  double precision,intent(in)   :: Fin(nBasin,nBasin)
+  double precision,intent(in)   :: ERIin(nBasin,nBasin,nBasin,nBasin)
+  double precision,intent(in)   :: E0
 
 ! Local variables
 
+  integer                       :: nBas
+  integer                       :: nC
+  integer                       :: nO
+  integer                       :: nV
+  integer                       :: nR
+
+  double precision,allocatable  :: F(:,:)
+  double precision,allocatable  :: sERI(:,:,:,:)
+  double precision,allocatable  :: ERI(:,:,:,:)
+
   logical                       :: dump_trans = .false.
   integer                       :: i,j,k,l
   integer                       :: a,b,c,d
-  integer                       :: ia,jb,iajb,kcld
+  integer                       :: ia,kc,iajb,kcld
   integer                       :: ishift,jshift
   integer                       :: ispin
   integer                       :: nS
   integer                       :: nD
-  integer                       :: nSD
+  integer                       :: nH
+  integer                       :: maxH 
   double precision,external     :: Kronecker_delta
-  double precision,allocatable  :: H(:,:),Omega(:)
+  double precision,allocatable  :: H(:,:)
+  double precision,allocatable  :: ECISD(:)
+
+  double precision              :: tmp
 
 ! Hello world
 
@@ -35,170 +54,257 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF)
   write(*,*)'******************************************************'
   write(*,*)
 
-! Compute CIS matrix
+! Spatial to spin orbitals
 
-  if(singlet_manifold) then
+  nBas = 2*nBasin
+  nC   = 2*nCin
+  nO   = 2*nOin
+  nV   = 2*nVin
+  nR   = 2*nRin
 
-    ispin = 1
+  allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas))
 
-    ! Dimensions
+  call spatial_to_spin_fock(nBasin,Fin,nBas,F)
+  call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
 
-    nS  = (nO - nC)*(nV - nR)
+! Antysymmetrize ERIs
-    nD  = (nO - nC)*(nO - nC + 1)/2*(nV - nR)*(nV - nR + 1)/2
-    nSD = 1 + nS + nD
 
-    print*,'nS  = ',nS 
+  allocate(ERI(nBas,nBas,nBas,nBas))
-    print*,'nD  = ',nD 
-    print*,'nSD = ',nSD 
 
-    ! Memory allocation
+  call antisymmetrize_ERI(2,nBas,sERI,ERI)
 
-    allocate(H(nSD,nSD),Omega(nSD))
+  deallocate(sERI)
 
-    ! 0D block
+! Compute CISD matrix
 
-    ishift = 0
+  nS = (nO - nC)*(nV - nR)
-    jshift = 1 + nS
+  nD = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2
+  nH = 1 + nS + nD
 
-    iajb = 0
+  write(*,*) 'nS = ',nS 
+  write(*,*) 'nD = ',nD 
+  write(*,*) 'nH = ',nH
+  write(*,*)
 
-    do i=nC+1,nO
+  maxH = min(nH,21)
-      do a=1,nV-nR
-        do j=i,nO
-          do b=a,nV-nR
 
-            iajb = iajb + 1
+  ! Memory allocation
-            H(ishift+1,jshift+iajb) = ERI(i,j,nO+a,nO+b)
+
+  allocate(H(nH,nH),ECISD(nH))
+ 
+  ! 00 block
+
+  ishift = 0 
+  jshift = 0
+
+  H(ishift+1,jshift+1) = E0
+ 
+  print*,'00 block done...'
+
+  ! 0S blocks
+
+  ishift = 0 
+  jshift = 1
+
+  ia = 0
+  do i=nC+1,nO
+    do a=1,nV-nR
+
+      ia = ia + 1
+      tmp = F(i,nO+a) 
+      H(ishift+1,jshift+ia) = tmp
+      H(jshift+ia,ishift+1) = tmp
+
+    end do
+  end do
+
+  print*,'0S blocks done...'
+
+  ! 0D blocks
+
+  ishift = 0
+  jshift = 1 + nS
+
+  iajb = 0
+  do i=nC+1,nO
+    do a=1,nV-nR
+      do j=i+1,nO
+        do b=a+1,nV-nR
+
+          iajb = iajb + 1
+          tmp = ERI(i,j,nO+a,nO+b)
+
+          H(ishift+1,jshift+iajb) = tmp
+          H(jshift+iajb,ishift+1) = tmp
 
-          end do
         end do
       end do
     end do
+  end do
   
-    ! SS block
+  print*,'0D blocks done...'
 
-    ishift = 1
+  ! SS block
-    jshift = 1
 
-    ia = 0
+  ishift = 1
-    jb = 0
+  jshift = 1
 
-    do i=nC+1,nO
+  ia = 0
-      do a=1,nV-nR
+  do i=nC+1,nO
+    do a=1,nV-nR
 
-        ia = ia + 1
+      ia = ia + 1
+      kc = 0
+      do k=nC+1,nO
+        do c=1,nV-nR
 
-        do j=nC+1,nO
+          kc = kc + 1
-          do b=1,nV-nR
+          tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
+              - F(i,k)*Kronecker_delta(a,c)                  &
+              + F(nO+a,nO+c)*Kronecker_delta(i,k)            &
+              - ERI(nO+a,k,nO+c,i)
 
-            jb = jb + 1
+          H(ishift+ia,jshift+kc) = tmp
 
-            H(ishift+ia,jshift+jb) & 
-              = Kronecker_delta(i,j)*Kronecker_delta(a,b)*(eHF(nO+a) - eHF(i)) & 
-              + ERI(nO+a,j,i,nO+b) -  ERI(nO+a,j,nO+b,i)
-
-          end do
         end do
-
       end do
+
     end do
+  end do
 
-    ! SD block
+  print*,'SS block done...'
 
-    ishift = 1 
+  ! SD blocks
-    jshift = 1 + nS
 
-    ia = 0
+  ishift = 1 
-    kcld = 0
+  jshift = 1 + nS
 
-    do i=nC+1,nO
+  ia = 0
-      do a=1,nV-nR
+  do i=nC+1,nO
+    do a=1,nV-nR
 
-        ia = ia + 1
+      ia = ia + 1
+      kcld = 0
 
-        do k=nC+1,nO
+      do k=nC+1,nO
-          do c=1,nV-nR
+        do c=1,nV-nR
-            do l=k,nO
+          do l=k+1,nO
-              do d=c,nV-nR
+            do d=c+1,nV-nR
 
-                kcld = kcld + 1
+              kcld = kcld + 1
+              tmp = - F(l,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+                    + F(l,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
+                    - F(k,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+                    + F(k,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
+                    - ERI(k,l,nO+d,i)*Kronecker_delta(a,c)                &
+                    + ERI(k,l,nO+c,i)*Kronecker_delta(a,d)                &
+                    - ERI(nO+a,l,nO+c,nO+d)*Kronecker_delta(i,k)          &
+                    + ERI(nO+a,k,nO+c,nO+d)*Kronecker_delta(i,l)           
                   
-                H(ishift+ia,jshift+kcld)                                               & 
+              H(ishift+ia,jshift+kcld) = tmp
-                = Kronecker_delta(i,k)*(ERI(nO+a,l,nO+c,nO+d) - ERI(nO+a,l,nO+d,nO+c)) &
+              H(jshift+kcld,ishift+ia) = tmp
-                - Kronecker_delta(i,l)*(ERI(nO+a,k,nO+c,nO+d) - ERI(nO+a,k,nO+d,nO+c)) &
-                - Kronecker_delta(a,c)*(ERI(k,l,i,nO+d)       - ERI(k,l,nO+d,i))       &
-                + Kronecker_delta(a,d)*(ERI(k,l,i,nO+c)       - ERI(k,l,nO+c,i))
 
-              end do
             end do
           end do
         end do
-
       end do
+
     end do
+  end do
 
-    ! DD block
+  print*,'SD blocks done...'
 
-    ishift = 1 + nS
+  ! DD block
-    jshift = 1 + nS
 
-    iajb = 0
+  ishift = 1 + nS
-    kcld = 0
+  jshift = 1 + nS
 
-    do i=nC+1,nO
+  iajb = 0
-      do a=1,nV-nR
+  do i=nC+1,nO
-        do j=i,nO
+    do a=1,nV-nR
-          do b=a,nV-nR
+      do j=i+1,nO
+        do b=a+1,nV-nR
 
-            iajb = iajb + 1
+          iajb = iajb + 1
 
-            do k=nC+1,nO
+          kcld = 0
-              do c=1,nV-nR
+          do k=nC+1,nO
-                do l=k,nO
+            do c=1,nV-nR
-                  do d=c,nV-nR
+              do l=k+1,nO
+                do d=c+1,nV-nR
        
-                    kcld = kcld + 1
+                  kcld = kcld + 1
+                  tmp = &
+                        E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & 
+                      + F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
+                      - F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
+                      - F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+                      + F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
+                      - F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+                      + F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+                      + F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
+                      - F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+                      + F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      - F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      - F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      + F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      - F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                      - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & 
+                      - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+                      + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+                      + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & 
+                      - ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
+                      - ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+                      + ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
+                      - ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
+                      - ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+                      - ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
+                      + ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
+                      + ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
+                      - ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+                      + ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) &
+                      - ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) &
+                      - ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) &
+                      + ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) &
+                      - ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+                      + ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l)   
 
-!                   H(ishift+iajb,jshift+kcld)                            & 
+                  H(ishift+iajb,jshift+kcld) = tmp
-!                   = Kronecker_delta(i,k)*(ERI(a,l,c,d) -  ERI(a,l,d,c)) &
-!                   - Kronecker_delta(i,l)*(ERI(a,k,c,d) -  ERI(a,k,d,c)) &
-!                   - Kronecker_delta(a,c)*(ERI(k,l,i,d) -  ERI(k,l,d,i)) &
-!                   + Kronecker_delta(a,d)*(ERI(k,l,i,c) -  ERI(k,l,c,i))
        
-                  end do
                 end do
               end do
             end do
-
           end do
+
         end do
       end do
     end do
+  end do
 
-    call diagonalize_matrix(nSD,H,Omega)
+  print*,'DD block done...'
-    call print_excitation('CISD  ',ispin,nS,Omega)
 
-    if(dump_trans) then
+  write(*,*)
-      print*,'Singlet CISD transition vectors'
+  write(*,*) 'Diagonalizing CISD matrix...'
-      call matout(nSD,nSD,H)
+  write(*,*)
-      write(*,*)
-    endif
 
+  call diagonalize_matrix(nH,H,ECISD)
+
+  print*,'CISD energies (au)'
+  call matout(maxH,1,ECISD)
+  write(*,*)
+
+  print*,'CISD excitation energies (eV)'
+  call matout(maxH-1,1,(ECISD(2:maxH)-ECISD(1))*HaToeV)
+  write(*,*)
+
+  if(dump_trans) then
+    print*,'Singlet CISD transition vectors'
+    call matout(nH,nH,H)
+    write(*,*)
   endif
 
-! if(triplet_manifold) then
-
-!   ispin = 2
-!
-!   call diagonalize_matrix(nSD,H,Omega)
-!   call print_excitation('CISD  ',ispin,nSD,Omega)
-
-!   if(dump_trans) then
-!     print*,'Triplet CIS transition vectors'
-!     call matout(nSD,nSD,H)
-!     write(*,*)
-!   endif
-
-! endif
-
 end subroutine CISD

src/GT/dynamic_Tmatrix_A.f90

@@ -59,12 +59,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
 
           do cd=1,nVV
             eps = + Omega1(cd) 
-            chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
+            chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2)
           end do
 
           do kl=1,nOO
             eps = - Omega2(kl)
-            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+            chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
           end do
 
           A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
@@ -73,12 +73,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
 
           do cd=1,nVV
             eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
-            chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) 
+            chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2) 
           end do
 
           do kl=1,nOO
             eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
-            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+            chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
           end do
 
           A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi
@@ -87,12 +87,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
 
           do cd=1,nVV
             eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j))
-            chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+            chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
           end do
 
           do kl=1,nOO
             eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b))
-            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+            chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
           end do
 
           ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 1d0*lambda*chi

src/GT/static_Tmatrix_A.f90

@@ -47,13 +47,13 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
           do cd=1,nVV
             eps = + Omega1(cd)
 !           chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
-            chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
+            chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2)
           enddo
 
           do kl=1,nOO
             eps = - Omega2(kl)
 !           chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
-            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
+            chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2)
           enddo
 
           TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi

src/GT/static_Tmatrix_B.f90

@@ -47,13 +47,13 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh
           do cd=1,nVV
             eps = + Omega1(cd)
 !           chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
-            chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2)
+            chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
           enddo
 
           do kl=1,nOO
             eps = - Omega2(kl)
 !           chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
-            chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
+            chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
           enddo
 
           TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi

src/HF/RHF.f90

@@ -1,4 +1,4 @@
-subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx)
+subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx)
 
 ! Perform restricted Hartree-Fock calculation
 
@@ -45,7 +45,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
   double precision,allocatable  :: J(:,:)
   double precision,allocatable  :: K(:,:)
   double precision,allocatable  :: cp(:,:)
-  double precision,allocatable  :: F(:,:)
   double precision,allocatable  :: Fp(:,:)
   double precision,allocatable  :: ON(:)
 
@@ -56,6 +55,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
   double precision,intent(out)  :: c(nBas,nBas)
   double precision,intent(out)  :: P(nBas,nBas)
   double precision,intent(out)  :: Vx(nBas)
+  double precision,intent(out)  :: F(nBas,nBas)
 
 ! Hello world
 
@@ -72,7 +72,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
 ! Memory allocation
 
   allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),        &
-           cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), &
+           cp(nBas,nBas),Fp(nBas,nBas),ON(nBas),              &
            error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
 
 ! Guess coefficients and eigenvalues

src/QuAcK/QuAcK.f90

@@ -64,6 +64,8 @@ program QuAcK
   double precision,allocatable  :: dipole_int_MO(:,:,:)
   double precision,allocatable  :: dipole_int_aa(:,:,:)
   double precision,allocatable  :: dipole_int_bb(:,:,:)
+  double precision,allocatable  :: F_AO(:,:)
+  double precision,allocatable  :: F_MO(:,:)
   double precision,allocatable  :: ERI_AO(:,:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: ixyz
@@ -247,7 +249,7 @@ program QuAcK
 
   allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
            S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), &
-           dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin))
+           dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin),F_AO(nBas,nBas))
 
 ! Read integrals
 
@@ -291,7 +293,7 @@ program QuAcK
 
     call cpu_time(start_HF)
     call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
-             nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
+             nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
     call cpu_time(end_HF)
 
     t_HF = end_HF - start_HF
@@ -433,6 +435,7 @@ program QuAcK
       ! Memory allocation
      
       allocate(ERI_MO(nBas,nBas,nBas,nBas))
+      allocate(F_MO(nBas,nBas))
  
       ! Read and transform dipole-related integrals
     
@@ -448,7 +451,8 @@ program QuAcK
       ket1 = 1
       ket2 = 1
       call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO)
-!     call AOtoMO_transform(nBas,cHF,T+V)
+      F_MO(:,:) = F_AO(:,:)
+      call AOtoMO_transform(nBas,cHF,F_MO)
     end if
 
   end if
@@ -734,7 +738,7 @@ program QuAcK
   if(doCID) then
 
     call cpu_time(start_CID)
-!   call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
+    call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF)
     call cpu_time(end_CID)
 
     t_CID = end_CID - start_CID
@@ -750,7 +754,7 @@ program QuAcK
   if(doCISD) then
 
     call cpu_time(start_CISD)
-!   call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF)
+    call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF)
     call cpu_time(end_CISD)
 
     t_CISD = end_CISD - start_CISD

src/utils/spatial_to_spin_fock.f90

@@ -0,0 +1,29 @@
+subroutine spatial_to_spin_fock(nBas,F,nBas2,sF)
+
+! Convert Fock matrix from spatial to spin orbitals
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nBas,nBas2
+  double precision,intent(in)   :: F(nBas,nBas)
+
+! Local variables
+
+  integer                       :: p,q
+  double precision,external     :: Kronecker_delta
+
+! Output variables
+
+  double precision,intent(out)  :: sF(nBas2,nBas2)
+
+    do p=1,nBas2
+      do q=1,nBas2
+
+        sF(p,q) = Kronecker_delta(mod(p,2),mod(q,2))*F((p+1)/2,(q+1)/2)
+
+      enddo
+    enddo
+
+end subroutine spatial_to_spin_fock