CSH angle incorrect with GROMOS force field and virtual sites
The GROMOS force field uses the same atom type for H in COH (serine) and CSH (cysteine). Since constraint parameters, in this case for constraining the angle with virtual sites, are looked up based on the two atoms involved in the constraint, C and H in this case, we can't distinguish COH and CSH and the COH geometry is used. All other force fields supplied with Gromacs have an HS atom type for CSH, so there there is no issue.
The "simple" solution is to add HS to GROMOS, but that would modifiy the force field, which we should not do.
We could auto-generate constraint distances for angles (when parameters are not found), based on the lengths of the two other bonds and the equilibrium angle. This would get rid of some derived parameters for all force fields. But that requires some coding and I don't know if such auto-generation would be undesirable for certain other cases.
A third solution would be to require constraints=h-bonds with virtual sites, this effectively gives the previous solution without extra code. But currently we can't detect from a topology file if the user requested replacement of hydrogens by virtual sites.
#1 Updated by David van der Spoel over 6 years ago
Isn't running with virtual sites modifying the force field anyway?
In other words, how bad is it to add an identical atom type to the force field?
On the other hand, this does show that the implementation of the code is incorrect since it makes assumptions about the force field. In fact the grompp option
constraints = h-angles would generate the correct constraints. Could that be a solution?
#2 Updated by Berk Hess over 6 years ago
I meant to say "h-angles" iso "h-bonds" in my reply.
I think that technically that's the best solution, since it would also force h-angle constraints in any other molecule definitions, which the user might not have processed with pdb2gmx -vsite h. The problem with that is that we currently can't detect if (part of) a topology was generated with -vsite h and asking the user to set this is asking for mistakes. But maybe a check equivalent for the one I implemented for too fast bond vibrations for the time step could catch a constraints=h-bonds setting and suggest the user to use h-angles with vsites.