So I am working on incorporating pseudocontact shifts into the Relax protocol and I need to implement the derivatives for minimization. As I understand it, to do this I need to implement the eval_atom_derivative function. However, I do not understand the vector parameters F1 and F2. Could you explain them to me? Are they just the vector along the gradient for the atom?
Thank you in advance.
The F1/F2 formulation of the derivitives comes from the formulation of Abe, Braun, Noguti and Go (1984)
"Rapid Calculation of First and Second Derivatives of Conformational Energy with Respect to Dihedral Angles for Proteins. General Recurrent Equations" Computers & Chemistry 8(4) pp. 239-247. (link) F1 and F2 correspond roughly to Fa and Ga, respectively, of equations 7a & 7b in that paper.
The key point to remember is that the Rosetta derivatives are written assuming torsional space minimization, so things are optimized to minimize over torsional degrees of freedom, rather than Cartesian degrees of freedom. This leads to derivatives which are a little difficult to understand, to say the least. Reading the Abe paper may or may not help you understand them.
My suggestion is to take one of the other existing terms as an example, and see how the derivatives were handled in that case. If you're looking for atom-wise derivatives, fa_atr and fa_rep would be obvious choices, but for optimization purposes they code is not obvious. fa_elec is possibly a better choice, although in the weekly releases it's been optimized, so it's less clear than it could be. You might want to take a look at the version in Rosetta3.5 or earlier, when it was called hack_elec (rosetta_source/src/core/scoring/hackelec/HackElecEnergy.cc). You can take a look at other terms, too, to see how they handle derivatives for their particular setup. The pair energy in src/core/scoring/methods/PairEnergy.cc and src/core/scoring/PairEPotential.cc might be helpful in this regard, as well.
One thing you definitely want to do when playing around with energy functions is to check to make sure that the analytical and calculated numeric derivatives match up. There are a number of functions in main/source/test/util/deriv_funcs.hh which are intended to help out with this, particularly the AtomDerivValidator and its simple_deriv_check() method. (These are used in test programs, and wouldn't be used in the actual calculations.)
Thank you for your answer.
You mentioned that the functions in main/source/test/util/deriv_funcs.hh to test the derivative calculations. What is the best way to use these functions?
You can take a look at some of the unit tests under the main/source/test/ directory to get an example of usage, but generally you would just instantiate an AtomDerivValidator object, passing an example Pose, a ScoreFunction with (just) your term turned on, and a MoveMap specifying which degrees of freedom you want to allow to move during minimization. Then you would simply call the simple_deriv_check() method of the AtomDerivValidator object. It will then check to make sure that the numeric and analytic derivatives are comparable around the current state of the Pose. If you want to get really thorough, you may want to try different Poses and different MoveMaps, to check if the derivitives make sense in different situations.
I tried instantiating an AtomDerivValidator by first including #include <../test/util/deriv_funcs.hh> However when I add this line into my code the compiler imidietely gives errors about TS_ASSERT and TS_ASSERT_DELTA not being declared in the scope. How can I fix this?
Those are defined in the CxxTest headers, usually included from the cxxtest/TestSuite.h header in main/source/external/cxxtest/cxxtest/TestSuite.h. You probably just want to add the main/source/external/cxxtest/ directory to your include path when you compile.
If that gets to be too complicated, you might just want to write your testing code in and run things through the Rosetta unit testing framework. See https://www.rosettacommons.org/docs/latest/development_documentation/test/run-unit-test and linked pages for an overview.
I am still working on implementing derivatives for pseudocontact shifts. I have been testing it by running a relax using only the PCS term in the energy function. However, I have always been getting an "Inaccurate G!" error after each time Rosetta calculates the derivatives. I do not fully understand what this means. Is it possible you can explain the reason for this error and would you know what might be causing it?
Thank you in advance and let me know if you need any aditional information.
We get "Inaccurate G!" errors with the standard Rosetta scorefunction, so it might not be anything to worry about.
What the message is indicating is that the calculated derivatives seems to be at odds with the change in energy values being observed during minimization. Some of this is down to inaccuracies/approximations in the minimizer itself (although we hope recent versions of Rosetta have fixed this), some of it is due to the functional form of the energy terms. (Histogram-based statistical potentials require smoothing/interpolation, and can have regions where they don't behave well for minimizers.)
As long as the numeric and analytical derivitives for your term are behaving well, I wouldn't worry about it too much. The one concern is that it's your derivatives which are causing the inaccuracy. I'd recommend doing parallel runs with/without the term, and seeing if the "Inaccurate G!" errors are exclusive (or are significantly more prevelent) to the runs with your term. If they are, you might want to look at the derivative implementation more closely, and conduct further tests.