A thermodynamically consistent, large-strain, multi-phase field approach (with consequent interface stresses) is generalized for the case with anisotropic interface (gradient) energy (e.g. an energy density that depends both on the magnitude and direction of the gradients in the phase fields). Such a generalization, if done in the “usual” manner, yields a theory that can be shown to be manifestly unphysical. These theories consider the gradient energy as anisotropic in the deformed configuration, and, due to this supposition, several fundamental contradictions arise. First, the Cauchy stress tensor is non-symmetric and, consequently, violates the moment of momentum principle, in essence the Herring (thermodynamic) torque is imparting an unphysical angular momentum to the system. In addition, this non-symmetric stress implies a violation of the principle of material objectivity. These problems in the formulation can be resolved by insisting that the gradient energy is an isotropic function of the gradient of the order parameters in the deformed configuration, but depends on the direction of the gradient of the order parameters (is anisotropic) in the undeformed configuration. We find that for a propagating nonequilibrium interface, the structural part of the interfacial Cauchy stress is symmetric and reduces to a biaxial tension with the magnitude equal to the temperature- and orientation-dependent interface energy. Ginzburg–Landau equations for the evolution of the order parameters and temperature evolution equation, as well as the boundary conditions for the order parameters are derived. Small strain simplifications are presented. Remarkably, this anisotropy yields a first order correction in the Ginzburg–Landau equation for small strains, which has been neglected in prior works. The next strain-related term is third order. For concreteness, specific orientation dependencies of the gradient energy coefficients are examined, using published molecular dynamics studies of cubic crystals. In order to consider a fully specified system, a typical sixth order polynomial phase field model is considered. Analytical solutions for the propagating interface and critical nucleus are found, accounting for the influence of the anisotropic gradient energy and elucidating the distribution of components of interface stresses. The orientation-dependence of the nonequilibrium interface energy is first suitably defined and explicitly determined analytically, and the associated width is also found. The developed formalism is applicable to melting/solidification and crystal-amorphous transformation and can be generalized for martensitic and diffusive phase transformations, twinning, fracture, and grain growth, for which interface energy depends on interface orientation of crystals from either side.