High order perturbation energies are computed for excited (1)A(1) states of BeH2 at geometries near the Be-->H-2 symmetric insertion transition state. The equations of multireference perturbation theory are solved through 30th order to study the difficulties in selecting the appropriate zeroth order Hamiltonian, orbitals, orbital energies, and reference functions for the computations of smooth molecular potential energy surfaces. The origin of the perturbative divergence produced by Moller-Plesset and Epstein-Nesbet partitionings is analyzed using a conceptually simple two-state model constructed using one state each from the reference and orthogonal spaces. The optimized zeroth order partitioning scheme (OPT) for double reference space computations with configurations 1a(1)(2)2a(1)(2)3a(1)(2) and 1a(1)(2)2a(1)(2)1b(2)(2) produces a truly convergent perturbation expansion through 30th order. The OPT energies are accurate in low orders as compared to the exact (197 dimensional) solution within the basis. The forced valence orbital degeneracy partitioning method (FD) also generates a truly convergent expansion for the same double reference space calculation, with slightly poorer low order energies than the OPT scheme. The BeH2 system facilitates the consideration of larger reference spaces (constructed using three through six orbitals) where the FD method produces highly accurate energies in low orders despite the asymptotic nature of the FD perturbation expansion. The "delayed" perturbative divergence behavior with the FD partitioning scheme (for large reference spaces) is shown to occur due to the incorrect ordering between the zeroth order energies of some reference and complementary space levels.