def sto(zeta,N=1,L=0,M=0,origin=(0,0,0)):
    """
    Use Stewarts STO-6G fits to create a contracted Gaussian approximation to a
    Slater function. Fits of other expansion lengths (1G, 3G, etc) are in the paper.

    Reference: RF Stewart, JCP 52, 431 (1970)

    >>> s = sto(1)
    >>> np.isclose(s(0,0,0),0.530121)
    True
    """
    nlm2powers = {
        (1,0,0) : (0,0,0,0),   # x,y,z,r
        (2,0,0) : (0,0,0,1),
        (3,0,0) : (0,0,0,2),
        (2,1,0) : (1,0,0,0),
        (2,1,1) : (0,1,0,0),
        (2,1,-1) : (0,0,1,0),
        (3,1,0) : (1,0,0,1),
        (3,1,1) : (0,1,0,1),
        (3,1,-1) : (0,0,1,1)
        }

    gexps_1s = [2.310303149e01,4.235915534e00,1.185056519e00,
                4.070988982e-01,1.580884151e-01,6.510953954e-02]
    gcoefs_1s = [9.163596280e-03,4.936149294e-02,1.685383049e-01,
                 3.705627997e-01,4.164915298e-01,1.303340841e-01]
    
    gexps_2s = [2.768496241e01,5.077140627e00,1.426786050e00,
                2.040335729e-01,9.260298399e-02,4.416183978e-02]
    gcoefs_2s = [-4.151277819e-03,-2.067024148e-02,-5.150303337e-02,
                 3.346271174e-01,5.621061301e-01,1.712994697e-01]
    
    gexps_2p = [5.868285913e00,1.530329631e00,5.475665231e-01,
                2.288932733e-01,1.046655969e-01,4.948220127e-02]
    gcoefs_2p = [7.924233646e-03,5.144104825e-02,1.898400060e-01,
                 4.049863191e-01,4.012362861e-01,1.051855189e-01]
    
    gexps_3s = [3.273031938e00,9.200611311e-01,3.593349765e-01,
                8.636686991e-02,4.797373812e-02,2.724741144e-02]
    gcoefs_3s = [-6.775596947e-03,-5.639325779e-02,-1.587856086e-01,
                 5.534527651e-01,5.015351020e-01,7.223633674e-02]
    
    gexps_3p = [5.077973607e00,1.340786940e00,2.248434849e-01,
                1.131741848e-01,6.076408893e-02,3.315424265e-02]
    gcoefs_3p = [-3.329929840e-03,-1.419488340e-02,1.639395770e-01,
                 4.485358256e-01,3.908813050e-01,7.411456232e-02]
    gexps_3d = [2.488296923,7.981487853e-1,3.311327490e-1,
                1.559114463e-1,7.877734732e-2,4.058484363e-2]
    gcoefs_3d = [7.283828112e-3,5.386799363e-2,2.072139149e-1,
                 4.266269092e-1,3.843100204e-1,8.902827546e-2]
    
    gexps_4s = [3.232838646,3.605788802e-1,1.717902487e-1,
                5.277666487e-2,3.163400284e-2,1.874093091e-2]
    gcoefs_4s = [1.374817488e-3,-8.666390043e-2,-3.130627309e-1,
                 7.812787397e-1,4.389247988-1,2.487178756e-2]
    gexps_4p = [2.389722618, 7.960947826e-1,3.415541380e-1,
                8.847434525e-2,4.958248334e-2,2.816929784e-2]
    gcoefs_4p = [-1.665913575e-3,-1.657464971e-2,-5.958513378e-2,
                 4.053115554e-1,5.433958189e-1,1.20970491e-1]
    
    gexps = { # indexed by N,s_or_p:
        (1,0) : gexps_1s,
        (2,0) : gexps_2s,
        (2,1) : gexps_2p,
        (3,0) : gexps_3s,
        (3,1) : gexps_3p
        }

    gcoefs = {  # indexed by N,s_or_p:
        (1,0) : gcoefs_1s,
        (2,0) : gcoefs_2s,
        (2,1) : gcoefs_2p,
        (3,0) : gcoefs_3s,
        (3,1) : gcoefs_3p
        }

    I,J,K,R = nlm2powers[(N,L,M)]
    exps = [zeta**2*expn for expn in gexps[(N,L)]]
    coefs = gcoefs[N,L]
    return cgbf(origin,(I,J,K),exps,coefs)
Use Stewarts STO-6G fits to create a contracted Gaussian approximation to a
    Slater function. Fits of other expansion lengths (1G, 3G, etc) are in the paper.

    Reference: RF Stewart, JCP 52, 431 (1970)

    >>> s = sto(1)
    >>> np.isclose(s(0,0,0),0.530121)
    True