# Source code for qcdb.libmintsgshell

    """Class to

    """
    import math

    #MAX_IOFF = 30000
    #extern size_t ioff[MAX_IOFF];
    #
    #MAX_DF = 500
    #extern double df[MAX_DF];
    #
    #MAX_BC = 20
    #extern double bc[MAX_BC][MAX_BC];
    #
    #MAX_FAC = 100
    #extern double fac[MAX_FAC];
    #
    #
    #MAX_DF = 500
    #extern double df[MAX_DF];
    #
    ## Globals
    #size_t ioff[MAX_IOFF];
    #double df[MAX_DF];
    #double bc[MAX_BC][MAX_BC];
    #double fac[MAX_FAC];
    #
    #def Wavefunction_initialize_singletons():
    #    done = False
    #
    #    if done:
    #        return
    #
    #    ioff[0] = 0;
    #    for (size_t i=1; i<MAX_IOFF; ++i)
    #        ioff[i] = ioff[i-1] + i;
    #
    #    df[0] = 1.0;
    #    df[1] = 1.0;
    #    df[2] = 1.0;
    #    for (int i=3; i<MAX_DF; ++i)
    #        df[i] = (i-1)*df[i-2];
    #
    #    for (int i=0; i<MAX_BC; ++i)
    #        for (int j=0; j<=i; ++j)
    #            bc[i][j] = combinations(i, j);
    #
    #    fac[0] = 1.0;
    #    for (int i=1; i<MAX_FAC; ++i)
    #        fac[i] = i*fac[i-1];
    #
    #    done = True

    def df(n):
        """Gives the double factorial of *i*"""
        return 1.0 if n <= 0 else 1.0 * n * df(n - 2)

    def INT_NCART(am):
        """Gives the number of cartesian functions for an angular momentum.
        define INT_NCART(am) ((am>=0) ? ((((am)+2)*((am)+1))>>1) : 0)

        """
        return (((am + 2) * (am + 1)) >> 1) if (am >= 0) else 0

    def INT_NPURE(am):
        """Gives the number of spherical functions for an angular momentum.
        #define INT_NPURE(am) (2*(am)+1)

        """
        return 2 * am + 1

    def INT_NFUNC(pu, am):
        """Gives the number of functions for an angular momentum based on pu.
        #define INT_NFUNC(pu,am) ((pu)?INT_NPURE(am):INT_NCART(am))

        """
        if pu == 'Cartesian' or pu == False:
            return INT_NCART(am)
        else:
            return INT_NPURE(am)

    def INT_CARTINDEX(am, i, j):
        """Computes offset index for cartesian function.
        #define INT_CARTINDEX(am,i,j) (((i) == (am))? 0 : (((((am) - (i) + 1)*((am) - (i)))>>1) + (am) - (i) - (j)))

        """
        return 0 if (i == am) else ((((am - i + 1) * (am - i)) >> 1) + am - i - j)

    def INT_ICART(a, b, c):
        """Given a, b, and c compute a cartesian offset.
        #define INT_ICART(a, b, c) (((((((a)+(b)+(c)+1)<<1)-(a))*((a)+1))>>1)-(b)-1)

        """
        return ((((((a + b + c + 1) << 1) - a) * (a + 1)) >> 1) - b - 1)

    def INT_IPURE(l, m):
        """Given l and m compute a pure function offset.
        #define INT_IPURE(l, m) ((l)+(m))

        """
        return l + m

    # Lookup array that when you index the angular momentum it returns the corresponding letter
    PrimitiveType = ['Normalized', 'Unnormalized']
    GaussianType = ['Cartesian', 'Pure']  # Cartesian = 0, Pure = 1

    class ShellInfo(object):
        """This class has the same behavior as GaussianShell, but implements everything using
        slower data structures, which are easier to construct. These are used to build the
        basis set, which builds more efficient pointer-based GaussianShell objects.
        *  @param e An array of exponent values.
        *  @param am Angular momentum.
        *  @param pure Pure spherical harmonics, or Cartesian.
        *  @param c An array of contraction coefficients.
        *  @param nc The atomic center that this shell is located on. Must map
            back to the correct atom in the owning BasisSet molecule_. Used
            in integral derivatives for indexing.
        *  @param center The x, y, z position of the shell. This is passed to
            reduce the number of calls to the molecule.
        *  @param start The starting index of the first function this shell
            provides. Used to provide starting positions in matrices.
        *  @param pt Is the shell already normalized?

        """

        def __init__(self, am, c, e, pure, nc, center, start, pt='Normalized'):
            # Angular momentum
            self.l = am
            # Flag for pure angular momentum (Cartesian = 0, Pure = 1)
            self.puream = pure
            # Exponents (of length nprimitives_)
            self.PYexp = e
            # Contraction coefficients (of length nprimitives_)
            self.PYcoef = c
            # ERD normalized contraction coefficients (of length nprimitives_)
            self.PYerd_coef = []
            # Original (un-normalized) contraction coefficients (of length nprimitives)
            self.PYoriginal_coef = [c[n] for n in range(len(c))]
            # Atom number this shell goes to. Needed when indexing integral derivatives.
            self.nc = nc
            # Atomic center number in the Molecule
            self.center = center
            #
            self.start = start
            # How many cartesian functions? (1=s, 3=p, 6=d, ...)
            self.PYncartesian = INT_NCART(self.l)
            # How many functions? (1=s, 3=p, 5/6=d, ...) * Dependent on the value of puream_
            self.PYnfunction = INT_NFUNC(self.puream, self.l)

            # Compute the normalization constants
            if pt == 'Unnormalized':
                self.normalize_shell()
                self.erd_normalize_shell()

        def primitive_normalization(self, p):
            """Normalizes a single primitive.
            @param p The primitive index to normalize.
            @return Normalization constant to be applied to the primitive.

            """
            tmp1 = self.l + 1.5
            g = 2.0 * self.PYexp[p]
            z = pow(g, tmp1)
            return math.sqrt((pow(2.0, self.l) * z) / (math.pi * math.sqrt(math.pi) * df(2 * self.l)))

        def contraction_normalization(self):
            """Normalizes an entire contraction set. Applies the normalization to the coefficients
            *  @param gs The contraction set to normalize.

            """
            e_sum = 0.0
            for i in range(self.nprimitive()):
                for j in range(self.nprimitive()):
                    g = self.PYexp[i] + self.PYexp[j]
                    z = pow(g, self.l + 1.5)
                    e_sum += self.PYcoef[i] * self.PYcoef[j] / z

            tmp = ((2.0 * math.pi / (2.0 / math.sqrt(math.pi))) * df(2 * self.l)) / pow(2.0, self.l)
            try:
                norm = math.sqrt(1.0 / (tmp * e_sum))
            except ZeroDivisionError:
                self.PYcoef[i] = [1.0 for i in range(self.nprimitive())]
            # Set the normalization
            for i in range(self.nprimitive()):
                self.PYcoef[i] *= norm

        def normalize_shell(self):
            """Handles calling primitive_normalization and
            contraction_normalization for you.

            """
            for i in range(self.nprimitive()):
                normalization = self.primitive_normalization(i)
                self.PYcoef[i] *= normalization
            self.contraction_normalization()

        def erd_normalize_shell(self):
            """

            """
            self.PYerd_coef = []
            tsum = 0.0
            for j in range(self.nprimitive()):
                for k in range(j + 1):
                    a1 = self.PYexp[j]
                    a2 = self.PYexp[k]
                    temp = self.PYoriginal_coef[j] * self.PYoriginal_coef[k]
                    temp2 = self.l + 1.5
                    temp3 = 2.0 * math.sqrt(a1 * a2) / (a1 + a2)
                    temp3 = pow(temp3, temp2)
                    temp *= temp3
                    tsum += temp
                    if j != k:
                        tsum += temp
            prefac = 1.0
            if self.l > 1:
                prefac = pow(2.0, 2 * self.l) / df(2 * self.l)
            norm = math.sqrt(prefac / tsum)
            for j in range(self.nprimitive()):
                self.PYerd_coef.append(self.PYoriginal_coef[j] * norm)

        def copy(self, nc=None, c=None):
            """Make a copy of the ShellInfo"""
            if nc is not None and c is not None:
                return ShellInfo(self.l, self.PYoriginal_coef, self.PYexp,
                    self.puream, nc, c,
                    self.start, 'Unnormalized')
            else:
                return ShellInfo(self.l, self.PYoriginal_coef, self.PYexp,
                    self.puream, self.nc, self.center,
                    self.start, 'Unnormalized')
            # better to just deepcopy?

        def nprimitive(self):
            """The number of primitive Gaussians"""
            return len(self.PYexp)

        def nfunction(self):
            """Total number of basis functions"""
            return INT_NFUNC(self.puream, self.l)

        def ncartesian(self):
            """Total number of functions if this shell was Cartesian"""
            return self.PYncartesian

        def am(self):
            """The angular momentum of the given contraction"""
            return self.l

        def amchar(self):
            """The character symbol for the angular momentum of the given contraction"""
            return 'spdfghiklmnopqrtuvwxyz'[self.l]

        def AMCHAR(self):
            """The character symbol for the angular momentum of the given contraction (upper case)"""
            return 'SPDFGHIKLMNOPQRTUVWXYZ'[self.l]

        def is_cartesian(self):
            """Returns true if contraction is Cartesian"""
            return True if self.puream == 'Cartesian' else False

        def is_pure(self):
            """Returns true if contraction is pure"""
            return True if self.puream == 'Pure' else False

        def center(self):
            """Returns the center of the Molecule this shell is on"""
            return self.center

        def ncenter(self):
            """Returns the atom number this shell is on. Used by integral derivatives for indexing."""
            return self.nc

        def exp(self, prim):
            """Returns the exponent of the given primitive"""
            return self.PYexp[prim]

        def coef(self, pi):
            """Return coefficient of pi'th primitive"""
            return self.PYcoef[pi]

        def erd_coef(self, pi):
            """Return ERD normalized coefficient of pi'th primitive"""
            return self.PYerd_coef[pi]

        def original_coef(self, pi):
            """Return unnormalized coefficient of pi'th primitive"""
            return self.PYoriginal_coef[pi]

        def exps(self):
            """Returns the exponent of the given primitive"""
            return self.PYexp

        def coefs(self):
            """Return coefficient of pi'th primitive and ci'th contraction"""
            return self.PYcoef

        def original_coefs(self):
            """Return unnormalized coefficient of pi'th primitive and ci'th contraction"""
            return self.PYoriginal_coef

        def pyprint(self, outfile=None):
            """Print out the shell"""
            text = """    %c %3d 1.00\n""" % (self.AMCHAR(), self.nprimitive())
            for K in range(self.nprimitive()):
                text += """               %20.8f %20.8f\n""" % (self.PYexp[K], self.PYoriginal_coef[K])

            if outfile is None:
                return text
            else:
                with open(outfile, mode='w') as handle:
                    handle.write(text)

        def __str__(self):
            """String representation of shell"""
            return self.pyprint(outfile=None)

        def normalize(self, l, m, n):
            """Normalize the angular momentum component"""
            return 1.0

        def function_index(self):
            """Basis function index where this shell starts."""
            return self.start

        def set_function_index(self, i):
            """Set basis function index where this shell starts."""
            self.start = i

    class GaussianShell(ShellInfo):
        """Class with same information as :py:class:`ShellInfo`. In C++,
        class uses more efficient data structures, but in Python differences
        minimal.

        """

        def __init__(self, am, nprimitive, oc, c, ec, e, pure, nc, center, start):
            """
            *  @param am Angular momentum.
            *  @param pure Pure spherical harmonics, or Cartesian.
            *  @param oc An array of contraction coefficients.
            *  @param c An array of normalized contraction coefficients.
            *  @param ec An array of ERD normalized contraction coefficients.
            *  @param e An array of exponent values.
            *  @param pure an enum describing whether this shell uses pure or Cartesian functions.
            *  @param nc The atomic center that this shell is located on. Must map back to the correct atom in the owning BasisSet molecule_. Used in integral derivatives for indexing.
            *  @param center The x, y, z position of the shell. This is passed to reduce the number of calls to the molecule.
            *  @param start The starting index of the first function this shell provides. Used to provide starting positions in matrices.
            *  @param pt Is the shell already normalized?

            """
            self.l = am
            self.PYnprimitive = nprimitive
            self.puream = pure
            self.PYexp = e
            self.PYoriginal_coef = oc
            self.PYcoef = c
            self.PYerd_coef = ec
            self.nc = nc
            self.center = center
            self.start = start
            self.PYncartesian = INT_NCART(self.l)
            self.PYnfunction = INT_NFUNC(self.puream, self.l)

        def nprimitive(self):
            """The number of primitive Gaussians"""
            return self.PYnprimitive

    #GaussianShell(0, nprimitive_,
    #    uoriginal_coefficients_, ucoefficients_, uerd_coefficients_,
    #    uexponents_, GaussianType(0), 0, xyz_, 0)
    #
    #GaussianShell(am, shell_nprim,
    #    &uoriginal_coefficients_[ustart+atom_nprim], &ucoefficients_[ustart+atom_nprim], &uerd_coefficients_[ustart+atom_nprim],
    #    &uexponents_[ustart+atom_nprim], puream, n, xyz_ptr, bf_count)
    #
    #GaussianShell(am, shell_nprim,
    #    &uoriginal_coefficients_[prim_count], &ucoefficients_[prim_count], &uerd_coefficients_[prim_count],
    #    &uexponents_[prim_count], puream, center, xyz_, bf_count)
    #
    #ShellInfo(am, contractions, exponents, gaussian_type, 0, center, 0, Unnormalized)