Question about gradient of energy function

[Anonymous]

Hi,

Is it possible to get the gradient (i.e. the derivative with respect to each independent degree of freedom) of the energy score evaluated on a given protein pose? It seems like the Minmover subroutine moves the pose in a small step in the direction of this gradient, meaning the gradient in implicit in the Rosetta code. I’d like to access the value of the gradient itself at any fixed pose, for example in terms of dihedral coordinates or Cartesian coordinates. Technically one could numerically approximate the gradient by perturbing each degree of freedom (say dihedral angle) one at a time and seeing how the energy function changes, but this is a very inefficient solution given the high dimensionality of the situation.

Since the formula for energy is very complicated and probably a nightmare to differentiate by hand, I’d be quite content to have the energy function in a complete explicit form so that I could enter it into an automatic differentiation software package (of which there are many good ones these days). So far most of the references I’ve seen in the Rosetta literature give a high-level explanation of each term but not quite an explicit form (or the various constants and so on needed to get that explicit form).

Thanks!

Kyler

[Anonymous]

I’m not going to actually answer your question, but I’ll explain a little of how the code works.

The score function is a container full of EnergyMethods.  Each EnergyMethod will provide a score AND the derivatives.  So, Rosetta already has analytic derivatives baked in for most scorefunction terms.  (There are a few that either can’t have derivatives, or don’t bother to define them; those obviously perform poorly in minimization).  You don’t need any help outside of Rosetta to get the gradients.

As to how to access the derivatives in PyRosetta – that I don’t know.  I would crack open the minimizer in the C++ code to start figuring it out, which clearly isn’t helpful to you.  There are wiki pages on https://www.rosettacommons.org/docs/wiki/rosetta_basics/structural_concepts/minimization-overview and https://www.rosettacommons.org/docs/wiki/rosetta_basics/scoring/scoring-explained but they don’t really help at the code leve.

[Anonymous]

This isn’t something that has easy access in Rosetta. Part of the reason for that is that Rosetta’s minimization machinery uses the formalism of Abe, Braun, Noguti and Go (1984 Computers & Chemistry 8(4) pp. 239-247) for the efficient calculations of derivatives with respect to torsional minimization. This optimization means that the native derivatives of EnergyMethod objects (F1/F2 vectors) are less than easily interpretable.

I’m guessing that to get some of the information you need, you’re going to have to pull apart some of the minimization machinery.

The first thing you want to do is set up the minimization:

scorefxn = get_score_function()



move_map = core.kinematics.MoveMap()  # Adjust heavily to pull out just the degrees of freedom you want

minoptions = core.optimization.MinimizerOptions( "linmin", 0.001, True ) # Other settings may be of interest



start_score = scorefxn(pose) # Important to make sure the scorefunction and pose are coordinated before proceeding



min_map = core.optimization.MinimizerMap() # Need to change this to CartesianMinimizerMap for cartmin

min_map.setup( pose, move_map )

scorefxn.setup_for_minimizing( pose, min_map );



# AtomTreeMultifunc is for internal-coordinate (lengths/angles/dihedrals) degrees of freedom.

# There's also a CartesianMultifunc for Cartesian (x/y/z) degrees of freedom.

multifunc = core.optimization.AtomTreeMultifunc( pose, min_map, scorefxn )



dofs = core.optimization.Multivec( min_map.nangles() ) # ndofs() for cartmin

min_map.copy_dofs_from_pose( pose, dofs );



start_func_val = multifunc( dofs )

To this point this is basically a minimial alteration of how the standard minimizer sets things up. At this point we skip the minimization process, and just pull out the derivatives.

vars = core.optimization.Multivec( dofs )

dE_dvars = core.optimization.Multivec( dofs )



multifunc.dfunc( vars, dE_dvars )

At this point, dE_dvars (which is actually just a utility::vector1< Real>) should contain the analytical derivatives for each degree of freedom. If you have more than one DOF enabled in the minimization, you should be able to call min_map.dof_nodes() to get a parallel list of DOF details. This object should have details about the DOF (including residue&atom number and the DOF_type [torsion/angle/length])

