def METHOD_NAME(s1, s2, s12, mo):
    """Make IAOs for a molecule or single k-point"""
    s21 = s12.conj().T
    # s2 is overlap in minimal reference basis and should never be singular:
    s2cd = scipy.linalg.cho_factor(s2)
    ctild = scipy.linalg.cho_solve(s2cd, numpy.dot(s21, mo))
    try:
        s1cd = scipy.linalg.cho_factor(s1)
        p12 = scipy.linalg.cho_solve(s1cd, s12)
        ctild = scipy.linalg.cho_solve(s1cd, numpy.dot(s12, ctild))
    # s1 can be singular in large basis sets: Use canonical orthogonalization in this case:
    except numpy.linalg.LinAlgError:
        x = scf.addons.canonical_orth_(s1, lindep_threshold)
        p12 = numpy.linalg.multi_dot((x, x.conj().T, s12))
        ctild = numpy.dot(p12, ctild)
    # If there are no occupied orbitals at this k-point, all but the first term will vanish:
    if mo.shape[-1] == 0:
        return p12
    ctild = vec_lowdin(ctild, s1)
    ccs1 = numpy.linalg.multi_dot((mo, mo.conj().T, s1))
    ccs2 = numpy.linalg.multi_dot((ctild, ctild.conj().T, s1))
    #a is the set of IAOs in the original basis
    a = (p12 + 2*numpy.linalg.multi_dot((ccs1, ccs2, p12))
         - numpy.dot(ccs1, p12) - numpy.dot(ccs2, p12))
    return a