text: """
  Helper Classes for Configuration Interaction methods
  
  References:
  - Equations from [Szabo:1996]
  """
  
  __authors__ = "Tianyuan Zhang"
  __credits__ = ["Tianyuan Zhang", "Jeffrey B. Schriber", "Daniel G. A. Smith"]
  
  __copyright__ = "(c) 2014-2018, The Psi4NumPy Developers"
  __license__ = "BSD-3-Clause"
  __date__ = "2017-05-26"
  
  from itertools import combinations
  
  
  class Determinant:
      """
      A class for a bit-Determinant.
      """
  
      def __init__(self, alphaObtBits=0, betaObtBits=0, alphaObtList=None, betaObtList=None):
          """
          Constructor for the Determinant
          """
  
          if alphaObtBits == 0 and alphaObtList != None:
              alphaObtBits = Determinant.obtIndexList2ObtBits(alphaObtList)
          if betaObtBits == 0 and betaObtList != None:
              betaObtBits = Determinant.obtIndexList2ObtBits(betaObtList)
          self.alphaObtBits = alphaObtBits
          self.betaObtBits = betaObtBits
  
      def getNumOrbitals(self):
          """
          Return the number of orbitals (alpha, beta) in this determinant
          """
  
          return Determinant.countNumOrbitalsInBits(self.alphaObtBits), Determinant.countNumOrbitalsInBits(
              self.betaObtBits)
  
      def getOrbitalIndexLists(self):
          """
          Return lists of orbital index
          """
  
          return Determinant.obtBits2ObtIndexList(self.alphaObtBits), Determinant.obtBits2ObtIndexList(self.betaObtBits)
  
      def getOrbitalMixedIndexList(self):
          """
          Return lists of orbital in mixed spin index alternating alpha and beta
          """
  
          return Determinant.obtBits2ObtMixSpinIndexList(self.alphaObtBits, self.betaObtBits)
  
      @staticmethod
      def countNumOrbitalsInBits(bits):
          """
          Return the number of orbitals in this bits
          """
  
          count = 0
          while bits != 0:
              if bits & 1 == 1:
                  count += 1
              bits >>= 1
          return count
  
      @staticmethod
      def countNumOrbitalsInBitsUpTo4(bits):
          """
          Return the number of orbitals in this bits
          """
  
          count = 0
          while bits != 0 and count < 4:
              if bits & 1 == 1:
                  count += 1
              bits >>= 1
          return count
  
      @staticmethod
      def obtBits2ObtIndexList(bits):
          """
          Return the corresponding list of orbital numbers from orbital bits
          """
  
          i = 0
          obts = []
          while bits != 0:
              if bits & 1 == 1:
                  obts.append(i)
              bits >>= 1
              i += 1
          return obts
  
      @staticmethod
      def mixIndexList(alphaList, betaList):
          """
          Mix the alpha and beta orbital index list to one mixed list
          """
  
          return [elem * 2 for elem in alphaList] + [elem * 2 + 1 for elem in betaList]
  
      @staticmethod
      def obtBits2ObtMixSpinIndexList(alphaBits, betaBits):
          """
          Return the corresponding list of orbital numbers of orbital bits
          """
  
          alphaList, betaList = Determinant.obtBits2ObtIndexList(alphaBits), Determinant.obtBits2ObtIndexList(betaBits)
          return Determinant.mixIndexList(alphaList, betaList)
  
      @staticmethod
      def obtIndexList2ObtBits(obtList):
          """
          Return the corresponding orbital bits of list from orbital numbers
          """
  
          if len(obtList) == 0:
              return 0
          obtList = sorted(obtList, reverse=True)
          iPre = obtList[0]
          bits = 1
          for i in obtList:
              bits <<= iPre - i
              bits |= 1
              iPre = i
          bits <<= iPre
          return bits
  
      @staticmethod
      def getOrbitalPositions(bits, orbitalIndexList):
          """
          Return the position of orbital in determinant
          """
  
          count = 0
          index = 0
          positions = []
          for i in orbitalIndexList:
              while index < i:
                  if bits & 1 == 1:
                      count += 1
                  bits >>= 1
                  index += 1
              positions.append(count)
              continue
          return positions
  
      def getOrbitalPositionLists(self, alphaIndexList, betaIndexList):
          """
          Return the positions of indexes in lists
          """
  
          return Determinant.getOrbitalPositions(self.alphaObtBits, alphaIndexList), Determinant.getOrbitalPositions(
              self.betaObtBits, betaIndexList)
  
      def addAlphaOrbital(self, orbitalIndex):
          """
          Add an alpha orbital to the determinant
          """
  
          self.alphaObtBits |= 1 << orbitalIndex
  
      def addBetaOrbital(self, orbitalIndex):
          """
          Add an beta orbital to the determinant
          """
  
          self.betaObtBits |= 1 << orbitalIndex
  
      def removeAlphaOrbital(self, orbitalIndex):
          """
          Remove an alpha orbital from the determinant
          """
  
          self.alphaObtBits &= ~(1 << orbitalIndex)
  
      def removeBetaOrbital(self, orbitalIndex):
          """
          Remove an beta orbital from the determinant
          """
  
          self.betaObtBits &= ~(1 << orbitalIndex)
  
      def numberOfCommonOrbitals(self, another):
          """
          Return the number of common orbitals between this determinant and another determinant
          """
  
          return Determinant.countNumOrbitalsInBits(self.alphaObtBits &
                                                    another.alphaObtBits), Determinant.countNumOrbitalsInBits(
                                                        self.betaObtBits & another.betaObtBits)
  
      def getCommonOrbitalsInLists(self, another):
          """Return common orbitals between this determinant and another determinant in lists"""
          return Determinant.obtBits2ObtIndexList(self.alphaObtBits &
                                                  another.alphaObtBits), Determinant.obtBits2ObtIndexList(
                                                      self.betaObtBits & another.betaObtBits)
  
      def getCommonOrbitalsInMixedSpinIndexList(self, another):
          alphaList, betaList = self.getCommonOrbitalsInLists(another)
          return Determinant.mixIndexList(alphaList, betaList)
  
      def numberOfDiffOrbitals(self, another):
          """
          Return the number of different alpha and beta orbitals between this determinant and another determinant
          """
  
          diffAlpha, diffBeta = Determinant.countNumOrbitalsInBits(
              self.alphaObtBits ^ another.alphaObtBits), Determinant.countNumOrbitalsInBits(
                  self.betaObtBits ^ another.betaObtBits)
          return diffAlpha / 2, diffBeta / 2
  
      def numberOfTotalDiffOrbitals(self, another):
          """
          Return the number of different orbitals between this determinant and another determinant
          """
  
          diffAlpha, diffBeta = self.numberOfDiffOrbitals(another)
          return diffAlpha + diffBeta
  
      def diff2OrLessOrbitals(self, another):
          """
          Return true if two determinants differ 2 or less orbitals
          """
  
          diffAlpha, diffBeta = Determinant.countNumOrbitalsInBitsUpTo4(
              self.alphaObtBits ^ another.alphaObtBits), Determinant.countNumOrbitalsInBitsUpTo4(
                  self.betaObtBits ^ another.betaObtBits)
          return (diffAlpha + diffBeta) <= 4
  
      @staticmethod
      def uniqueOrbitalsInBits(bits1, bits2):
          """
          Return the unique bits in two different bits
          """
  
          common = bits1 & bits2
          return bits1 ^ common, bits2 ^ common
  
      @staticmethod
      def uniqueOrbitalsInLists(bits1, bits2):
          """
          Return the unique bits in two different bits
          """
  
          bits1, bits2 = Determinant.uniqueOrbitalsInBits(bits1, bits2)
          return Determinant.obtBits2ObtIndexList(bits1), Determinant.obtBits2ObtIndexList(bits2)
  
      def getUniqueOrbitalsInLists(self, another):
          """
          Return the unique orbital lists in two different determinants
          """
  
          alphaList1, alphaList2 = Determinant.uniqueOrbitalsInLists(self.alphaObtBits, another.alphaObtBits)
          betaList1, betaList2 = Determinant.uniqueOrbitalsInLists(self.betaObtBits, another.betaObtBits)
          return (alphaList1, betaList1), (alphaList2, betaList2)
  
      def getUnoccupiedOrbitalsInLists(self, nmo):
          """
          Return the unoccupied orbitals in the determinants
          """
  
          alphaBits = ~self.alphaObtBits
          betaBits = ~self.betaObtBits
          alphaObts = []
          betaObts = []
          for i in range(nmo):
              if alphaBits & 1 == 1:
                  alphaObts.append(i)
              alphaBits >>= 1
              if betaBits & 1 == 1:
                  betaObts.append(i)
              betaBits >>= 1
          return alphaObts, betaObts
  
      def getSignToMoveOrbitalsToTheFront(self, alphaIndexList, betaIndexList):
          """
          Return the final sign if move listed orbitals to the front
          """
  
          sign = 1
          alphaPositions, betaPositions = self.getOrbitalPositionLists(alphaIndexList, betaIndexList)
          for i in range(len(alphaPositions)):
              if (alphaPositions[i] - i) % 2 == 1:
                  sign = -sign
          for i in range(len(betaPositions)):
              if (betaPositions[i] - i) % 2 == 1:
                  sign = -sign
          return sign
  
      def getUniqueOrbitalsInListsPlusSign(self, another):
          """
          Return the unique orbital lists in two different determinants and the sign of the maximum coincidence determinants
          """
  
          alphaList1, alphaList2 = Determinant.uniqueOrbitalsInLists(self.alphaObtBits, another.alphaObtBits)
          betaList1, betaList2 = Determinant.uniqueOrbitalsInLists(self.betaObtBits, another.betaObtBits)
          sign1, sign2 = self.getSignToMoveOrbitalsToTheFront(alphaList1,
                                                              betaList1), another.getSignToMoveOrbitalsToTheFront(
                                                                  alphaList2, betaList2)
          return (alphaList1, betaList1), (alphaList2, betaList2), sign1 * sign2
  
      def getUniqueOrbitalsInMixIndexListsPlusSign(self, another):
          """
          Return the unique orbital lists in two different determinants and the sign of the maximum coincidence determinants
          """
  
          (alphaList1, betaList1), (alphaList2, betaList2), sign = self.getUniqueOrbitalsInListsPlusSign(another)
          return Determinant.mixIndexList(alphaList1, betaList1), Determinant.mixIndexList(alphaList2, betaList2), sign
  
      def toIntTuple(self):
          """
          Return a int tuple
          """
  
          return (self.alphaObtBits, self.betaObtBits)
  
      @staticmethod
      def createFromIntTuple(intTuple):
          return Determinant(alphaObtBits=intTuple[0], betaObtBits=intTuple[1])
  
      def generateSingleExcitationsOfDet(self, nmo):
          """
          Generate all the single excitations of determinant in a list
          """
  
          alphaO, betaO = self.getOrbitalIndexLists()
          alphaU, betaU = self.getUnoccupiedOrbitalsInLists(nmo)
          dets = []
  
          for i in alphaO:
              for j in alphaU:
                  det = self.copy()
                  det.removeAlphaOrbital(i)
                  det.addAlphaOrbital(j)
                  dets.append(det)
  
          for k in betaO:
              for l in betaU:
                  det = self.copy()
                  det.removeBetaOrbital(k)
                  det.addBetaOrbital(l)
                  dets.append(det)
  
          return dets
  
      def generateDoubleExcitationsOfDet(self, nmo):
          """
          Generate all the double excitations of determinant in a list
          """
  
          alphaO, betaO = self.getOrbitalIndexLists()
          alphaU, betaU = self.getUnoccupiedOrbitalsInLists(nmo)
          dets = []
  
          for i in alphaO:
              for j in alphaU:
                  for k in betaO:
                      for l in betaU:
                          det = self.copy()
                          det.removeAlphaOrbital(i)
                          det.addAlphaOrbital(j)
                          det.removeBetaOrbital(k)
                          det.addBetaOrbital(l)
                          dets.append(det)
  
          for i1, i2 in combinations(alphaO, 2):
              for j1, j2 in combinations(alphaU, 2):
                  det = self.copy()
                  det.removeAlphaOrbital(i1)
                  det.addAlphaOrbital(j1)
                  det.removeAlphaOrbital(i2)
                  det.addAlphaOrbital(j2)
                  dets.append(det)
  
          for k1, k2 in combinations(betaO, 2):
              for l1, l2 in combinations(betaU, 2):
                  det = self.copy()
                  det.removeBetaOrbital(k1)
                  det.addBetaOrbital(l1)
                  det.removeBetaOrbital(k2)
                  det.addBetaOrbital(l2)
                  dets.append(det)
          return dets
  
      def generateSingleAndDoubleExcitationsOfDet(self, nmo):
          """
          Generate all the single and double excitations of determinant in a list
          """
  
          return self.generateSingleExcitationsOfDet(nmo) + self.generateDoubleExcitationsOfDet(nmo)
  
      def copy(self):
          """
          Return a deep copy of self
          """
  
          return Determinant(alphaObtBits=self.alphaObtBits, betaObtBits=self.betaObtBits)
  
      def __str__(self):
          """
          Print a representation of the Determinant
          """
          a, b = self.getOrbitalIndexLists()
          return "|" + str(a) + str(b) + ">"
  
  
  import numpy as np
  
  
  class HamiltonianGenerator:
      """
      class for Full CI matrix elements
      """
  
      def __init__(self, H_spin, mo_spin_eri):
          """
          Constructor for MatrixElements
          """
  
          self.Hspin = H_spin
          self.antiSym2eInt = mo_spin_eri
  
      def generateMatrix(self, detList):
          """
          Generate CI Matrix
          """
  
          numDet = len(detList)
          matrix = np.zeros((numDet, numDet))
          for i in range(numDet):
              for j in range(i + 1):
                  matrix[i, j] = self.calcMatrixElement(detList[i], detList[j])
                  matrix[j, i] = matrix[i, j]
          return matrix
  
      def calcMatrixElement(self, det1, det2):
          """
          Calculate a matrix element by two determinants
          """
  
          numUniqueOrbitals = None
          if det1.diff2OrLessOrbitals(det2):
              numUniqueOrbitals = det1.numberOfTotalDiffOrbitals(det2)
              if numUniqueOrbitals == 0:
                  return self.calcMatrixElementIdentialDet(det1)
              if numUniqueOrbitals == 2:
                  return self.calcMatrixElementDiffIn2(det1, det2)
              elif numUniqueOrbitals == 1:
                  return self.calcMatrixElementDiffIn1(det1, det2)
              else:
                  # 
                  return 0.0
          else:
              return 0.0
  
      def calcMatrixElementDiffIn2(self, det1, det2):
          """
          Calculate a matrix element by two determinants where the determinants differ by 2 spin orbitals
          """
  
          unique1, unique2, sign = det1.getUniqueOrbitalsInMixIndexListsPlusSign(det2)
          return sign * self.antiSym2eInt[unique1[0], unique1[1], unique2[0], unique2[1]]
  
      def calcMatrixElementDiffIn1(self, det1, det2):
          """
          Calculate a matrix element by two determinants where the determinants differ by 1 spin orbitals
          """
  
          unique1, unique2, sign = det1.getUniqueOrbitalsInMixIndexListsPlusSign(det2)
          m = unique1[0]
          p = unique2[0]
          Helem = self.Hspin[m, p]
          common = det1.getCommonOrbitalsInMixedSpinIndexList(det2)
          Relem = 0.0
          for n in common:
              Relem += self.antiSym2eInt[m, n, p, n]
          return sign * (Helem + Relem)
  
      def calcMatrixElementIdentialDet(self, det):
          """
          Calculate a matrix element by two determinants where they are identical
          """
  
          spinObtList = det.getOrbitalMixedIndexList()
          Helem = 0.0
          for m in spinObtList:
              Helem += self.Hspin[m, m]
          length = len(spinObtList)
          Relem = 0.0
          for m in range(length - 1):
              for n in range(m + 1, length):
                  Relem += self.antiSym2eInt[spinObtList[m], spinObtList[n], spinObtList[m], spinObtList[n]]
          return Helem + Relem
  