That’s for the standard (internal coordinate) MinimizerMap. The CartesianMinimizerMap works slightly differently, where you can call min_map.get_atom(n) to get the core.id.AtomID (residue&atom number) for the atom which corresponds to the  3N-2, 3N-1 and 3Nth entry (x, y, and z, respecitively, I believe) in the dE_dvars Multivec.

 

Note that I haven’t tested *any* of this in a PyRosetta context, and am only marginally sure of how this work in the C++ context. There’s likely to be a large amount of fiddling to get things to work. I’d also highly recommend that you run a fair number of tests (e.g. by comparing simple systems with manually computed numerical derivatives) to make sure that the numbers you’re getting out match what you think they should be and there isn’t a bug in your implementation or in the descriptions I gave to you.

[Anonymous]

Thanks for the replies! @rmoretti, I wasn’t able to get any of the code with Multivec to work. Do you know how I can import it or otherwise get it to run?

[Anonymous]

What sort of errors are you getting?   If it’s simply because the core.optimization.Multivec name isn’t recognized, you can juse use the PyRosetta equivalent of `utility::vector1< Real >`, which I think would be as simple as

dofs = rosetta.Vector1( [0.0]* min_map.nangles() )



...



vars = rosetta.Vector1( list( dofs ) )

dE_dvars = rosetta.Vector1( list( dofs ) )

 

[Anonymous]

Great, that works! I checked at least the numerial partial derivatives in each coordinate and they seem to be consistent with what dE_dvars outputs.

[Anonymous]

Actually, while Rosetta’s derivatives seem to agree with the numerical derivatives for the full-atom score function, they don’t seem to agree anymore when I use the centroid version. Any reason why that might be? Here’s the pyrosetta code I’m trying.

from pyrosetta import *

init()

p = pose_from_pdb(‘1qlq.clean.pdb’)

switch = SwitchResidueTypeSetMover(‘centroid’)

switch.apply(p)

from pyrosetta.teaching import *

sf = get_cen_scorefxn()

start_score = sf(p)

movemap = MoveMap()

movemap.set_bb(True)

minmap = rosetta.core.optimization.MinimizerMap()

minmap.setup(p,movemap)

sf.setup_for_minimizing(p,minmap)

multifunc = rosetta.core.optimization.AtomTreeMultifunc(p,minmap,sf)

num_dofs = minmap.nangles()

x = Vector1([0.0]*num_dofs)

minmap.copy_dofs_from_pose(p,x)

E = multifunc(x)

dE_dx = Vector1([0.0]*num_dofs)

multifunc.dfunc(x,dE_dx)

print ‘rosetta dE_dx: ‘

print dE_dx

eps = 0.001

numerical_derivs = []

for i in range(100):

    x_pert = Vector1(list(x))

    x_pert[1+i] += eps

    E_pert = multifunc(x_pert)

    numerical_derivs.append((E_pert-E)/eps)

numerical_derivs = Vector1(numerical_derivs)

print ‘first 100 numerical derivs: ‘

print numerical_derivs

 

[Anonymous]

The standard centroid scorefunction (and the one you’re getting with get_cen_scorefxn()) is set up for point evaluation, not minimization. Generally speaking, most of the time that you’re using a centroid scoring function you’re doing a coarse Monte Carlo sampling and don’t need to do derivatives, so many of the centroid scoreterms are set up as a binned lookup evaluation, rather than as a continuous function. As such, there’s discontinuities and other issues which mean that the derivatives aren’t sensible.

That said, there are “smooth” versions of centroid scoreterms which are set up for minimization. Basically, the smoothed versions takes the binned statistical data the regular centroid scorefunctions are based on and do things like fit splines to them, such that you get reasonable derivative behavior. 

If you’re looking for derivatives in centroid mode, I’d recommend using something like the `score4_smooth` weights file (`create_score_function(“score4_smooth”)`). Or, if you’re interested in getting something reasonably close to what you get with get_cen_scorefxn(), use the `cen_std_smooth` weights with the `score4L` patch. (`create_score_function(“cen_std_smooth”,”score4L”)`).

[Author]

Posts