The mean field approximation to the Ising model is a canonical variational\ntool that is used for analysis and inference in Ising models. We provide a\nsimple and optimal bound for the KL error of the mean field approximation for\nIsing models on general graphs, and extend it to higher order Markov random\nfields. Our bound improves on previous bounds obtained in work in the graph\nlimit literature by Borgs, Chayes, Lov\\'asz, S\\'os, and Vesztergombi and\nanother recent work by Basak and Mukherjee. Our bound is tight up to lower\norder terms. Building on the methods used to prove the bound, along with\ntechniques from combinatorics and optimization, we study the algorithmic\nproblem of estimating the (variational) free energy for Ising models and\ngeneral Markov random fields. For a graph $G$ on $n$ vertices and interaction\nmatrix $J$ with Frobenius norm $\\| J \\|_F$, we provide algorithms that\napproximate the free energy within an additive error of $\\epsilon n \\|J\\|_F$ in\ntime $\\exp(poly(1/\\epsilon))$. We also show that approximation within $(n\n\\|J\\|_F)^{1-\\delta}$ is NP-hard for every $\\delta > 0$. Finally, we provide\nmore efficient approximation algorithms, which find the optimal mean field\napproximation, for ferromagnetic Ising models and for Ising models satisfying\nDobrushin's condition.\n