<commit_before>module module_cc
  implicit none

! Energies           
  double precision                            :: Eccd    = 0.0d0
  double precision                            :: Ecc1    = 0.0d0
  double precision                            :: Ecc2    = 0.0d0
  double precision, allocatable               :: t2(:,:,:,:),t2old(:,:,:,:)
  double precision, allocatable               :: SO_Occ(:)

contains

!#############################################
!#         Calculate CC doubles
!#############################################
subroutine calc_Eccd(CCD,frac_occ)
  use module_data,      only : Spins,dim_1e,nOccA,nOccB,Eps,Occ
  use module_data,      only : DropMO,DoDrop,DoSingles,dim_1e,AMO
  use module_data,      only : DampCC,DoCCLS,CCMax
  use module_data,      only : Eps_SO,F_SO,DoReNorm,OmegaReg
  implicit none
  logical                   :: CCD 
  logical                   :: frac_occ
  logical                   :: MBPT3 = .false.
  integer                   :: iSpin,MOZero,iT2,nfocca,nfoccb,nfocc
  integer                   :: i,j,a,b 
  integer                   :: k,l,c,d

  integer                   :: nOcc,nmo
  integer                   :: maxT = 1000
  double precision          :: occ_factor,Dijab
  double precision          :: Ecctol = 1.0d-6
  double precision          :: LSmult = 1.0d4
  double precision          :: Eold   = 0.0d0
  double precision          :: deltaE = 99.0d0
  double precision          :: in1 ! intermediate during T Amp equations

  maxT = CCMax

  if (.not. DoCCLS) then
      LSmult=1.0d0
  else
      LSmult=1.0d4
  endif

  Eccd  = 0.0d0

  if(spins==1)then
    write(*,*) "    "
    write(*,*) "    Calculating RHF-(L)CCD correlation energy"
    write(*,*) "    "

    nmo = dim_1e*2
    MOZero = 1

    if(DoDrop)then
      MOZero = DropMO*2+1
    else
      MOZero = 1
    endif
    nOcc =  nOccA*2

    if(.not.allocated(t2))then
      allocate(t2(1:nmo,1:nmo,1:nmo,1:nmo),  &
            t2old(1:nmo,1:nmo,1:nmo,1:nmo))
    endif

    t2 = 0.0d0
    t2old = 0.0d0
    Dijab = 1.0d0
    Eccd = 0.0d0


  do iT2=1,maxT ! T amplitude loop
    Eold  = Eccd
    t2old  = t2 !R0.5d0*(t2+t2old) ! uses mixing

    !T2 equations
!$OMP PARALLEL PRIVATE(i,j,a,b,k,l,c,d,Dijab,in1)
!$OMP DO
    do i=MOZero,nOcc
      do j=MOZero,nOcc 
        do a=nOcc+1,nmo
          do b=nOcc+1,nmo
            Dijab = (Eps_SO(i) + Eps_SO(j) - Eps_SO(a) - Eps_SO(b))

!           [ + 1.0 ] * v ( p3 p4 h1 h2 )
            in1 = 1.0d0 * AMO(a,b,i,j)
            !write(*,*) i,j,a,b,in1

!           [ - 1.0 + 1.0 * P( p4 p3 h1 h2 => p3 p4 h1 h2 ) ] * Sum ( p5 ) * f ( p4 p5 ) * t ( p5 p3 h1 h2 ) done
            do d = nOcc+1,nmo
              if ((d == a) .or. (d == b)) then
                  cycle 
              endif
              in1 = in1 - 1.0d0 * F_SO(b,d) * t2old(a,d,i,j)
              in1 = in1 + 1.0d0 * F_SO(a,d) * t2old(d,b,i,j)
            enddo

!           [ - 1.0 + 1.0 * P( p3 p4 h1 h2 => p3 p4 h2 h1 ) ] * Sum ( h5 ) * f ( h5 h1 ) * t ( p3 p4 h5 h2 ) done
            do k = MOZero,nOcc
              if ((k == i) .or. (k == j)) then
                  cycle
              endif
              in1 = in1 - 1.0d0 * F_SO(k,j) * t2old(a,b,i,k)
              in1 = in1 + 1.0d0 * F_SO(k,i) * t2old(a,b,k,j)
            enddo

!           [ + 0.5 ] * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 ) * v ( p3 p4 p5 p6 ) done
            do c = nOcc+1,nmo
              do d = nOcc+1,c-1
                in1 = in1 + 1.0d0 * t2old(c,d,i,j) * AMO(a,b,c,d)
              enddo
            enddo

!           [ + 0.5 ] * Sum ( h5 h6 ) * t ( p3 p4 h5 h6 ) * v ( h5 h6 h1 h2 ) done
            do k = MOZero,nOcc
              do l = MOZero,k-1
                in1 = in1 + 1.0d0 * t2old(a,b,k,l) * AMO(k,l,i,j)
              enddo
            enddo

