def METHOD_NAME(mf, disp=1e-4, mo_rep=False):
    if not mf.converged: mf.METHOD_NAME()
    mo_coeff = np.asarray(mf.mo_coeff)
    RESTRICTED= (mo_coeff.ndim==3)
    cell = mf.cell
    cells_a, cells_b = gen_cells(cell, disp/2.0) # generate a bunch of cells with disp/2 on each cartesion coord
    mfset = run_mfs(mf, cells_a, cells_b) # run mean field calculations on all these cells
    vmat = get_vmat(mf, mfset, disp) # extracting <u|dV|v>/dR
    hmat = run_hess(mfset, disp)
    omega, vec = solve_hmat(cell, hmat)
    mass = cell.atom_mass_list() * MP_ME
    vec = _freq_mass_weighted_vec(vec, omega, mass)
    if mo_rep:
        if RESTRICTED:
            vmat = np.einsum('xuv,up,vq->xpq', vmat, mo_coeff[0].conj(), mo_coeff[0])
        else:
            vmat = np.einsum('xsuv,sup,svq->xspq', vmat, mo_coeff[:,0].conj(), mo_coeff[:,0])
    if RESTRICTED:
        mat = np.einsum('xJ,xpq->Jpq', vec, vmat)
    else:
        mat = np.einsum('xJ,xspq->sJpq', vec, vmat)
    return mat, omega