Skip to content

Commit

Permalink
Merge pull request #186 from Becksteinlab/finalize-071
Browse files Browse the repository at this point in the history
finalize 0.7.1
  • Loading branch information
orbeckst authored Aug 6, 2021
2 parents 8afad4b + 8462759 commit 9c5b966
Show file tree
Hide file tree
Showing 7 changed files with 236 additions and 54 deletions.
2 changes: 1 addition & 1 deletion .readthedocs.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ formats:

# Optionally set the version of Python and requirements required to build your docs
python:
version: 2.7
version: 3.8
install:
- requirements: doc/requirements.txt
system_packages: true
5 changes: 3 additions & 2 deletions CHANGES
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@ Add summary of changes for each release. Use ISO dates. Reference
GitHub issues numbers and PR numbers.


2021-xx-xx 0.7.1
orbeckst
2021-08-04 0.7.1
orbeckst, iorga

Fixes

Expand All @@ -19,6 +19,7 @@ Fixes
error like "results.DeltaA does not contain 'Gibbs'") (PR #185)
* ensure that pickle files created under Python 2 are also readable under
Python 3
* documentation fixes (PR #176)


2021-08-02 0.7.0
Expand Down
139 changes: 139 additions & 0 deletions doc/examples/README.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
.. -*- coding: utf-8 -*-
========================
MDPOW Example: Benzene
========================

The ``examples`` directory contains input file for the solvation free
energy calculation of benzene in solvent (water and octanol), using
the OPLS-AA force field and the TIP4P water model.

There are two ways to use MDPOW:

1. commandline interface (CLI) through scripts
2. Python API


The online documentation goes into greater details. This README just
gives the briefest of summaries.



Included files
==============

The files in the examples directory are ::

benzene.yml
benzene/
├── README.txt
├── benzene.itp
├── benzene.pdb
├── benzene.pdb.png
├── dVdl.pdf
└── session.py

Run input file ``benzene.yml``
------------------------------

For running the simulations with the CLI, the run input file
``benzene.yml`` and the GROMACS input files ``benzene/benzene.itp``
and ``benzene/benzene.pdb`` are needed. The path to these files is
already in the runinput file so leave all files in place.


Benzene directory
-----------------

The ``benzene`` directory contains input files for GROMACS

* ``benzene.pdb``: coordinates in PDB format
* ``benzene.itp``: topology and parameters

an image of the chemical structure ``benzene.pdb.png``, output from
plotting dH/dlambda (for the MDPOW TI estimator) ``dVdl.pdf`` for
benzene in water.

The ``README.txt`` file explains how to try out the Python API, with
the commands recorded in ``session.py``.



Using the CLI
=============

Make sure that GROMACS commands can be run (e.g., ``source GMXRC`` or
``module load gromacs``). All commands are run from the directory that
contains the run input file ``benzene.yml``.

To calculate the water solvation free energy::

mdpow-equilibrium --solvent water benzene.yml
mdpow-fep --solvent water benzene.yml
mdpow-solvationenergy --solvent water benzene

See output in file ``energies.txt`` where all energies are in kJ/mol::

ITP mol solvent dG err_dG Coulomb err_Cou VDW err_VDW directory
---------- ---------------- ------ -------- --------- -------- --------- -------- ---------
BNZ water -8.32 1.49 +7.53 0.51 +0.78 1.40 benzene


To calculate the octanol solvation free energy::

mdpow-equilibrium --solvent octanol benzene.yml
mdpow-fep --solvent octanol benzene.yml
mdpow-solvationenergy --solvent octanol benzene
The octanol energies are appended to the ``energies.txt`` file::

BNZ water -8.32 1.49 +7.53 0.51 +0.78 1.40 benzene
BNZ octanol -17.03 1.24 +1.93 0.25 +15.10 1.22 benzene

Based on these *very* short simulations, the solvation free energy of benzene
in water of -8.32±1.49 kJ/mol is less favorable than the solvation free energy
in octanol of –17.03±1.24 kJ/mol.

Using these numbers, one can directly calculate the water-octanol partition
coefficient as

.. math::
\log P_{OW} = (\Delta G_W - \Delta G_O)/kT \log e
or use a script ::

mdpow-pow benzene/

which outputs ::

mdpow.fep : INFO [BNZ] Free energy of transfer water --> octanol: -8.712 (1.941) kJ/mol
mdpow.fep : INFO [BNZ] log P_ow: 1.517 (0.338)

and appends a line to the file ``pow.txt`` with content ::


ITP mol transferFE error logPow err_logPow directory
---------- --------- ------ ------ -------- ---------- ---------
BNZ -8.71 1.94 logPow +1.52 0.34 benzene

(It also writes output to ``energies.txt``. See online docs.)

Instead of running ``mdpow-solvationenergy`` for different solvents,
one can also just run ``mdpow-pow DIRECTORY`` and it will run the
separate free energy calculations automatically and then compute the
octanol-water partition coefficient.


The example CLI runs execute GROMACS immediately and run the FEP calculations
serially. For production calculations you will likely want to run these in
parallel on a computer cluster or HPC. The online documentation explains how to
change the runinput file to stop at various steps instead of running MD
simulations immediately.


Using the Python API
====================

See the ``benzene/README.txt`` and ``benzene/session.py`` files as well as the
online documentation.
2 changes: 1 addition & 1 deletion doc/requirements.txt
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,4 @@ pandas
pyyaml
GromacsWrapper>=0.5.1
alchemlyb
mdanalysis<2
mdanalysis
17 changes: 16 additions & 1 deletion doc/sphinx/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,8 +61,23 @@

# The short X.Y version.
version = '.'.join(packageversion.split('.')[:2])

# The full version, including alpha/beta/rc tags.
release = packageversion
##release = packageversion
#
# clean up for RTD:
# 0.7.0+5.g8afad4b.dirty --> 0.7.0+5.g8afad4b
# 0.7.0+0.g43485dd.dirty --> 0.7.0
try:
ver, rc = packageversion.split("+")
except ValueError:
ver, rc = packageversion, None

if not rc or rc.startswith('0'):
release = ver
else:
release = ver + "+" + rc.replace(".dirty", "")


# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
Expand Down
106 changes: 61 additions & 45 deletions doc/sphinx/source/scripts.txt
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,49 @@ input file (see `Equilibrium simulations`_) is the value of the
variable ``name`` in the ``[setup]`` section.


.. _runinput-file:

Run input file
--------------

The ``mdpow-*`` scripts are the **commandline interface** (CLI) to the
Python API functionality in the :mod:`mdpow` package. Whereas the
Python API requires passing of parameters to classes and functions and
therefore exposes the full flexibility of MDPOW, the CLI works with a
narrower set of options which are collected in a **run input
file**. This run input file or *RUNFILE* is named :file:`runinput.yml`
by default.

A template *RUNFILE* can be generated with
:program:`mdpow-get-runinput`

.. code-block:: bash

mdpow-get-runinput runinput.yml

which will copy the default run input file bundled with MDPOW and put
it in the current directory under the name :file:`runinput.yml`.

For an example of a *RUNFILE* see :download:`benzene.yml
<../../examples/benzene.yml>`. The comments in the file serve as the
documentation.

The run input file uses `YAML`_ syntax (and is parsed by :mod:`yaml`).

.. Note::

It is recommended to use absolute paths to file names.

.. Note::

It is recommended to enclose all strings in the input file in
quotes, especially if they can be interpreted as numbers. For
example, a name "005" would be interpreted as the number 5 unless
explicitly quoted.

.. _YAML: http://yaml.org/


.. _mdpow-equilibrium-label:

Equilibrium simulations
Expand All @@ -50,10 +93,10 @@ The :program:`mdpow-equilibrium` script
:program:`mdrun` herself (e.g. on a cluster) or let the script do it
locally on the workstation

The script runs essentially the same steps as described in the tutorial
:ref:`example_octanol-label` but it gathers all required parameters from a run
input file and it allows one to stop and continue and the protocol
transparently.
The script runs essentially the same steps as described in the
tutorial :ref:`example_octanol-label` but it gathers all required
parameters from the :ref:`run input file<runinput-file>` and it allows
one to stop and continue and the protocol transparently.

It requires as at least Gromacs 4.6.5 ready to run (check that all commands can
be found). The required **input** is
Expand All @@ -62,23 +105,6 @@ be found). The required **input** is
2. a structure file (PDB or GRO) for the compound
3. a Gromacs ITP file for the compound (OPLS/AA force field)

A template *RUNFILE* can be generated with :program:`mdpow-get-runinput`
``runinput.yml``, which will copy the default run input file bundled with
MDPOW and put it in the current directory under the name
``runinput.yml``.

For an example of a *RUNFILE* see :download:`benzene.yml
<../../examples/benzene.yml>`. It is recommended to use absolute paths
to file names. The run input file uses `YAML`_ syntax (and is parsed
by :mod:`yaml`).

.. _YAML: http://yaml.org/

.. Note:: It is recommended to enclose all strings in the input file
in quotes, especially if they can be interpreted as
numbers. For example, a name "005" would be interpreted as
the number 5 unless explicitly quoted.

The script keeps track of the stages of the simulation protocol (in
the state files :file:`water.simulation`, :file:`octanol.simulation`
etc) and allows the user to **restart from the last completed
Expand Down Expand Up @@ -141,7 +167,7 @@ input file.

You will require:

1. at least Gromacs 4.0.5 ready to run (check that all commands can
1. at least Gromacs 4.6.5 ready to run (check that all commands can
be found)
2. the results from a previous complete run of :program:`mdpow-equilibrium`

Expand Down Expand Up @@ -590,16 +616,6 @@ compounds as above.

error on **VDW**

**Vdp**

correction if the FEP are run at constant volume but the
desired solvation free energy is the Gibbs energy (**currently
neglected** and set to 0)

**errVdp**

error on **Vdp**

**w2oct**

transfer free energy from water to octanol (difference between
Expand Down Expand Up @@ -627,19 +643,19 @@ compounds as above.
.. _table-energies-label:

.. Table:: Solvation free energies for compounds in water and octanol (in kJ/mol).

========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================
molecule solvent DeltaG0 errDG0 coulomb errCoul VDW errVDW Vdp errVdp w2oct errw2oct logPOW errlogP directory
========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================
BNZ water -2.97 0.21 +7.65 0.07 -4.68 0.20 +0.00 0.00 -12.87 0.43 +2.24 0.07 benzene
BNZ octanol -15.84 0.37 +1.40 0.19 +14.44 0.32 +0.00 0.00 -12.87 0.43 +2.24 0.07 benzene
OC9 water -16.03 0.32 +27.35 0.09 -11.32 0.31 +0.00 0.00 -16.24 1.12 +2.83 0.20 octanol
OC9 octanol -32.28 1.08 +11.32 0.92 +20.96 0.56 +0.00 0.00 -16.24 1.12 +2.83 0.20 octanol
URE water -53.52 0.17 +56.94 0.11 -3.41 0.14 +0.00 0.00 +4.66 1.13 -0.81 0.20 urea
URE octanol -48.86 1.12 +35.75 1.09 +13.11 0.25 +0.00 0.00 +4.66 1.13 -0.81 0.20 urea
902 water -25.46 0.11 +34.93 0.10 -9.48 0.06 +0.00 0.00 +9.35 1.06 -1.63 0.18 water_TIP4P
902 octanol -16.11 1.05 +21.16 1.05 -5.05 0.09 +0.00 0.00 +9.35 1.06 -1.63 0.18 water_TIP4P
========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================
========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================
molecule solvent DeltaG0 errDG0 coulomb errCoul VDW errVDW w2oct errw2oct logPOW errlogP directory
========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================
BNZ water -2.97 0.21 +7.65 0.07 -4.68 0.20 -12.87 0.43 +2.24 0.07 benzene
BNZ octanol -15.84 0.37 +1.40 0.19 +14.44 0.32 -12.87 0.43 +2.24 0.07 benzene
OC9 water -16.03 0.32 +27.35 0.09 -11.32 0.31 -16.24 1.12 +2.83 0.20 octanol
OC9 octanol -32.28 1.08 +11.32 0.92 +20.96 0.56 -16.24 1.12 +2.83 0.20 octanol
URE water -53.52 0.17 +56.94 0.11 -3.41 0.14 +4.66 1.13 -0.81 0.20 urea
URE octanol -48.86 1.12 +35.75 1.09 +13.11 0.25 +4.66 1.13 -0.81 0.20 urea
902 water -25.46 0.11 +34.93 0.10 -9.48 0.06 +9.35 1.06 -1.63 0.18 water_TIP4P
902 octanol -16.11 1.05 +21.16 1.05 -5.05 0.09 +9.35 1.06 -1.63 0.18 water_TIP4P
========== ========== ======== ======== ======== ======== ======== ======== ======== ======== ======== ======== ====================


House-keeping scripts
Expand Down
19 changes: 15 additions & 4 deletions mdpow/forcefields.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,9 @@
.. autodata:: SPECIAL_WATER_COORDINATE_FILES
.. autodata:: GROMACS_WATER_MODELS
:noindex:
.. autodata:: GROMACS_SOLVENT_MODELS
:noindex:
Internal classes and functions
------------------------------
Expand All @@ -61,8 +63,9 @@
import logging
logger = logging.getLogger("mdpow.forecefields")

#: Default force field. At the moment, only OPLS-AA is directly
#: supported.
#: Default force field. At the moment, OPLS-AA, CHARMM/CGENFF, and AMBER/GAFF
#: are directly supported. However, it is not recommended to change the
#: default here as this behavior is not tested.
DEFAULT_FORCEFIELD = "OPLS-AA"

#------------------------------------------------------------
Expand Down Expand Up @@ -127,8 +130,16 @@ def _create_water_models(watermodelsdat):
description=description)
return models

#: Use the default water model unless another water model is chosen in the runinput
#: file (``setup.watermodel``).
#: Use the default water model unless another water model is chosen in the
#: :ref:`run input file <runinput-file>` file by setting the
#: ``setup.watermodel`` parameter.
#:
#: .. warning::
#: Select the native water model **manually** and do not rely on the default
#: set here. For CHARMM/GAFF the CHARMM TIP3P model is recommended.
#: For AMBER/GAFF the TIP3P mdeol is often used. Choosing the correct water model
#: is a *scientific* decision that *you* have to make conscientiously.
#:
DEFAULT_WATER_MODEL = "tip4p"

#: Dictionary of :class:`GromacsSolventModel` instances, one for each Gromacs water
Expand Down

0 comments on commit 9c5b966

Please sign in to comment.