I am in need of running molecular dynamics simulation for different protonation states of 1cwp.pdb (CCMV) over the range 4 to 7.

 One way of making the .psf file will be to alias the residues with protonation states given in the topology file and use guesscoord to add protons. But this is completely independent of the local structure of the protein.

I was thinking of making a .pqr file ( using pdb2pqr that has debumping algorithm and H-bond optimization, that makes protonation sensitive to the protein structure) and make a .psf out of it. This means aliasing  for a large number of  the added protons every time I do this, for different systems.

Please let me know if there is an efficient way of doing this.