code: def getPoly(rCut, nMax):
      
      rCutVeryHard = rCut+5.0
      rx = 0.5*rCutVeryHard*(x + 1)
  
      basisFunctions = []
      for i in range(1, nMax + 1):
          basisFunctions.append(lambda rr, i=i, rCut=rCut: (rCut - np.clip(rr, 0, rCut))**(i+2))
  
      # Calculate the overlap of the different polynomial functions in a
      # matrix S. These overlaps defined through the dot product over the
      # radial coordinate are analytically calculable: Integrate[(rc - r)^(a
      # + 2) (rc - r)^(b + 2) r^2, {r, 0, rc}]. Then the weights B that make
      # the basis orthonormal are given by B=S^{-1/2}
      S = np.zeros((nMax, nMax))
      for i in range(1, nMax+1):
          for j in range(1, nMax+1):
              S[i-1, j-1] = (2*(rCut)**(7+i+j))/((5+i+j)*(6+i+j)*(7+i+j))
      betas = sqrtm(np.linalg.inv(S))
  
      # If the result is complex, the calculation is currently halted.
      if (betas.dtype == np.complex128):
          raise ValueError(
              "Could not calculate normalization factors for the polynomial basis"
              " in the domain of real numbers. Lowering the number of radial "
              "basis functions is advised."
          )
  
      fs = np.zeros([nMax, len(x)])
      for n in range(1, nMax+1):
          fs[n-1, :] = (rCut-np.clip(rx, 0, rCut))**(n+2)
  
      gss = np.dot(betas, fs)
  
      return nMax, rx, gss
docs: Used to calculate discrete vectors for the polynomial basis functions.
  
      Args:
          rCut(float): Radial cutoff
          nMax(int): Number of polynomial radial functions