  Constructing matrix product operators (MPO) is at the core of the modern
density matrix renormalization group (DMRG) and its time dependent formulation.
For DMRG to be conveniently used in different problems described by different
Hamiltonians, in this work we propose a new generic algorithm to construct the
MPO of an arbitrary operator with a sum-of-products form based on the bipartite
graph theory. We show that the method has the following advantages: (i) It is
automatic in that only the definition of the operator is required; (ii) It is
symbolic thus free of any numerical error; (iii) The complementary operator
technique can be fully employed so that the resulting MPO is globally optimal
for any given order of degrees of freedom; (iv) The symmetry of the system
could be fully employed to reduce the dimension of MPO. To demonstrate the
effectiveness of the new algorithm, the MPOs of Hamiltonians ranging from the
prototypical spin-boson model and Holstein model to the more complicated ab
initio electronic Hamiltonian and the anharmonic vibrational Hamiltonian with
sextic force field are constructed. It is found that for the former three
cases, our automatic algorithm can reproduce exactly the same MPOs as the
optimally hand-crafted ones already known in the literature.
