The three-dimensional (3-D) VMEC code has been modified to model an energetic species with a variant of a Bi-Maxwellian distribution function that satisfies the constraint B . del F-h = 0, and the 3-D TERPSICHORE stability code has been extended to investigate the effects of pressure anisotropy in two limits. The lower limit is based on a purely fluid Kruskal-Oberman (KO) energy principle (ignoring the stabilizing kinetic integral), and the upper limit is obtained from an energy principle in which the hot particle pressure and current density refrain from interacting with the dynamics of the instability because their diamagnetic drift frequency is considered much larger than the dominant growth rate. We have specifically investigated the instability properties of a Heliotron device with a major radius of 3.9 m and total (beta) similar or equal to 3.9%, where the energetic particle contribution (beta(h)) varies from 0 to 1.3% for T-parallel to/T-perpendicular to = 4. Both models demonstrate a significant stabilization of a global m/n = 4/2 mode as (beta(h))/(beta) approaches 1/3, with the noninteracting (NI) hot particle model virtually reaching the point of marginal stability when (beta(h))/(beta) similar or equal to 1/4. A variation of T-parallel to/T-perpendicular to at fixed (beta(h))/(beta) = 1/3 shows a noticeable decrease in the absolute magnitude of the unstable eigenvalue in the range 1 <= T-parallel to/T-perpendicular to <= 2 in the fluid KO model. The NI energetic particle model remains near marginality in the range 1 <= T-parallel to/T-perpendicular to <= 4. The Mercier stability is consistent with the global mode calculations. Performance benchmarks of the TERPSICHORE code on several different computers are presented.