ialibay / mdrestraintsgenerator Goto Github PK

So ideally we would want to not have to have a user defined ligand selection.

Ideas for ligand picking strategies could be:

1. Calculate an RMSF on all ligand atoms.
   a) We could take two heavy atoms with the lowest RMSF and one with the highest. This should help reduce the flexibility of the ligand.

2. Use ligand sub-graphs to pick out ring atoms.
   a) At the moment the human side of things usually picks out ring one ring atom (at least for fragments) and it's next two connected heavy atom. We could probably just use sub-graphs or gaff atom naming to pick those out automatically.

As per the title.

Doing this should allow for a cleaner API + custom ligand/protein picking.

The circmean & circvar methods have been changed in scipy, with the latter providing a normalised variance. Technically these now better reflect the way in which you calculate circular statistics, however this changes the weighting we had towards the combined variance (which in some ways highlights an issue with the maybe naive way in which we did things). This needs further investigation but I believe that the bond variance now has a greater contribution towards the summed variance.

For now we will pin to scipy<1.8, and then investigate how to move from here.

As evidenced by the recent gmx redmine issue: https://redmine.gromacs.org/issues/2953, it may be prudent to implement a means of generating Boresch restraints that does not rely on the intermolecular_restraints code. This would allow us to leverage the pull code and help bring the .mdp files closer to those of the hardwall restraint simulations. It would also probably make the restraints easier to understand for new users.

This should be rather simple:

1. Generate a set of index groups using MDAnalysis which defines each of the 6 atoms involved in a Boresch restraint.
2. Generate a set of default .mdp options for the restraints.. something like:
   pull-group1-name = atom1
   pull-group2-name = atom2
   ... etc...
   pull-coord1-type = umbrella
   pull-coord2-type = umbrella
   .... etc ...
   pull-coord1-geometry = distance
   pull-coord2-geometry = angle
   ... for pull-coord3
   pull-coord4-geometry = dihedral
   ... for pull-coord[5-6]
   pull-coord1-init = [nm]
   pull-coord2-init = [deg]
   ... etc...

Ideally we can just provide this as a printed output for pasting into future mdp files and/or we can just create idealised template mdp files directly using this protocol.

See #21

I can't seem to be able to reproduce this locally though, and MacOS runs fine.. similarly to the MDA user guide's issues, this one's going to take time to solve..

This was brought up thanks to @xiki-tempula

The code expects specific topology attributes that are only contained in more detailed topology formats (tpr, parm7, etc...), but does not make it clear that these are necessary & that these topology attributes are missing. Proper checks w/ error messages should be added.

As per the title.

As per title, we need tests and examples to make sure this is reproducible.

It works but it's horribly inefficient. We can just use one of the intersection class methods for this.

Due to the way gromacs calculations dihedrals we end up with issues when crossing from 180 to -180 on a dihedral. To avoid this we should just make sure that angles that are within a given cutoff of that boundary are given really high variances (so that they won't get picked).

To do - this is documented in one of Vytas' recent papers on NEQ absolute binding free energy calculations.

At the moment the histograms plot the mean value for each distribution, ideally we would also want to show the values of the closest frame, this would ensure that we are picking something sensible.

The code update for this will require implementing the rmsd calculation as part of the class object, but it shouldn't cause large increases in computational costs (unless many frames are involved).

Probably needs a full refactor, we should be allowing folks to analyse whatever portion of the trajectory they are looking to run.

On top of that, we probably should avoid doing a transfer to memory where we can avoid it - this is horribly memory intensive.

This is a lot more efficient and ensures ordering.

As a means to test the AMBER GPU-TI code for Boresch restrained TI simulations, it would be useful to use the code to dump out an AMBER topology that has fictitious bonds added via parmed (see the setBond, setAngle and addDihedral + timerge commands). This should allow us to do a set of transformations where we:

A) Transform the system into a bonded one.
B) Do a one or two step charge/vdw perturbation.

Some of the issues here is that the MBAR code may not work properly, but it is worth trying just in case (if anything it'd make a nice tutorial).

Unfortunately these kinds of scenarios can happen...

This is problematic because the lowest RMSD position ends up picking up a rather unfavourable dihedral angle (I'll update if this simulation doesn't fail...).

For cases such as these, it would be useful to implement an interactive mode where the users can limit the sampling selection to a region of interest.

This is going to be rather large thing to implement, it'll probably a need a complete refactoring of the code.

One way to fix this in the short term is to look at the range in the values and see if they exceed 180 degrees, then shift all values accordingly.

This might be an issue for the mean, so it will have to be treated accordingly.

Minimal usage documentation should be added temporarily so that people at least know where to look for things.

Already partially started, the plot() and write() methods should allow users to define a path.

As per title, in some cases resid+1 doesn't exist or even worse, in some cases you end up chain hopping.

Either way this is generally a terrible idea to try to get a restraint for, but hopefully it should get filtered out by the restraint picker.

From @jbluck
It would be useful to have this code automatically write out the analytical restraint energy.

Right now the plots are unexplained and need legends.

Somehow I managed to make it so that the Angle and Dihedral class objects all have the same filename for plotting, so they just overwrite each other.