It's no different than any simulation of a protein. pdb2gmx will handle multiple chains. The only thing you have to do is find a force field that supports phosphorylated residues. Most do (AMBER, CHARMM, GROMOS9 43A1p and likely some of the newer parameter sets).