Issue when perturbing molecule about biosimspace HOT 23 CLOSED

from biosimspace.

Hi,

I think the problem might be in the type labels at least for Gromacs. The atoms that are unique to molecule one probably should be labelled as dummy atoms at lambda=1 and the atoms that are unique to to molecule two probably should be labelled as dummy atoms at lambda=0 [du for type0 or type1 in the atoms directive of morph.top].

from biosimspace.

Yes, I am doing FEP with GROMACS. Everything you said is right and I think the problem is that for atoms unique to molecule two the atom type is not changed properly to "du" in BioSimSpace. As a result, you get one atom type for two different LJ parameters, and the left one is written out in GROMACS, in this case the one with zero charge and LJ terms.

EDIT: I just realised I posted at the same time as skfegan and we seem to reach the same conclusions...

from biosimspace.

P.S. If an atom is in both molecules but is changing type, for example from hydrogen to carbon, it will also need two different atom type labels.

from biosimspace.

from biosimspace.

Just an update to say that I am working on a fix. Having looked at this and some of the example topology files included in the supplementary materials I can see what the problem is. However, it turns out to be more complicated to implement a simple fix to the existing GROMACS parser within Sire than I had anticipated. I'll post an update when I make more progress.

Thanks again for flagging this.

from biosimspace.

I've added a fix to the feature-mapping branch of Sire that hopefully resolves the issues with defining dummy atoms for perturbable molecules.

Comparing with the approach described here I still see several differences:

• We assign dummy atoms a pytpe (particle type) of D in the atomtypes section whereas the Mobley scheme sets them as A.
• We assign dummy atoms to have a value of zero for both Lennard-Jones terms, sigma and epsilon. The Mobley scheme only zeros the epsilon term.

Having not run free energy perturbation simulations using GROMACS I'm not sure of the consequences of these differences. If you get a chance to test the fix and could report back, that would be great. Now that I've fixed the basic issue (and updated the parser logic) it should be easy for me to make additional tweaks.

Cheers.

from biosimspace.

I've just run your input files through FESetup and it produces different output yet again:

• There is only a single dummy atom type in the atomtypes section labeled with a pytpe of A.
• The dummy atom has zero for both Lennard-Jones terms.
• All dummy atoms are given a mass of 1.0000 in the atoms section, rather than being assigned the mass of the real atom at the lambda end point (as is done with BioSimSpace the Mobley scheme).

I can't image that all three formats are valid. It would be good to clear this up if anyone knows better.

from biosimspace.

The mass of 1.0000 does not seem right (I have opened an FESetup issue to remind myself to fix that). I am not sure about the difference between pytpe A and pytpe D.

from biosimspace.

Hi again,

Many thanks for the fix, I will test it when I get the chance, I am currently experiencing technical problems and trying to fix them in the next few days.

