An energy conservative field directional scheme to approximate the solution of Gross-Pitaevskii systems of nonlinear equations in 2D is derived from the general idea of the Lie-Trotter splitting technique. In addition, the method preserves a numerical invariant which can be interpreted as a time step perturbation of the total mass. Under certain conditions, e.g., when the Josephson junction is neglected or when the initial conditions satisfy certain algebraic properties, mass conservation is also achieved. Using a suitable class of high-order symmetric finite difference approximations of the Laplacian operator in 2D, we prove that, for the thermodynamically stable regime, the directional splitting scheme is of order tau + h(2p) in the l(2,h)-norm, for p = 1, 2, 3, 4. High-order accuracy in time can be achieved by combining the proposed basic time step with well-known composition methods, as a result conservative methods of order tau(2q) + h(2p) with q = 1, 2, 3, 4 are obtained. Conservation and accuracy are numerically validated for model problems with and without internal atomic Josephson junction. Our performance study shows that the proposed technique is suitable for long-term simulations.