We describe how to efficiently construct the quantum chemical Hamiltonian operator in matrix product form. We present its implementation as a density matrix renormalization group (DMRG) algorithm for quantum chemical applications in a purely matrix product based framework. Existing implementations of DMRG for quantum chemistry are based on the traditional formulation of the method, which was developed from a viewpoint of Hilbert space decimation and attained a higher performance compared to straightforward implementations of matrix product based DMRG. The latter variationally optimizes a class of ansatz states known as matrix product states (MPS), where operators are correspondingly represented as matrix product operators (MPO). The MPO construction scheme presented here eliminates the previous performance disadvantages while retaining the additional flexibility provided by a matrix product approach; for example, the specification of expectation values becomes an input parameter. In this way, MPOs for different symmetries - abelian and non-abelian - and different relativistic and non-relativistic models may be solved by an otherwise unmodified program.