As for the ptype, I think the right value should be A (if you look at the GROMACS documentation in https://ftp.gromacs.org/pub/manual/manual-2016.1.pdf, page 118, you can see that A stands for atom and D stands for virtual site).

I also think that the atom mass should be related to the perturbed atom rather than 1.000. Thanks for noting that.

from biosimspace.

Thanks. It's not 100% clear, but it seems like you only label atoms with the D or V type if you are using the virtual_sites section, which isn't used for free energy perturbation topologies. I'll update the code to use A for all types.

I'll see if I can figure out the benefits of using each approach for the LJ parameters (zero well depth only, or zero for both terms).

from biosimspace.

I just tried it with the D type and simulations fail, because GROMACS doesn't let you have a virtual site with a non-zero mass.

from biosimspace.

Okay, that's good to know. I've just pushed an update to Sire that sets all atomtypes to A. Let me know if that works. (Or just edit the topology file manually.)

Cheers.

from biosimspace.

Yes, that fixed everything. I will update later when I run the complete simulations with the new dummy atoms, but from what I am seeing, the topology files look fine.

from biosimspace.

Fantastic, thanks for the update. As I mentioned, the only wrinkle is how to handle to LJ terms but I'll have to refer to an expert for that. It should be easy enough to try both approaches and compare results.

At some point, would you be okay with sharing an example configuration file for the GROMACS free energy perturbation simulations that you are running? Currently BioSimSpace only provides support for running calculations using SOMD, but I'd like to add in GROMACS support too. (Currently you can only use BioSimSpace to do basic MD with GROMACS.)

Thanks again for raising the issue. It's good to know that someone is using BioSimSpace for its intended purpose, especially some of the features that haven't been as well tested. (Mostly due to lack of GROMACS expertise on my part.)

from biosimspace.

from biosimspace.

Thanks, @jmichel80, that's really useful. I had read somewhere about the changes in GROMACS keywords. I'll make sure to take a look at that. I can easily modify BioSimSpace to check the GROMACS version and throw an error if it's too old.

Do you have any comment on the Lennard-Jones parameters for dummy atoms in GROMACS? I had a quick read of your "reproducibility" paper and it mentions possible stability issues due to discontinuities in the potential. As mentioned above, BioSimSpace (following the SOMD pert file writer) and FESetup zero both terms (sigma and epsilon), but the Mobley alchemical-setup script and the Chodera lab's mmtools only zero epsilon.

from biosimspace.

from biosimspace.

Thanks for the clear explanation. All of the files you linked to, e.g. here, use zero for both terms, so I'll stick with that.

I'll work on getting the GROMACS driver up and running so we can re-run some of the D3R data set.

from biosimspace.

I think Julien gave a very good answer about the LJ potentials. I also think that both parameters should be zero but I am not aware of any studies that address that either.

I am working on a piece of software that attempts to fully automate high-throughput relative FEP simulations so I have written quite a bit of code, some of which closely mimicking BioSimSpace in places where I need more functionality. Part of this is a very flexible class which is designed for creating GROMACS parameter files. This class can be readily extended to AMBER and other pieces of software (I haven't done that yet, though) and I currently need the green light from my industrial sponsors to share it, but I should know where I stand in two weeks' time. After that, I hope that I could contribute some code to BioSimSpace if you think it could be useful.

Until then, here are the input parameter files I am using for my simulations (where init-lambda-state runs from 0 to the number of lambdas - 1). These have been tested on GROMACS 2018. Since there are many ways in which you can choose many of these, these files are just an example, rather than a rule.
MDPs.zip
GROMACS' website has some very good information on which name corresponds to which parameter:
http://manual.gromacs.org/documentation/2018/user-guide/mdp-options.html

from biosimspace.

Thanks for the mdp files, they'll be extremely useful. Also, definitely let us know if you'd like to contribute code, or have suggestions for missing features or improvements.

With regards to writing parameter files...

At the moment we've provided the user with a rather limited set of options in our Protocol classes which cover the basic functionality that is common to all packages we support. You create a protocol, which is then translated into appropriate configuration files and command-line arguments by a Process object, e.g. Process.Amber will convert a Protocol.Minimisation into the input files needed for a minimisation simulation run with Amber.

If you want to customise the protocol at a package specific level, e.g. if you are an Amber expert, then you can do it through the Process.Amber object, but this would break interoperability. (You now have an Amber specific script.) Also note that you also often need to know something about the system in order to write the configuration files, e.g. is the system solvated or vacuum, what box size is needed, etc.

Things get a bit complicated for more complex simulations like free energy perturbation since there are protocols for each lambda leg, then a protocol for the overall simulation. I think it would be good to refactor things a bit to aid flexibility.

We were imagining that it would be the protocols that people would customise and share, but at the moment it's a little tricky to do so

from biosimspace.

It's been a while since I started testing this system but I ran into some discrepancies which turned out to be a bug in GROMACS! Otherwise I think that the file output is correct now (although I cannot be sure). In any case, I am considering this bug as fixed and I will reopen it if I find some other issue.

With regards to the Protocols, I agree that it's difficult to implement it for many pieces of software at the same time, since different packages support different things (although I am trying to go for the highest common subset of parameters - a nontrivial task!). Since my class is only adapted for GROMACS at present, I think I need to develop it further to make it readily usable for other pieces of software. I will keep you posted.

from biosimspace.

from biosimspace.