Support toluene solvent (#224)

* toluene solvent parameters for OPLS-AA and GAFF * add toluene as solvent (#223) * add toluene option for AMBER and OPLS-AA in forcefields.py * add Gtol and pTW for the solvation free energy of a solute in toluene. * Add TolueneSimulation for running equilibrium simulation for a solute in toluene. * add toluene as options for equilibrium and fep simulations * "[atomtypes]" in 1tol.itp are moved to ffnonbonded.itp and renamed to avoid conflict with atomtypes in solute parameters generated by GAFF * add tests for OPLS-AA and AMBER toluene solvent * add default box size for toluene solvent * move new gaff solvent parameters to a separate itp file * change 1tol_gaff to 1tol_amber * add script to calculate Ptw * add test for toluene fep analysis function * support gaff2 atomtypes * add descriptions of the parameters * update docs Co-authored-by: Bogdan Iorga <[email protected]> Co-authored-by: Oliver Beckstein <[email protected]>
Becksteinlab · Jan 11, 2023 · 9849104 · 9849104
1 parent e395ac7
commit 9849104
Show file tree

Hide file tree

Showing 16 changed files with 16,048 additions and 11 deletions.
diff --git a/mdpow/equil.py b/mdpow/equil.py
@@ -56,7 +56,7 @@
 # TODO: change to water distance 1.2 in the future (1.0 for
 #       compatibility with our SAMPL5 runs)
 #: minimum distance between solute and box surface (in nm)
-DIST = {'water': 1.0, 'octanol': 1.5, 'cyclohexane': 1.5, 'wetoctanol': 1.5}
+DIST = {'water': 1.0, 'octanol': 1.5, 'cyclohexane': 1.5, 'wetoctanol': 1.5, 'toluene': 1.5}
 
 class Simulation(Journalled):
     """Simple MD simulation of a single compound molecule in water.
@@ -123,7 +123,7 @@ def __init__(self, molecule=None, **kwargs):
           *forcefield*
               'OPLS-AA' or 'CHARMM' or 'AMBER'
           *solvent*
-              'water' or 'octanol' or 'cyclohexane' or 'wetoctanol'
+              'water' or 'octanol' or 'cyclohexane' or 'wetoctanol' or 'toluene'
           *solventmodel*
               ``None`` chooses the default (e.g, :data:`mdpow.forcefields.DEFAULT_WATER_MODEL`
               for ``solvent == "water"``. Other options are the models defined in
@@ -658,3 +658,8 @@ def  _setup_solvate(self, **kwargs):
         ionkwargs['struct'] = sol['struct']
         params = gromacs.setup.solvate_ion(**ionkwargs)
         return params
+
+class TolueneSimulation(Simulation):
+    """Equilibrium MD of a solute in a box of toluene."""
+    solvent_default = 'toluene'
+    dirname_default = os.path.join(Simulation.topdir_default, solvent_default)
diff --git a/mdpow/fep.py b/mdpow/fep.py
@@ -87,10 +87,15 @@
    :members:
    :inherited-members:
 
+.. autoclass:: Gtol
+   :members:
+   :inherited-members:
+
 .. autofunction:: pOW
 
 .. autofunction:: pCW
 
+.. autofunction:: pTW
 
 Developer notes
 ---------------
@@ -1325,6 +1330,13 @@ class Gwoct(Goct):
     dirname_default = os.path.join(Gsolv.topdir_default, solvent_default)
 
 
+class Gtol(Gsolv):
+    """Sets up and analyses MD to obtain the solvation free energy of a solute in toluene.
+    """
+    solvent_default = "toluene"
+    dirname_default = os.path.join(Gsolv.topdir_default, solvent_default)
+
+
 def p_transfer(G1, G2, **kwargs):
     """Compute partition coefficient from two :class:`Gsolv` objects.
 
@@ -1440,6 +1452,7 @@ def p_transfer(G1, G2, **kwargs):
     # lower case initials, in reverse order of transfer, e.g.
     # water -> octanol:      P_ow
     # water -> cyclohexane:  P_cw
+    # water -> toluene:      P_tw
     coefficient = "P_{0}{1}".format(
         G2.solvent_type.lower()[0], G1.solvent_type.lower()[0])
 
@@ -1542,3 +1555,48 @@ def pCW(G1, G2, **kwargs):
         logger.error(msg)
         raise ValueError(msg)
     return p_transfer(*args, **kwargs)
+
+def pTW(G1, G2, **kwargs):
+    """Compute water-toluene partition coefficient from two :class:`Gsolv` objects.
+
+    transfer free energy from water into toluene::
+
+            DeltaDeltaG0 = DeltaG0_tol - DeltaG0_water
+
+    toluene/water partition coefficient::
+
+            log P_tol/wat =  log [X]_tol/[X]_wat
+
+    :Arguments:
+       *G1*, *G2*
+           *G1* and *G2* should be a :class:`Ghyd` and a :class:`Gtol` instance,
+           but order does not matter
+       *force*
+           force rereading of data files even if some data were already stored [False]
+       *stride*
+           analyze every *stride*-th datapoint in the dV/dlambda files
+       *start*
+           Start frame of data analyzed in every fep window.
+       *stop*
+           Stop frame of data analyzed in every fep window.
+       *SI*
+           Set to ``True`` if you want to perform statistical inefficiency
+           to preprocess the data.
+       *estimator*
+           Set to ``alchemlyb`` if you want to use alchemlyb estimators,
+           or ``mdpow`` if you want the default TI method.
+       *method*
+           Use `TI`, `BAR` or `MBAR` method in `alchemlyb`, or `TI` in `mdpow`.
+
+    :Returns: (transfer free energy, log10 of the toluene/water partition coefficient = log Ptw)
+    """
+    if G1.solvent_type == "water" and G2.solvent_type == "toluene":
+        args = (G1, G2)
+    elif G1.solvent_type == "toluene" and G2.solvent_type == "water":
+        args = (G2, G1)
+    else:
+        msg = "For pTW need water and toluene simulations but instead got {0} and {1}".format(
+            G1.solvent_type, G2.solvent_type)
+        logger.error(msg)
+        raise ValueError(msg)
+    return p_transfer(*args, **kwargs)
diff --git a/mdpow/forcefields.py b/mdpow/forcefields.py
@@ -173,7 +173,9 @@ def get_water_model(watermodel=DEFAULT_WATER_MODEL):
         identifier="wetoctanol", itp="1octwet.itp", coordinates="1octwet.gro"),
     'wetoctanolnew': GromacsSolventModel(
         identifier="wetoctanol", itp="1octwetnew.itp", coordinates="1octwet.gro",
-        description=new_octanol)
+        description=new_octanol),
+    'toluene': GromacsSolventModel(
+        identifier="toluene", itp="1tol.itp", coordinates="1tol_oplsaa.gro"),
     }
 
 CHARMM_SOLVENT_MODELS = {
@@ -192,6 +194,8 @@ def get_water_model(watermodel=DEFAULT_WATER_MODEL):
         identifier="wetoctanol", itp="1octwet.itp", coordinates="1octwet_amber.gro"),
     'cyclohexane': GromacsSolventModel(
         identifier="cyclohexane", itp="1cyclo.itp", coordinates="1cyclo_amber.gro"),
