I am recording the positions of all atoms before I run a Rosetta protocol + external software perturbations, and using these positions to "reset" the system after the protocals have completed. The score is computed before and after the protocols, and after resetting the system. I noticed resetting causes a large increase in the score (-1000 to 25000).
Looking at the score by residue, the culprit was pro_close for Prolines...
-My first guess is that this is an improper ring angle, except that this structure was a starting structure in Rosetta.
-I am fairly sure all atoms are in place, as I have visualized the structure and compared atom lists in outputted PDBs.
-The pro_close term is the only large term in the score.
Are there some other structural changes of residues that occur with repacking during a Rosetta protocol? (e.g. HIS_D to HIS_E)
It's a known issue that straight-from-the-PDB structures tend to score abysmally poor under the Rosetta score function. The issue is that the Rosetta energy function is sensitive enough that even small (0.1 Ang) movements can result in large scorefunction values. So small differences in atom position which wouldn't make a difference for electron density do make a difference for Rosetta.
Usually it is recommended to relax the structure with constraints to the starting atom positions (http://www.rosettacommons.org/manuals/archive/rosetta3.2.1_user_guide/pr... in the "regular" Rosetta manual) - I'm not sure how best to do it through PyRosetta, though. You could play around with the FastRelax mover (protocols::relax::FastRelax)
A) Prolines have a virtual atom. I forget its name but I think it's NV. Anyway, your coordinate reset is probably not resetting the virtual atom, so it remains out of place. The only thing the virtual atom is used for is to check its overlap with the proper backbone atom (probably N), so the only score that goes haywire is pro_close. Prolines aren't rings in the atom_tree because the atom_tree is a directed acyclic graph (no rings!)
B) If you are resetting entire poses, why not just use an equality operator? Is there not one for poses? Perhaps python problematically forces shallow copies when you want deep copies?
Once again, Steven, your advice saved me hours of digging. Thank you.
For anyone interested, I found a quick and semi-clean solution. Rosetta.src.core.Conformation.cc:
//////////////////////////////////////////////////////////////////////////////////@details
/// @note This could be rewritten to avoid a refold using set_bond_angle, etc might be safer since
/// by-xyz ignores possibility of propagating dependencies from the moving atoms...
void
Conformation::rebuild_residue_connection_dependent_atoms( Size const seqpos, Size const connid )
{
// assumes no interpendencies among the atoms being built
update_residue_coordinates();
Residue const & rsd( residue_(seqpos) );
for ( Size i=1, ie=rsd.natoms(); i<= ie; ++i ) {
if ( rsd.icoor(i).depends_on_residue_connection( connid ) ) {
set_xyz( AtomID(i,seqpos), rsd.icoor(i).build( rsd, *this ) );
}
}
}
From this you can pull out the "set_xyz" line, which is creating an AtomID for whatever atom in the residue you'd like and then using the "build(residue,conformation())" method to get its ideal coordinate. To move a virtual atom, simply create an AtomID for it and grab a Residue.icoor(idx).build(residue,conformation()) vector. Both of these go into the pose.conformation().set_xyz(AtomID,xyzVector) method.
In my case I needed to move the virtual atom 'NV' (found in minirosetta_database/chemical/residue_type_sets/fa_standard/residue_types/l-caa/AMINOACID.param) on a Proline, after I had moved the rest of the atoms in the residue.