As demonstrated in the density matrix renormalization group (DMRG) method, approximating many-body wave function of electrons using a matrix product state (MPS) is a promising way to solve electronic structure problems. The expressibility of an MPS is determined by the size of the matrices or, in other words, the bond dimension, which unfortunately may be required to be very large in quantum chemistry simulations. In this study, we propose to calculate the ground state energies of molecular systems by variationally optimizing quantum circuit MPS (QCMPS) with a relatively small number of qubits. It is demonstrated that with carefully chosen circuit structure and orbital localization scheme, QCMPS can reach a similar accuracy as that achieved in DMRG with an exponentially large bond dimension. QCMPS simulation of a linear hydrogen molecular chain with 50 orbitals can reach the chemical accuracy using only 6 qubits at a moderate circuit depth. These results suggest that QCMPS is a promising wave function ansatz in the variational quantum eigensolver algorithm for molecular systems.