tmpchem
commented
6 years ago Hi Mario. The important thing is that these analytical formulas match the values given by the numerical gradient, which they indeed have been confirmed to do. This is in fact not an error, but by necessary by design to get the correct result. This was actually one of the primary hold ups on this project for many months, as I was trying to get this formula for torsions and outofplanes and couldn't solve the math by hand until I consulted with an expert in molecular coordinate optimizations who recommended the formulas referenced below. Once I implemented those, the analytical and numerical gradients matched perfectly.
These formulas come from section 4-1 of "Molecular Vibrations" by Wilson, Decius, and Cross. They give both the direction and relative magnitude for these atomic gradient vectors for the various internal coordinate structures (bonds, angles, torsions, outofplanes). Essentially the reason for this division is that we aren't conserving the unit length of this vector, but rather the unity of the product of this vector and the bond length. You can think about it like torque or angular momentum. The longer the bond is, the greater the torque it generates due to the longer lever arm, which we must compensate for by having a smaller vector pulling it on the end.
https://github.com/tmpchem/computational_chemistry/blob/608b09255b36979043dc698b3f15a6d9646a17fd/scripts/molecular_mechanics/mmlib/gradient.py#L162
Why is the radius showing up in the fraction when the vectors are supposed to be normalized already? It's skewing up the direction of max increasing angle for the middle atom in the direction of the shortest side. What's the interpretation behind it?
PS: This also appears to be present in torsion calculations and oop angle gradients.