def METHOD_NAME():
    global mol, molsym, m, msym, mc_ref
    mol = gto.M(
    verbose = 5,
    output = '/dev/null',
    atom = [
        ["O" , (0. , 0.     , 0.)],
        [1   , (0. ,-0.757  , 0.587)],
        [1   , (0. , 0.757  , 0.587)]],
    basis = '631g',
    )
    m = scf.RHF(mol)
    m.conv_tol = 1e-10
    m.scf()
    molsym = gto.M(
    verbose = 5,
    output = '/dev/null',
    atom = [
        ["O" , (0. , 0.     , 0.)],
        [1   , (0. ,-0.757  , 0.587)],
        [1   , (0. , 0.757  , 0.587)]],
    basis = '631g',
    symmetry = True
    )
    msym = scf.RHF(molsym)
    msym.conv_tol = 1e-10
    msym.scf()
    mc_ref = mcscf.CASSCF (m, 4, 4).state_average_([0.25,]*4)
    # SA-CASSCF may be stuck at a local minimum e_tot = -75.75381945 with the
    # default initial guess from HF orbitals. The initial guess below is closed
    # to the single state CASSCF orbitals which can lead to a lower SA-CASSCF
    # energy e_tot = -75.762754627
    mo = mc_ref.sort_mo([4,5,6,10], base=1)
    mc_ref.kernel (mo)