It could be done in PyMOL (step 1 needs to be performed only once):

1. Build a library of phosphorylated serine (3 rotamers), threonine (3), and tyrosine (6) residues.

2. Align the proper rotamer of the phosphorylated residue to the target residue for phosphorylation (align command).

3. Adjust the torsional angle involving the phosphate group if needed (set_dihedral).

4. Combine the coordinates from the target protein and the phosphorylated residue into PDB file (save).

The whole procedure can be automated with Python scripts, with or without PyMOL.

i would expect at least some of the commercial molecular graphics/modeling packages to have such a functionality, including the prediction of phosphorylation sites.

Dear Michael Bakker ,

No personal experience but you might consider Chimera: https://www.cgl.ucsf.edu/chimera/

Not sure how familiar you are with this program, but you might have a look at for example: http://plato.cgl.ucsf.edu/pipermail/chimera-users/2012-January/007068.html

Best regards.