There are plenty of different variants, particularly with regards to the conditions on the Hamiltonian. It's a bit of a game, for example, to try and find the simplest possible class of Hamiltonians for which simulation is still BQP-complete.

The statement will roughly be along the lines of: let $|psirangle$ be a (normalised) product state, $H$ be a Hamiltonian from some particular class (e.g. consisting only of nearest-neighbour couplings on a one-dimensional lattice), $hat O$ an observable comprising a tensor product of one-body operators such that $|hat O|leq 1$, and $t$ be a time. Given the promise that $langlepsi|e^{iHt}hat O e^{-iHt}|psirangle$ is either greater than $frac12+a$ or less than $frac12-a$ for some $a$ (e.g. $a=frac16$), decide which is the case.

---

Further Details

Hamiltonian Simulation is BQP-hard

The basic construction (originally due to Feynman, here tweaked a bit) basically shows how you can design a Hamiltonian that implements any quantum computation, including any BQP-complete computation. The observable you would measure is just $Z$ on a particular output qubit, the two measurement outcomes corresponding to 'yes' and 'no'.

The simplest sort of Hamiltonian you might think of is to consider a computation of $N-1$ sequential unitaries $U_n$ acting on $M$ qubits, starting from a state $|0rangle^{otimes M}$. Then you can introduce an extra $N$ qubits, and specify the Hamiltonian $$ H=frac{2}{N}sum_{n=1}^{N-1}sqrt{n(N-n)}left(|10ranglelangle 01|_{n,n+1}otimes U+|01ranglelangle 10|_{n,n+1}otimes U^daggerright). $$ If you prepare your initial state as $|1rangle|0rangle^{otimes(N-1)}|0rangle^{otimes M}$ then after a time $Npi/4$, it will be in a state $|0rangle^{otimes (N-1)}|1rangle|Phirangle$ where $|Phirangle$ is the output of the desired computation. The funny coupling strengths that I've used here, the $sqrt{n(N-n)}$ are chosen specifically to give deterministic evolution, and are related to the concept of perfect state transfer. Usually you'll see results stated with equal couplings, but probabilistic evolution.

To see how this works, you define a set of states $$ |psi_nrangle=|0rangle^{otimes(n-1)}|1rangle|0rangle^{otimes{N-n}}otimesleft(U_{n-1}U_{n-2}ldots U_1|0rangle^{otimes M}right). $$ The action of the Hamiltonian is then $$ H|psi_nrangle=frac2Nsqrt{(n-1)(N+1-n)}|psi_{n-1}rangle+frac2Nsqrt{n(N-n)}|psi_{n+1}rangle, $$ which proves that the evolution is restricted to an $Ntimes N$ subspace which is represented by a tridiagonal matrix (which is the specific thing studied in perfect state transfer).

Of course, this Hamiltonian doesn't have any particularly nice properties - it is highly non-local, for example. There are many tricks that can be played to simplify the Hamiltonian to being, for example, one-dimensional. It can even be translationally invariant if you want, at the cost of having to prepare a more complex initial product state (at that point, the computation is no longer encoded in the Hamiltonian, which is universal, but is encoded in the input state). See here, for example.

Hamiltonians can be simulated with a universal quantum computer

The evolution of any Hamiltonian which is local on some lattice, acting on an initial product state, for a time that is no more than polynomial in the system size, can be simulated by a quantum computer, and any efficiently implementable measurement can be applied to estimate an observable. In this sense, you can see that Hamiltonian simulation is no harder than a quantum computation, the counter-point to the previous statement that quantum computation is no harder than Hamiltonian simulation.

There are many ways to do this (and there have been some recent papers that show significant improvements in error scaling for certain classes of Hamiltonian). Hre's quite a simple one. Take the Hamiltonian $H$ that you want to simulate. Split it up into different parts, $H_i$, each of which commutes. For example, on a nearest-neighbour Hamiltonian on some graph, you don't need more pieces than the maximum degree of the graph. You then Trotterize the evolution, writing the approximation $$ e^{iHt}approx left(e^{-iH_1delta t}e^{-iH_2delta t}ldots e^{-iH_ndelta t}right)^{t/delta t} $$ So, you just have to construct a circuit that implements terms like $e^{-iH_1delta t}$, which is composed of commuting terms $H_1=sum_nh_n$, each of which acts only on a small number of qubits. $$ e^{-iH_1delta t}=prod_{n}e^{-ih_ndelta t} $$ Since this is just a unitary on a small number of terms, a universal quantum computer can implement it.