I am working with a protein with a bound small molecule (non-covalent). I've managed in creating parameters for the ligand and relaxing the complex following this workflow: https://new.rosettacommons.org/demos/latest/public/relax_around_chemically_bound_ligand/README
Next I'd like to do the same but including some ligand flexibility. I've generated conformers and new parameters, but when I try to relax the complex, the ligand coordinates are wrong. Not sure if there is something that I should do differently.
For context, I'm allowing a limited set of conformational flexibility to my ligand - just sample one small substituent exposed to solvent and generate aprox. 5 conformers. I would like my ligand core to stay in place and sample different conformations for the substituent. Is this possible? Is there an example of a workflow like this?
One thing I didn't see mentioned in that tutorial is checking the neighbor atom.
When substituting conformers, it's the neighbor atom which is kept fixed. (The confomers are all aligned on the neighbor atom and its environment.) Normally, the neighbor atom is chosen to be in the center of the ligand. This minimizes changes when the ligand is free, but if the ligand has a covalent attachement that's separated from the neighbor atom by a rotatable bond, this means that the connection point can move in 3D space when you get conformer substitution, leading to larger distances in the connecting bond.
What you would want to do is edit the params file to change the NBR_ATOM line to specify the connection point. (Or ideally, an atom with at least two heavy atom connections which is ridgidly attached to the connection point.) This way when the rotamers are substituted, the position and orientation of the connection point stays fixed, so you don't open up the bond of the ligand to the rest of the system.
Note that if you change the NBR_ATOM you would also want to change the NBR_RADIUS line. The NBR_RADIUS helps speed calculations, and is the furthest distance from the NBR_ATOM to any other heavy atom in any potential rotamer. (You can get inaccurate energies if this is too short, but there's little penalty beside being a bit slower if you go too long.)
Thanks for the quick response! I think I found the problem - my input structure was the same as my original single conformer test, and it looks like the new parameters are not compatible (different atom order?) leading to garbled structures. Works OK after rebuilding the input structure.
Quick follow-up question - how are ligand conformations explored in relaxation ? I am using the options below, I'm assuming this should constraint only the backbone, but I get very low diversity of ligand conformations. I don't know if most of them are found unfavorable for my protein conformation, or if I need to add something else to sample ligand conformations.
-nstruct 100
-relax:constrain_relax_to_start_coords
-relax:ramp_constraints false
-ex1
-ex2
-use_input_sc
-flip_HNQ
-no_optH false
-in:file:extra_res_fa RL.fa.params
-in:file:extra_res_fa TL.fa.params
All ligand conformers in the conformer file are considered. The first things that happens is that each conformer is realigned against the current conformation by aligning the environment around the neighbor atom. There's some additional diversification, mostly with hydroxyl rotamers. Then there's a rudimentary bump check, to make sure that atoms in the conformers aren't clashing with backbones or other atoms which can't move out of the way. Then packing occurs, where the conformers are sampled in a Monte Carlo process. The settings are such that you can assume that all ligand conformers which have made it to that point are tested to see if they help make a low energy structuture. All but one of the conformers are discarded by the end of this process, and only the one conformer which is the low energy one (or as close as a random Monte Carlo process can get it) will be taken forward.
Note that this process of generating, aligning and filtering happens anew each time the packer is called, so that the conformer assement is redone frequently during the relax process, every time the backbone changes.
Generally, if you don't have setting which are turning off ligand repacking entirely, all of the rotamers are considered, and only the lowest energy one shows up in your results. Depending on your starting point, relax is decently convergent, so it's not unexpected that you get more-or-less the same ligand conformation each time, especially if you're constraining the backbone.
The backbone constraints (-relax:constrain_relax_to_start_coords & -relax:ramp_constraints false) will only affect protein residues (residues with backbones), and won't affect non-polymeric ligands (which are all sidechain for Rosetta). But the context in which they're placed is rather constrained, so energy minimization considerations mean that it's not unexpected that they'll find more or less the same location, particularly if that location is well matched to the ligand structure.