+    'toluene': GromacsSolventModel(
+        identifier="toluene", itp="1tol.itp", coordinates="1tol_amber.gro"),
     }
 
 #: Solvents available in GROMACS; the keys of the dictionary
@@ -282,7 +286,8 @@ def get_top_template(identifier):
     """Return the topology file template suitable for the solvent model."""
 
     templates = {'water': 'system.top', 'octanol': 'system.top',
-                 'cyclohexane': 'system.top', 'wetoctanol': 'system_octwet.top'}
+                 'cyclohexane': 'system.top', 'wetoctanol': 'system_octwet.top',
+                 'toluene': 'system.top',}
     try:
         return templates[identifier]
     except KeyError:

diff --git a/mdpow/run.py b/mdpow/run.py
@@ -166,6 +166,7 @@ def equilibrium_simulation(cfg, solvent, **kwargs):
         'octanol': equil.OctanolSimulation,
         'wetoctanol': equil.WetOctanolSimulation,
         'cyclohexane':equil.CyclohexaneSimulation,
+        'toluene': equil.TolueneSimulation,
         }
     try:
         Simulation = Simulations[solvent]
@@ -277,19 +278,21 @@ def fep_simulation(cfg, solvent, **kwargs):
         'water': equil.WaterSimulation,
         'octanol': equil.OctanolSimulation,
         'wetoctanol': equil.WetOctanolSimulation,
