Energy Minimization

Minimization isn’t a miracle – it will sometimes get stuck in local minima. Repacking then minimizing is often recommended, and if you really want to reduce the energy, the ramping repulsive FastRelax protocol is suggested (though more for full proteins rather than peptides). Also, if you have an isolated four residue peptide, you might not get the same per-residue energies that a full protein will (because you’re missing a number of favorable protein-protein interactions). That’s probably not your problem, though (read on).

Regarding what “dfpmin” means, see https://www.rosettacommons.org/manuals/archive/rosetta3.5_user_guide/d7/df1/minimization_overview.html It’s simply the minimization algorithm that’s being used. It’s a decent choice, though “dfpmin_armijo_nonmonotone” is likely a better overall choice (though it will depend on your particular case).

The True is a designation regarding neighbor list updating during minimization. This is an implementation detail. There are some score terms which are environmental (neighbor) dependent. Minimization can change the number of neighbors, so constantly updating the neighborlist during minimization gives more accurate energy evaluations, but at the cost of a slower runtime. Leaving it as True will probably be fine. (Looking at the code, there’s some complications, and I don’t actually think that MinMover setting actually controls much.)

10 is the tolerance, and is probably where your issues are. The default is for relative tolerance, and this is typically a fractional value. Right now you’re telling the minimizer to stop when it reaches a value which it estimates to be +- 10x of the final minimum value. E.g., if you expect a final minimum at about -6 REU, you’re telling the minimizer you’re fine with a final value anywhere from -6 REU up to -6 + -6*-10 = +54 REU. Try using something smaller, like 0.01 or 0.0001. It’ll take longer, but you should get better energies.

You can also specify absolute tolerances, by using method names like “dfpmin_atol” and “dfpmin_armijo_nonmonotone_atol”, but then you’d with a setting of 10 you’d still be looking at a large range.

The difference in energies after minimization isn’t too surprising – the different residue identities likely support the entire hairpin differently, a change that might not be captured in the CCD refinement step. If you take a look at the structures, you might see changes like the sheets of the hairpin coming apart, or other such issues. To a limited extent this can be starting structure dependent – if you repeated the CCD closure multiple times you might get different values (as I believe that CCD is stochastic.) I’d recommend trying exact replicates and see the spread of result energies. I would go with the after-full-minimization energy as indicative of the appropriateness of the mutation, but with the knowledge that minimization can cause the structure to fall apart if the starting structure is not quite right.

Minimization tolerance is typically specified as a relative value. A value of 0.01 means that the minimization algorithm thinks that the final energy value is within 1% of the true value of the local minimum. How stringent you want to set this depends on how precise you want the results. Lower values will give more minimized structures, at the expense of longer runtimes.

Regarding aligning structures, there are a number of functions in the core.scoring namespace. Take a look at calpha_superimpose_pose(modified_pose, reference_pose), which assumes that the two poses are the same length, and superimpose_pose (mod_pose, ref_pose, atom_map) which allows you to specify which atoms to superimpose over.

I wouldn’t say that when you get a structure with respect to rounds of minimization tells you anything. It’s too sensitive to subtle changes in starting conformation.

The threshold isn’t a firm value – it’s simply the convergence criteria that is used for the DFP algorithm. Regarding what to set it at, that depends on how stringently you wish to minimize. The tolerances are set up as relative values. So a tolerance of 0.01 means that the final result should be within 1% of the energy value of the local minimum. For a large system giving you total energies in the range of -300 REU or so, that means the end result will likely be within 3 REU of the true local minimum. For a smaller system with total energies around -20 REU, that means the end result will likely be within 0.2 REU of the true local minimum. If you reduce that to 0.0001, the minimizer will execute extra cycles such that the value will be 0.03 REU for the aforementioned large system, or 0.002 REU for the smaller system. Smaller numbers give more precise results, at the cost of greater runtime. It all depends on how stringent you wish your results to be.

Regarding HBond detection, it’s typically set up in Rosetta to be energy based. How the energy is computed is a complex statistics-based spline-fit polynomial. While you could theoretically change the polynomial such that the inter-atomic distances are longer, such a change would be non-trivial. Probably easier is to change the threshold energy at which one counts something as a hydrogen bond. That should allow for larger distances, but would also allow for less-ideal angles and orientations. Additionally at some point (I don’t know which) the energy would drop to zero, and no more discrimination would be possible. — If you’re finding hydrogen bonds with the HBondSet class, I believe that any (backbone) pairing with non-zero energy is included in the set already, so you would simply adjust your exclusion criteria.

I’m not quite sure how you’re detecting hydrogen bonds currently. If you’re using the HBondSet class, that should contain all the potential hydrogen bonds with non-zero energy already, so to change the threshold you would just change at what level you filter them as you iterate through the set.

If you’re looking for weak salt bridge interactions, I might suggest forgoing the internal HBond code altogether, and making a custom scanner for the interactions. You can make a double loop through all the residues looking for potential interacting pairs, and then check the relevant atom geometries. This way you can tailor your criteria to whatever specifics you’re looking for.

Bharat, can you go ahead and create a new topic for this? Thanks!

You can get the energy of a certain subset of residues with the get_sub_score(pose, residues) method of ScoreFunction. Here residues is a vector1 of bools the same length as the pose which indicates which residues to include in the scoring. (You can specify residues to exclude with get_sub_score_exclude_res() method.)

If you want specific scoreterms, probably the easiest way in PyRosetta is to simply create a new ScoreFunction with just those energy terms you want.

