Free boundary three-dimensional anisotropic pressure magnetohydrodynamic equilibria with nested magnetic flux surfaces are computed through the minimisation of the plasma energy functional $W={\int}{V}{d^3}x\left[{B^2}/(2\mu_0)+p{||}/(\Gamma-1)\right]$. The plasma–vacuum interface is varied to guarantee the continuity of the total pressure $\left[{p}_{\perp}+{B^2}/(2\mu_0)\right]$ across it and the vacuum magnetic field must satisfy the Neumann boundary condition that its component normal to this interface surface vanishes. The vacuum magnetic field corresponds to that driven by the plasma current and external coils plus the gradient of a potential function whose solution is obtained using a Green's function method. The energetic particle contributions to the pressure are evaluated analytically from the moments of the variant of a bi-Maxwellian distribution function that satisfies the constraint ${\bf B\cdot\nabla}{\cal F}_h=0$. Applications to demonstrate the versatility and reliability of the numerical method employed have concentrated on high-β off-axis energetic particle deposition with large parallel and perpendicular pressure anisotropies in a 2-field period quasiaxisymmetric stellarator reactor system. For large perpendicular pressure anisotropy, the hot particle component of the pperpendicular distribution localises in the regions where the energetic particles are deposited. For large parallel pressure anisotropy, the pressures are more uniform around the flux surfaces.