  We present a new wavefunction ansatz that combines the strengths of spin
projection with the language of matrix product states (MPS) and matrix product
operators (MPO) as used in the density matrix renormalization group (DMRG).
Specifically, spin-projected matrix product states (SP-MPS) are constructed as
$|\Psi^{(N,S,M)}_{SP-MPS}\rangle=\mathcal{P}_S|\Psi_{MPS}^{(N,M)}\rangle$,
where $\mathcal{P}_S$ is the spin projector for total spin $S$ and
$|\Psi_{MPS}^{(N,M)}\rangle$ is an MPS wavefunction with a given particle
number $N$ and spin projection $M$. This new ansatz possesses several
attractive features: (1) It provides a much simpler route to achieve
spin-adaptation (i.e. to create eigenfunctions of $\hat{S}^2$) compared to
explicitly incorporating the non-Abelian SU(2) symmetry into the MPS. In
particular, since the underlying state in the SP-MPS uses only Abelian
symmetries, one does not need the singlet embedding scheme for non-singlet
states, as normally employed in spin-adapted DMRG, to achieve a single
consistent variationally optimized state. (2) Due to the use of
$|\Psi_{MPS}^{(N,M)}\rangle$ as its underlying state, the SP-MPS can be closely
connected to broken symmetry mean-field states. This allows to
straightforwardly generate the large number of broken symmetry guesses needed
to explore complex electronic landscapes in magnetic systems. Further, this
connection can be exploited in the future development of quantum embedding
theories for open-shell systems. (3) The sum of MPOs representation for the
Hamiltonian and spin projector $\mathcal{P}_S$ naturally leads to an
embarrassingly parallel algorithm for computing expectation values and
optimizing SP-MPS. (4) Optimizing SP-MPS belongs to the
variation-after-projection (VAP) class of spin projected theories.