-        'cyclohexane':equil.CyclohexaneSimulation
+        'cyclohexane': equil.CyclohexaneSimulation,
+        'toluene': equil.TolueneSimulation,
         }
     Simulations = {
         'water': fep.Ghyd,
         'octanol': fep.Goct,
         'wetoctanol': fep.Gwoct,
-        'cyclohexane':fep.Gcyclo
+        'cyclohexane':fep.Gcyclo,
+        'toluene': fep.Gtol,
         }
     try:
         EquilSimulation = EquilSimulations[solvent]
         Simulation = Simulations[solvent]
     except KeyError:
-        raise ValueError("solvent must be 'water', 'octanol', 'wetoctanol' or 'cyclohexane'")
+        raise ValueError("solvent must be 'water', 'octanol', 'wetoctanol', 'cyclohexane' or 'toluene'")
     # generate a canonical path under dirname
     topdir = kwargs.get("dirname", None)
     if topdir is None:

diff --git a/mdpow/tests/test_fep_analysis.py b/mdpow/tests/test_fep_analysis.py
@@ -76,3 +76,7 @@ def test_pOW_error(Ghyd, Ghyd_other):
 def test_pCW_error(Ghyd, Ghyd_other):
     with pytest.raises(ValueError):
         mdpow.fep.pCW(Ghyd, Ghyd_other)
+
+def test_pTW_error(Ghyd, Ghyd_other):
+    with pytest.raises(ValueError):
+        mdpow.fep.pTW(Ghyd, Ghyd_other)
diff --git a/mdpow/tests/test_forcefields.py b/mdpow/tests/test_forcefields.py
@@ -7,7 +7,7 @@
 
 # currently supported
 WATERMODELS = ('tip4p', 'tip3p', 'tip5p', 'spc', 'spce', 'm24', 'tip4pd')
-SOLVENTMODELS = ('water', 'cyclohexane', 'octanol')
+SOLVENTMODELS = ('water', 'cyclohexane', 'octanol', 'toluene')
 
 class TestIncludedForcefiels(object):
     @staticmethod
@@ -40,9 +40,13 @@ class TestIncludedSolvents(object):
             '1cyclo.gro': os.path.join('mdpow', 'top', '1cyclo.gro'),
             '1cyclo.itp': os.path.join('mdpow', 'top', 'oplsaa.ff', '1cyclo.itp')
         },
+        'toluene': {
+            '1tol_oplsaa.gro': os.path.join('mdpow', 'top', '1tol_oplsaa.gro'),
+            '1tol.itp': os.path.join('mdpow', 'top', 'oplsaa.ff', '1tol.itp')
+        },
     }
 
-    @pytest.mark.parametrize("solvent_name", ["tip4p", "octanol", "cyclohexane"])
+    @pytest.mark.parametrize("solvent_name", ["tip4p", "octanol", "cyclohexane", "toluene"])
     def test_solvent(self, solvent_name):
         solvent = self.solvents[solvent_name]
         for filename, path in solvent.items():
@@ -123,6 +127,12 @@ def test_get_solvent_wetoctanol(self, forcefield):
         assert (mdpow.forcefields.get_solvent_model(model, forcefield=forcefield) is
                 mdpow.forcefields.GROMACS_SOLVENT_MODELS[forcefield][model])
 
+    @pytest.mark.parametrize("forcefield", ['OPLS-AA', 'AMBER'])
+    def test_get_solvent_toluene(self, forcefield):
+        model = 'toluene'
+        assert (mdpow.forcefields.get_solvent_model(model, forcefield=forcefield) is
+                mdpow.forcefields.GROMACS_SOLVENT_MODELS[forcefield][model])
+
     @staticmethod
     def test_get_solvent_identifier_default_is_water():
         assert (mdpow.forcefields.get_solvent_identifier('water') is

diff --git a/mdpow/tests/test_solvation.py b/mdpow/tests/test_solvation.py
@@ -14,6 +14,7 @@
         "octanol" : equil.OctanolSimulation,
         "cyclohexane" : equil.CyclohexaneSimulation,
         "wetoctanol" : equil.WetOctanolSimulation,
+        "toluene": equil.TolueneSimulation,
         }
 
 test_file = {"OPLS-AA": 'benzene.itp',
@@ -52,6 +53,10 @@ def test_solvation_octanol(setup, ff):
 def test_solvation_cyclohexane(setup):
     solvation(setup, "cyclohexane")
 
+@pytest.mark.parametrize("ff", ['OPLS-AA', 'AMBER'])
+def test_solvation_toluene(setup, ff):
+    solvation(setup, "toluene", ff)
+
 @pytest.mark.xfail(gromacs.release.startswith('4')
                    or gromacs.release.startswith('5')
                    or gromacs.release.startswith('2016'),