michellab / BioSimSpace

Code and resources for the EPSRC BioSimSpace project.
https://biosimspace.org
GNU General Public License v3.0
77 stars 19 forks source link

SOMD #28

Closed lohedges closed 5 years ago

lohedges commented 6 years ago

I'm currently working on the SOMD molecular dynamics driver that will be part of the BioSimSpace.Process package. (This is separate to the SOMD free energy driver.) I've run into a few issues that I'll document here.

1) No access to thermodynamic time series data, e.g. energy, temperature, pressure, etc. As such, the only functionality of the Process.Somd class is getTrajectory and getSystem (which calls getTrajectory and grabs the last frame). This makes SOMD of limited use for interactive simulation, i.e. no ability to extract and monitor data as the simulation runs. (Perhaps data can be computed from within Sire using some of the restart files?) 2) You can only run simulations at fixed temperature, so no heating/cooling protocols are possible. 3) The OpenMMD.py script that is part of SOMD uses Sire.Units when setting default parameters, e.g. the integration time step. However, no units are used when reading parameters from a configuration file so the parser will have to be modified in order to support any user specified values that come from BioSimSpace.Protocol objects. At present, the parser just reads the configuration values as strings, then blindly passes them to any function in which they are needed (which then complain that the arguments are of the wrong type).

lohedges commented 6 years ago

Actually, it turns out that you need to use the exact Sire.Units names in the configuration file, for example:

equilibration timestep = 0.0005 picosecond

Not...

equilibration timestep = 0.0005 ps

as is generated by running

$HOME/sire.app/bin/somd --help-config

SOMD will fail if the unit name is incorrect, or when it is missing, i.e. no default unit set if the variable is just an integer or float. This should be okay for the purposes of driving SOMD via BioSimSpace, since the variables will be type checked and validated in the protocol objects and we have full control over writing the configuration file. However, this isn't particularly useful for someone wanting to run SOMD directly. I suggest updating the --help-config output to display the correct units, or implementing some kind of unit validation within the resolveParameters decorator from here.

jmichel80 commented 6 years ago

@ppxasjsm can you raise an issue on the sire git page to make sure this gets fixed ?

lohedges commented 6 years ago

I'm just doing this now :-)

lohedges commented 6 years ago

Another question for the SOMD experts: How do you implement solvent constraints, as used in, e.g. SHAKE? For the other MD drivers we default to a 2 femtosecond time step and make solvent bonds rigid. From looking at the output of $HOME/sire.app/bin/somd --help-config it looks like it's possible to use a harmonic restraint on the solvent atoms, e.g.

use restraints = True
restraint force constant = 100.0

Is this what I should be using?

jmichel80 commented 6 years ago

Hi Lester,

Bonds in water molecules are constrained by default unless you set constraint = none.

To avoid bugs, water should have residue name 'WAT' in the topology. This is especially important for water models that have extra point site (e.g. TIP4P).

best wishes,

Julien


Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Tue, Oct 9, 2018 at 11:22 AM Lester Hedges notifications@github.com<mailto:notifications@github.com> wrote:

Another question for the SOMD experts: How do you implement solvent constraints, as used in, e.g. SHAKE? For the other MD drivers we default to a 2 femtosecond time step and make solvent bonds rigid. From looking at the output of $HOME/sire.app/bin/somd --help-config it looks like it's possible to use a harmonic restraint on the solvent atoms, e.g.

use restraints = True restraint force constant = 100.0

Is this what I should be using?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/28#issuecomment-428140201, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALLd5L4dxuoOQJABy1eM_7znFpLg42IOks5ujHjXgaJpZM4XCW_w.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

lohedges commented 6 years ago

Thanks, that's good to know. I'll add some functionality to detect the water molecules and change the residue name accordingly during the setup stage of any SOMD simulation. For systems that are solvated by BioSimSpace, the default water name is SOL, since this is what is used by gmx solvate.

lohedges commented 6 years ago

From the --help-config output it looks like the option for enabling a (Andersen) thermostat is:

thermostat = True

However, looking at some example SOMD input files, e.g. here, I see the option:

andersen = True

Are these two options equivalent, or has one superseded the other?

jmichel80 commented 6 years ago

Hi Lester,

This looks like a bug in the cfg file at the link you pulled, one should write thermostat = True what happens is that thermostat = True is the default setting. setting andersen = True triggers no code. One annoying problem with the parser is that keywords that have no meaning are silently ignored. This can cause confusion when making typos writing config files.

Best wishes,

Julien


Dr. Julien Michel, Senior Lecturer Room 263, School of Chemistry University of Edinburgh David Brewster road Edinburgh, EH9 3FJ United Kingdom phone: +44 (0)131 650 4797 http://www.julienmichel.net/

On Thu, Oct 11, 2018 at 10:41 AM Lester Hedges notifications@github.com<mailto:notifications@github.com> wrote:

From the --help-config output it looks like the option for enabling a (Andersen) thermostat is:

thermostat = True

However, looking at some example SOMD input files, e.g. herehttps://github.com/michellab/cyp-trivector/tree/master/SOMD_input/Br-aryl_series/CypA/22_BM1%7EI05_BM1/bound/run001/sim/input, I see the option:

andersen = True

Are these two options equivalent, or has one superseded the other?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/michellab/BioSimSpace/issues/28#issuecomment-428890413, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ALLd5D6PHadk-AYcxHeyRk9X90edblg5ks5ujxI0gaJpZM4XCW_w.

The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.