Incorrectly translated equation in gaff

[verseon-it]

Apologies if this has already been addressed, since the code is very old

In antechamber/charmmgen.c routine empangle, (https://github.com/choderalab/ambermini/blob/master/antechamber/charmmgen.c#L1120) there are the following lines of code:

dparm = (bondlength1 - bondlength2) * (bondlength1 - bondlength2);
dparm =
	dparm / (bondlength1 + bondlength2) * (bondlength1 + bondlength2);
force =
	143.9 * zparm1 * cparm * zparm2 * exp(-2 * dparm) / (bondlength1 + bondlength2);
force /= sqrt(angle * 3.1415926 / 180.0);

First concern: you will notice that the second line of code is probably not performing the calculation intended by the original author.

Second concern: the last line is dividing by the square-root of the bend angle in radians. However, the correct operation (Eq. 5 in Gaff paper, J Comput Chem 25: 1157–1174, 2004) would be to divide by the square of the angle.

That said, employing the correct expressions results in values that don't appear correct. So perhaps this same lines of code were used in the original fit for the Z and C parameters. In which case correcting the code is not a good idea. But, perhaps a comment explaining the situation might be helpful for the next person that encounters this.