Solvating an MDMC Universe

For simulating systems in solution, MDMC’s inbuilt solvate method can be used to add solvent molecules of a desired density to your universe.

Create the Universe

See the Building a Universe wiki for details on how to create a Universe of your own specifications. For the purposes of this tutorial, we will be solvating a simple universe that contains only 4 Hydrogen molecules:

# Import the Atom, Molecule, and Universe classes
# Import the HarmonicPotential class (needed to create a Bond)
from MDMC.MD import Atom, Bond, HarmonicPotential, Molecule, Universe

# Initialise a Universe with dimensions in Ang
universe = Universe([10.0, 15.0, 20.0])

# Create a pair of Hydrogen Atoms
H1 = Atom('H')
H2 = H1.copy(position=[1., 1., 1.])

# Initialise a H-H Bond
HH_bond = Bond(H1, H2, function=HarmonicPotential(1., 100., interaction_type='bond'))

# Make a H2 Molecule
H2_mol_1 = Molecule(atoms=[H1, H2], interactions=[HH_bond])

# Create 3 copies of the Molecule at different positions
H2_mol_2 = H2_mol_1.copy(position=[3.5, 3.5, 3.5])
H2_mol_3 = H2_mol_1.copy(position=[6.5, 6.5, 6.5])
H2_mol_4 = H2_mol_1.copy(position=[9.5, 9.5, 9.5])

# Add the 4 Hydrogen Molecules to the Universe
for molecule in [H2_mol_1, H2_mol_2, H2_mol_3, H2_mol_4]:
    universe.add_structure(molecule)

Universe created with:
Dimensions [10. 15. 20.]

Solvating the Universe

The solvate method accepts 3 main parameters; density, tolerance, and solvent (as well as some **settings). An example call to solvate the Universe created above would therefore be:

universe.solvate(0.6, tolerance=1, solvent='SPCE')

Force field created by solvent SPCE

Parameters

Explanations of each parameter passed to solvate can be found below.

1. density

The desired density of the bulk solvent in your universe, in MDMC units of amu / Ang ^ 3. See Units in MDMC for instructions on how to convert your density into MDMC units.

Note: with this parameter you are specifying the bulk density of the solvent. If you have any solute molecules already present in your universe, the total density of your universe after solvation will be higher than the desired density you pass to solvate (plus or minus the tolerance you pass, see below).

In the above example, passing the density as 0.6 means that the universe will be solvated with SPCE water with a bulk density of 0.6 amu / Ang ^ 3 (+/- the tolerance). The density of the universe in total will be greater than 0.6 amu / Ang ^ 3, because of the 4 Hydrogen molecules.

2. tolerance

With this parameter you can specify the percentage tolerance of the density that you would like to be achieved for the bulk density of the solvent.

For the above example, passing a density of 0.6 and setting tolerance=1 will achieve a bulk solvent density of:

• (0.6 amu / Ang ^ 3) +/- 1 %

• equivalent to (0.6 +/- 0.006) amu / Ang ^ 3

Note: the tolerance has a default value of 1 %. Setting the tolerance to anything lower than increases the risk of solvate not converging to within the specified tolerance.

3. solvent

Using a solvent with pre-defined coordinates.

MDMC has a few in-built solvents that can be used to solvate your Universe. These have pre-defined atomic coordinates and interactions.

The in-built solvents are listed in the help method for solvate:

help(universe.solvate)

Help on method solvate in module MDMC.MD.simulation:

solvate(density: float, tolerance: float = 1.0, solvent: str = 'SPCE', max_iterations: int = 100, **settings: dict) -> None method of MDMC.MD.simulation.Universe instance
    Fills the ``Universe`` with solvent molecules according to pre-defined
    coordinates.

    Parameters
    ----------
    density : float
        The desired density of the ``Solvent`` that solvates the
        ``Universe``, in units of ``amu Ang ^ -3``
    tolerance : float, optional
        The +/- percentage tolerance of the density to be achieved.
        The default is 1 %. Tolerances of less than 1 % are at risk
        of not converging.
    solvent : str, optional
        A `str` specifying an inbuilt ``Solvent`` from the following:
        SPC, TIP3P, TIP3PFB, SPCE.
        The default is 'SPCE'.
    max_iterations: int, optional
        The maximum number of times to try to solvate the universe to
        within the required density before stopping. Defaults to 100.
    **settings
        ``constraint_algorithm`` (`ConstraintAlgorithm`)
            A ``ConstraintAlgorithm`` which is applied to the ``Universe``.
            If an inbuilt ``Solvent`` is selected (e.g. 'SPCE') and
            ``constraint_algorithm`` is not passed, the
            ``ConstraintAlgorithm`` will default to ``Shake(1e-4, 100)``.

    Raises
    ------
    ValueError
        If the ``Universe`` has already been solvated with a different
        density.

Example: solvating with SPCE water

universe.solvate(0.6, tolerance=1, solvent='SPCE')

4. **settings

constraint_algorithm

You can specify a ConstraintAlgorithm which is applied to the Universe. If specifying solvent as a string representing one of the in-built solvents (i.e. ‘SPCE’), then a Shake(1e-4, 100) constraint algorithm is applied by default.