Unfortunately, if you just want mainchain atoms, things get a little harder. Rosetta is residue-based, so there isn’t a simple way to pull out the sum. Instead you’d probably have to put together something yourself. One easy way is to make a copy of the pose and just mutate all the residues to glycine. This would leave you with just the backbone-backbone interactions. You could try to sum up just the backbone components of the score, but that’d probably be harder to do.

It somewhat depends on the size of the complex and how much noise you expect from the modeling. As a gauge, 0.7 to 1.3 REU would be about what one would expect for a single hydrogen bond.

It’s at the level where it’s worth paying attention to, but there’s certainly no guarantee that a difference in that range would necessarily translate to an experimentally discernible difference. I’d be especially wary if it’s an atom-atom energy. The Rosetta energy function is optimized more on a whole residue level (or whole system level), rather than a per atom residue, so you may get large favorable interaction between certain atoms which are balanced by disfavorable interactions between closely associated atoms. Also, atomwise electrostatics in Rosetta aren’t as well benchmarked as other scoreterms, so caution would be advised.

One possible way to get a handle on things is to look at sets of randomly generated interactions of a similar type. Most of those won’t be good interactions, so you can get a better handle on what the noise level in your particular case is.

To apply a constraint during minimization, you’d need to add a constraint object to the pose and make sure the appropriate score term is on in the scorefunction you’re using. In your case it would be a core.scoring.constraints.AtomPairConstraint, added with the pose.add_constraint() method. You also need to make sure the atom_pair_constraint scoreterm has a non-zero weight in the scorefunction you’re using during minimization.

If you want to scan through rotamers looking for the best one, I’d look at the RotamerTrialsMover (in protocols.simple_moves). This will randomly pick a variable position in the associated packer task, then scan through all the possible rotamers in the current context, replace the current rotamer with the best one, and then move on to the next position. If you give a packer task with only a single residue allowed to repack, it will optimize that one. You can also use the RotamerTrialsMinMover to minimize each rotamer in the current context prior to scoring.

Regarding +/-30, the random() function returns a number between 0 and 1. random()*50 would thus be between 0 and 50, and -25+random()*50 would be between -25 and 25. You then add that to the current phi/psi value, and reset the angle to that value. To get a different range, you’d adjust the multiplier and the offset to get a random number in the range you want. E.g. for +/-30, you’d use -30+random()*60

If you’re looking to modify the backbone coordinates of the loop while keeping the strands fixed, you’re going to need to reorder the foldtree.

Rosetta typically uses a foldtree where the effect of backbone changes propagate to the next residue. What you want to do instead is to set up a foldtree where the two sheets are connected to each other at the root of the foldtree, and the loop is represented by two “tails” which go out to meet each other.

The way this is typically done is with loop modeling. Set up a protocols.loops.Loop object defining the loop – the constructor takes the pose-numbered start residue, end residue and cutpoint residue for the loop. Then add it to a protocols.loops.Loops object with the add_loop() method. Then you can use the protocols.loops.FoldTreeFromLoops mover. Give it the Loops object with the loops() method, and then call the apply() function on your pose. That should set up a foldtree where changes to the loop residues don’t propagate to non-loop residues.

The other caveat you need to be aware of is that changing the backbone coordinates of the loop will likely result in chainbreaks (non-natural bond lengths, angles and dihedrals at the cutpoint. You’ll want to turn on the chainbreak scoreterm in the scorefunction to tell you how bad it is, and then you can try various loop modeling techniques to attempt to close the chainbreak. (Regular minimization, CCD closure, kinematic closure, cartesian minimization, etc.)

First off, when you print the foldtree, what does it print out?

Secondly, you need to apply the reordered foldtree to the pose *before* you do your backbone manipulation. As you have it now, you change the dihedrals with the regular linear foldtree, and then impose your loop foldtree to the now-separated pose. Pose manipulation happens with the foldtree that’s currently in the pose.

Hmm, best I can tell, calpha_superimpose_pose returns the rmsd prior to superposition. Try inserting calls to the non-modifying function core::scoring::CA_rmsd(pose1,pose2) before and after superposition to check.

What do you mean by “reducing the value” for dfp_min? The dfpmin settings are for setting the type of minimizer – there’s really not any value that can be increased/reduced there.

I’m wondering if what you mean is that you’re changing the tolerance associated with the “dfpmin” minimizer. If that’s the case, I’m not quite sure what’s going on in your situation. In general, starting with the exact same pose, minimizing with a smaller tolerance shouldn’t give you higher energies – it should generally be the same or lower (though internal implementation details may make it such that there’s a very slight anomalous difference – especially if you’re using the nonmonotone variants). For different poses, though, I don’t think there’s necessarily any consistency. I wouldn’t be surprised if two very similar structures showed an energy difference at one tolerance setting that disappeared or inverted when the tolerance was lowered.

I think you may be putting too much though/effort into finding an “optimal” tolerance setting. My general recommendation is to set the tolerance such that the associated delta is in the noise level for your particular application. (e.g. for the default relative tolerance, a setting of 0.001 on a -300 REU protein will have an associated nominal delta of ~0.3 REU; if you’re using the “atol” versions of the minimizer types, then the set tolerance is the delta value.) If you’re really looking to find the lowest energy, other techniques are likely more worth your time, such as using the repulsive ramping of the relax protocol (I believe you can use a movemap to do a fixed backbone relax), doing loop remodeling, or using things like cartesian or flexible bond length/angle minimization.