Comments (23)
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) ofD
in theatomtypes
section whereas the Mobley scheme sets them asA
. - 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 apytpe
ofA
. - 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.
Related Issues (20)
- Storage quota exceeded on the Anaconda cloud HOT 1
- HILLS file not updating during restarting simulation HOT 4
- Is it possible to add ions to an already solvated system, without adding any more water molecules? HOT 4
- Problem when converting Gro/Top to Rst7 HOT 12
- Making HMR work with dummies HOT 11
- Make _removeDummies() recognise dummy bonds HOT 10
- Would it be possible to include HiMap as an alternative to LoMap in the generateNetwork method? HOT 2
- Unable to import BioSimSpace after mamba installation HOT 12
- BioSimSpace.Align.viewMapping: py3Dmol.view not rendering in Jupyter Notebook HOT 11
- `_toRegularMolecule` failing with `convert_amber_dummies=True` HOT 15
- Absolute Binding Free Energy Calculations with GROMACS HOT 7
- Pickling a parametrised molecule from an SDF file with GAFF2 HOT 4
- Are atomtypes meant to be updated when using repartitionHydrogenMass? HOT 16
- [BUG] OSError "It looks like you failed to include a topology file." when reading SDF which doesn't contain redundant bond order information HOT 10
- [BUG] Cannot import BioSimSpace after installation of biosimspace 2023.1.2 using mamba HOT 3
- AnalysisError: SOMD free-energy analysis failed! HOT 5
- [BUG] HOT 10
- Reporting a vulnerability
- [BUG] Water molecule coordinates are setting to (0,0,0) HOT 4
- Issue in reproducing the funnel methadynamics tutorial HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from biosimspace.