We present the first of a two-part Mathematica notebook collection that implements a symbolic approach for the application of the density matrix expansion (DME) to the Hartree–Fock (HF) energy from a chiral effective field theory (EFT) three-nucleon interaction at N 2LO. The final output from the notebooks is a Skyrme-like energy density functional that provides a quasi-local approximation to the non-local HF energy. In this paper, we discuss the derivation of the HF energy and its simplification in terms of the scalar/vector–isoscalar/isovector parts of the one-body density matrix. Furthermore, a set of steps is described and illustrated on how to extend the approach to other three-nucleon interactions. Program summary Program title: SymbHFNNN Catalogue identifier: AEGC_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEGC_v1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 96 666 No. of bytes in distributed program, including test data, etc.: 378 083 Distribution format: tar.gz Programming language: Mathematica 7.1 Computer: Any computer running Mathematica 6.0 and later versions Operating system: Windows Xp, Linux/Unix RAM: 256 Mb Classification: 5, 17.16, 17.22 Nature of problem: The calculation of the HF energy from the chiral EFT three-nucleon interaction at N 2LO involves tremendous spin–isospin algebra. The problem is compounded by the need to eventually obtain a quasi-local approximation to the HF energy, which requires the HF energy to be expressed in terms of scalar/vector–isoscalar/isovector parts of the one-body density matrix. The Mathematica notebooks discussed in this paper solve the latter issue. Solution method: The HF energy from the chiral EFT three-nucleon interaction at N 2LO is cast into a form suitable for an automatic simplification of the spin–isospin traces. Several Mathematica functions and symbolic manipulation techniques are used to obtain the result in terms of the scalar/vector–isoscalar/isovector parts of the one-body density matrix. Running time: Several hours