!          [ + 1.0 
!            - 1.0 * P( p3 p4 h2 h1 => p4 p3 h2 h1 ) 
!            - 1.0 * P( p3 p4 h2 h1 => p3 p4 h1 h2 ) 
!            + 1.0 * P( p3 p4 h2 h1 => p4 p3 h1 h2 ) ] 
!          * Sum ( p5 h6 ) * t ( p5 p3 h6 h2 ) * v ( h6 p4 h1 p5 ) done
            do c = nOcc+1,nmo
              do k = MOZero,nOcc
                in1 = in1 - 1.0d0 * t2old(a,c,i,k) * AMO(k,b,j,c)
                in1 = in1 + 1.0d0 * t2old(b,c,i,k) * AMO(k,a,j,c)
                in1 = in1 + 1.0d0 * t2old(a,c,j,k) * AMO(k,b,i,c)
                in1 = in1 - 1.0d0 * t2old(b,c,j,k) * AMO(k,a,i,c)
              enddo
            enddo

            if(CCD)then ! do full CCD
            !Terms that are quadratic in T2
              do c = nOcc+1,nmo
                do d = nOcc+1,nmo
                  do k = MOZero,nOcc
                    do l = MOZero,nOcc

                      in1 = in1 + 0.5d0 *  t2old(a,c,i,k)*t2old(d,b,l,j)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(b,c,i,k)*t2old(d,a,l,j)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(b,c,j,k)*t2old(d,a,l,i)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(a,c,j,k)*t2old(d,b,l,i)   * AMO(k,l,c,d)

                      in1 = in1 + 0.25d0*  t2old(c,d,i,j)*t2old(a,b,k,l)   * AMO(k,l,c,d)

                      in1 = in1 - 0.5d0 *  t2old(a,c,i,j)*t2old(b,d,k,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(b,c,i,j)*t2old(a,d,k,l)   * AMO(k,l,c,d)

                      in1 = in1 - 0.5d0 *  t2old(a,b,i,k)*t2old(c,d,j,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(a,b,j,k)*t2old(c,d,i,l)   * AMO(k,l,c,d)

                    enddo
                  enddo
                enddo
              enddo
            endif
<commit_msg>added factorization for closed shell cc<commit_after>module module_cc
  implicit none

! Energies           
  double precision                            :: Eccd    = 0.0d0
  double precision                            :: Ecc1    = 0.0d0
  double precision                            :: Ecc2    = 0.0d0
  double precision, allocatable               :: t2(:,:,:,:),t2old(:,:,:,:)
  double precision, allocatable               :: SO_Occ(:)

contains

!#############################################
!#         Calculate CC doubles
!#############################################
subroutine calc_Eccd(CCD,frac_occ)
  use module_data,      only : Spins,dim_1e,nOccA,nOccB,Eps,Occ
  use module_data,      only : DropMO,DoDrop,DoSingles,dim_1e,AMO
  use module_data,      only : DampCC,DoCCLS,CCMax
  use module_data,      only : Eps_SO,F_SO,DoReNorm,OmegaReg
  implicit none
  logical                   :: CCD 
  logical                   :: frac_occ
  logical                   :: MBPT3 = .false.
  integer                   :: iSpin,MOZero,iT2,nfocca,nfoccb,nfocc
  integer                   :: i,j,a,b 
  integer                   :: k,l,c,d

  integer                   :: nOcc,nmo
  integer                   :: maxT = 1000
  double precision          :: occ_factor,Dijab
  double precision          :: Ecctol = 1.0d-6
  double precision          :: LSmult = 1.0d4
  double precision          :: Eold   = 0.0d0
  double precision          :: deltaE = 99.0d0
  double precision          :: in1 ! intermediate during T Amp equations

  maxT = CCMax

  if (.not. DoCCLS) then
      LSmult=1.0d0
  else
      LSmult=1.0d4
  endif

  Eccd  = 0.0d0

  if(spins==1)then
    write(*,*) "    "
    write(*,*) "    Calculating RHF-(L)CCD correlation energy"
    write(*,*) "    "

    nmo = dim_1e*2
    MOZero = 1

    if(DoDrop)then
      MOZero = DropMO*2+1
    else
      MOZero = 1
    endif
    nOcc =  nOccA*2

    if(.not.allocated(t2))then
      allocate(t2(1:nmo,1:nmo,1:nmo,1:nmo),  &
            t2old(1:nmo,1:nmo,1:nmo,1:nmo))
    endif

    t2 = 0.0d0
    t2old = 0.0d0
    Dijab = 1.0d0
    Eccd = 0.0d0


  do iT2=1,maxT ! T amplitude loop
    Eold  = Eccd
    t2old  = t2 !R0.5d0*(t2+t2old) ! uses mixing

    !TODO:construct intermediates here
       if(factorized)then
       if(.not.allocated(X1))then
         allocate(X1(MOZero:nOcc,nOcc+1:nmo,nOcc+1:nmo,nOcc+1:nmo),  &
                  X2(MOZero:nOcc,MOZero:nOcc,MOZero:nOcc,MOZero:nOcc), &
                  X3(nOcc+1:nmo,nOcc+1:nmo), &  
                  X4(MOZero:nOcc,MOZero:nOcc))
       endif       

       do l=MOZero,nOcc
         do d=nOcc+1,nmo 
           do a=nOcc+1,nmo
             do i=nOcc+1,nmo
               Xin1 = 0.0d0
               do k=MOZero,nOcc
                 do c=nOcc+1,nmo 
                   Xin1 = Xin1 + t2old(a,c,i,k)*AMO(k,l,c,d)
                 enddo
               enddo
               X1(l,d,a,i) = Xin1
             enddo
           enddo
         enddo
       enddo

       do k=MOZero,nOcc
         do l=MOZero,nOcc
           do i=MOZero,nOcc
             do j=MOZero,nOcc
               Xin2 = 0.0d0
               do c=nOcc+1,nmo 
                 do d=nOcc+1,nmo 
                   Xin2 = Xin2 + t2old(c,d,i,j)*AMO(k,l,c,d)
                 enddo
               enddo
               X2(k,l,i,j) = Xin2
             enddo
           enddo
         enddo
       enddo

       do b=nOcc+1,nmo
         do c=nOcc+1,nmo
           Xin3 = 0.0d0
           do k=MOZero,nOcc
             do l=MOZero,nOcc
               do d=nOcc+1,nmo
                 Xin3 = Xin3 + t2old(b,d,k,l)*AMO(k,l,c,d)
               enddo
             enddo
           enddo
           X3(b,c) = Xin3
         enddo
       enddo


       do k=MOZero,nOcc 
         do j=MOZero,nOcc
           Xin4 = 0.0d0
           do l=MOZero,nOcc
             do c=nOcc+1,nmo 
               do d=nOcc+1,nmo
                 Xin4 = Xin4 + t2old(c,d,j,l)*AMO(k,l,c,d) 
               enddo
             enddo
           enddo
           X4(k,j) = Xin4
         enddo
       enddo

    endif

    !T2 equations
!$OMP PARALLEL PRIVATE(i,j,a,b,k,l,c,d,Dijab,in1)
!$OMP DO
    do i=MOZero,nOcc
      do j=MOZero,nOcc 
        do a=nOcc+1,nmo
          do b=nOcc+1,nmo
            Dijab = (Eps_SO(i) + Eps_SO(j) - Eps_SO(a) - Eps_SO(b))

!           [ + 1.0 ] * v ( p3 p4 h1 h2 )
            in1 = 1.0d0 * AMO(a,b,i,j)
            !write(*,*) i,j,a,b,in1

!           [ - 1.0 + 1.0 * P( p4 p3 h1 h2 => p3 p4 h1 h2 ) ] * Sum ( p5 ) * f ( p4 p5 ) * t ( p5 p3 h1 h2 ) done
            do d = nOcc+1,nmo
              if ((d == a) .or. (d == b)) then
                  cycle 
              endif
              in1 = in1 - 1.0d0 * F_SO(b,d) * t2old(a,d,i,j)
              in1 = in1 + 1.0d0 * F_SO(a,d) * t2old(d,b,i,j)
            enddo

!           [ - 1.0 + 1.0 * P( p3 p4 h1 h2 => p3 p4 h2 h1 ) ] * Sum ( h5 ) * f ( h5 h1 ) * t ( p3 p4 h5 h2 ) done
            do k = MOZero,nOcc
              if ((k == i) .or. (k == j)) then
                  cycle
              endif
              in1 = in1 - 1.0d0 * F_SO(k,j) * t2old(a,b,i,k)
              in1 = in1 + 1.0d0 * F_SO(k,i) * t2old(a,b,k,j)
            enddo

!           [ + 0.5 ] * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 ) * v ( p3 p4 p5 p6 ) done
            do c = nOcc+1,nmo
              do d = nOcc+1,c-1
                in1 = in1 + 1.0d0 * t2old(c,d,i,j) * AMO(a,b,c,d)
              enddo
            enddo

!           [ + 0.5 ] * Sum ( h5 h6 ) * t ( p3 p4 h5 h6 ) * v ( h5 h6 h1 h2 ) done
            do k = MOZero,nOcc
              do l = MOZero,k-1
                in1 = in1 + 1.0d0 * t2old(a,b,k,l) * AMO(k,l,i,j)
              enddo
            enddo

!          [ + 1.0 
!            - 1.0 * P( p3 p4 h2 h1 => p4 p3 h2 h1 ) 
!            - 1.0 * P( p3 p4 h2 h1 => p3 p4 h1 h2 ) 
!            + 1.0 * P( p3 p4 h2 h1 => p4 p3 h1 h2 ) ] 
!          * Sum ( p5 h6 ) * t ( p5 p3 h6 h2 ) * v ( h6 p4 h1 p5 ) done
            do c = nOcc+1,nmo
              do k = MOZero,nOcc
                in1 = in1 - 1.0d0 * t2old(a,c,i,k) * AMO(k,b,j,c)
                in1 = in1 + 1.0d0 * t2old(b,c,i,k) * AMO(k,a,j,c)
                in1 = in1 + 1.0d0 * t2old(a,c,j,k) * AMO(k,b,i,c)
                in1 = in1 - 1.0d0 * t2old(b,c,j,k) * AMO(k,a,i,c)
              enddo
            enddo

            if(CCD .and. .not. factorized)then ! do full CCD
            !Terms that are quadratic in T2
              do c = nOcc+1,nmo
                do d = nOcc+1,nmo
                  do k = MOZero,nOcc
                    do l = MOZero,nOcc

                      in1 = in1 + 0.5d0 *  t2old(a,c,i,k)*t2old(d,b,l,j)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(b,c,i,k)*t2old(d,a,l,j)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(b,c,j,k)*t2old(d,a,l,i)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(a,c,j,k)*t2old(d,b,l,i)   * AMO(k,l,c,d)

                      in1 = in1 + 0.25d0*  t2old(c,d,i,j)*t2old(a,b,k,l)   * AMO(k,l,c,d)

                      in1 = in1 - 0.5d0 *  t2old(a,c,i,j)*t2old(b,d,k,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(b,c,i,j)*t2old(a,d,k,l)   * AMO(k,l,c,d)

                      in1 = in1 - 0.5d0 *  t2old(a,b,i,k)*t2old(c,d,j,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(a,b,j,k)*t2old(c,d,i,l)   * AMO(k,l,c,d)

                    enddo
                  enddo
                enddo
              enddo
            elseif(CCD .and. factorized)then 
               ! Sum over l,d
                      in1 = in1 + 0.5d0 *  t2old(d,b,l,j) * X1(l,d,a,i)
                     !in1 = in1 + 0.5d0 *  t2old(a,c,i,k)*t2old(d,b,l,j)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(d,a,l,j) * X1(l,d,b,i)
                     !in1 = in1 - 0.5d0 *  t2old(b,c,i,k)*t2old(d,a,l,j)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(d,a,l,i) * X1(l,d,b,j) 
                     !in1 = in1 + 0.5d0 *  t2old(b,c,j,k)*t2old(d,a,l,i)   * AMO(k,l,c,d)
                      in1 = in1 - 0.5d0 *  t2old(d,b,l,i) * X1(l,d,a,j)
                     !in1 = in1 - 0.5d0 *  t2old(a,c,j,k)*t2old(d,b,l,i)   * AMO(k,l,c,d)

               ! Sum over k,l
                      in1 = in1 + 0.25d0*  t2old(a,b,k,l) * X2(k,l,i,j)
                     !in1 = in1 + 0.25d0*  t2old(c,d,i,j)*t2old(a,b,k,l)   * AMO(k,l,c,d)

               ! Sum over c
                      in1 = in1 - 0.5d0 *  t2old(a,c,i,j) * X3(b,c) 
                     !in1 = in1 - 0.5d0 *  t2old(a,c,i,j)*t2old(b,d,k,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(b,c,i,j) * X3(a,c)
                     !in1 = in1 + 0.5d0 *  t2old(b,c,i,j)*t2old(a,d,k,l)   * AMO(k,l,c,d)

               ! Sum over k
                      in1 = in1 - 0.5d0 *  t2old(a,b,i,k) * X4(k,j)
                     !in1 = in1 - 0.5d0 *  t2old(a,b,i,k)*t2old(c,d,j,l)   * AMO(k,l,c,d)
                      in1 = in1 + 0.5d0 *  t2old(a,b,j,k) * X4(k,i)
                     !in1 = in1 + 0.5d0 *  t2old(a,b,j,k)*t2old(c,d,i,l)   * AMO(k,l,c,d)

            endif

{'fraction_non_alphanumeric': 0.15623506927854755, 'fraction_numerical': 0.058767319636884856, 'mean_word_length': 1.2758620689655173, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4060222', 'n_tokens_mistral': 6709, 'n_tokens_neox': 6375, 'n_words': 1164}