Selecting which parameters to fit
For complex simulations, there may be a lot of parameters; a simple argon simulation has two parameters, and later on in this guide we will see a water simulation with a total of 8 - this only gets bigger as the molecules and forces become more complex.
Thus, a user may only want to refine a subset of all their parameters. In MDMC this can be done in two ways:
Fixing a Parameter
Parameter objects can be fixed. This can either be set when they are initialised (created) or changed for an existing Parameter. By default, an initialised Parameter is not fixed.
[1]:
from MDMC.MD.parameters import Parameter
charge = Parameter(value=0.5, name='charge', fixed=True, unit='e')
Attempting to change a fixed parameter produces a warning and does not change the parameter.
[2]:
charge.value = 12
print(f"The value of the charge parameter is: {charge.value}")
The value of the charge parameter is: 0.5 e
For parameters which already exist (for example, a parameter created by a potential), the parameter can be fixed by setting their fixed attribute to True.
[3]:
sigma = Parameter(1.0, name='sigma', unit='Ang')
print('Is sigma fixed: {}'.format(sigma.fixed))
sigma.fixed = True
print('Is sigma fixed: {}'.format(sigma.fixed))
Is sigma fixed: False
Is sigma fixed: True
Filtering parameters
Only parameters that are passed to Control (as fit_parameters) will be refined. While it is simplest to pass all parameters in a Universe to Control, it is also possible to filter out a subset. To demonstrate this, we use a universe filled with water under a SPCE force field, as used in running a simulation.
[4]:
from MDMC.MD import *
universe = Universe(dimensions=21.75, constraint_algorithm=Shake(1e-4, 100), electrostatic_solver=PPPM(accuracy=1e-5))
H1 = Atom('H')
H2 = Atom('H', position=(0., 1.63298, 0.))
O = Atom('O', position=(0., 0.81649, 0.57736))
H_coulombic = Coulombic(atoms=[H1, H2], cutoff=10.)
O_coulombic = Coulombic(atoms=O, cutoff=10.)
water_mol = Molecule(position=(0, 0, 0),
velocity=(0, 0, 0),
atoms=[H1, H2, O],
interactions=[Bond((H1, O), (H2, O), constrained=True),
BondAngle(H1, O, H2, constrained=True)],
name='water')
universe.fill(water_mol, num_density=0.03356718472021752)
O_dispersion = Dispersion(universe, [O.atom_type, O.atom_type], cutoff=10., vdw_tail_correction=True)
universe.add_force_field('SPCE')
Universe created with:
Dimensions [21.75 21.75 21.75]
There are 8 parameters in the universe:
[5]:
print(universe.parameters)
{'equilibrium_state (#23)': <Parameter
{ID: 23,
type: 'equilibrium_state',
value: 1.0 Ang,
unit: 'Ang',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>, 'potential_strength (#24)': <Parameter
{ID: 24,
type: 'potential_strength',
value: 4637.0 kJ / mol Ang ^ 2,
unit: 'kJ / mol Ang ^ 2',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>, 'equilibrium_state (#17)': <Parameter
{ID: 17,
type: 'equilibrium_state',
value: 109.47 deg,
unit: 'deg',
fixed: False,
constraints: None,
interactions_name: 'BondAngle',
functions_name: 'HarmonicPotential',
tied: False}>, 'potential_strength (#18)': <Parameter
{ID: 18,
type: 'potential_strength',
value: 383.0 kJ / mol rad ^ 2,
unit: 'kJ / mol rad ^ 2',
fixed: False,
constraints: None,
interactions_name: 'BondAngle',
functions_name: 'HarmonicPotential',
tied: False}>, 'epsilon (#37)': <Parameter
{ID: 37,
type: 'epsilon',
value: 0.6502 kJ / mol,
unit: 'kJ / mol',
fixed: False,
constraints: None,
interactions_name: 'Dispersion',
functions_name: 'LennardJones',
tied: False}>, 'sigma (#38)': <Parameter
{ID: 38,
type: 'sigma',
value: 3.166 Ang,
unit: 'Ang',
fixed: False,
constraints: None,
interactions_name: 'Dispersion',
functions_name: 'LennardJones',
tied: False}>, 'charge (#4)': <Parameter
{ID: 4,
type: 'charge',
value: 0.4238 e,
unit: 'e',
fixed: False,
constraints: None,
interactions_name: 'Coulombic',
functions_name: 'Coulomb',
tied: False}>, 'charge (#27)': <Parameter
{ID: 27,
type: 'charge',
value: -0.8476 e,
unit: 'e',
fixed: False,
constraints: None,
interactions_name: 'Coulombic',
functions_name: 'Coulomb',
tied: False}>}
So while all 8 parameters can be passed when initiliasing Control, they can also be filtered. Parameters objects have a number of convenience methods to assist with this:
[6]:
parameters = universe.parameters
help(parameters)
Help on Parameters in module MDMC.MD.parameters object:
class Parameters(builtins.dict)
| Parameters(init_parameters: "Optional[Union[Parameter, 'list[Parameter]']]" = None)
|
| A `dict-like` object where every element is a ``Parameter`` indexed by name,
| which contains a number of helper methods for filtering.
|
| Although ``Parameters`` is a `dict`, it should be treated like a `list` when writing to it;
| i.e. initialise it using a `list` and use `append` to add to it. These parameters can
| then be accessed by their name as a key.
|
| In short; Parameters writes like a list and reads like a dict.
|
| Parameters
| ----------
| init_parameters: ``Parameter`` or `list` of ``Parameter``s, optional, default None
| The initial ``Parameter`` objects that the ``Parameters`` object contains.
|
| Attributes
| ----------
| array: np.ndarray
| An alphabetically-sorted numpy array of the ``Parameter``s stored in this object.
|
| Method resolution order:
| Parameters
| builtins.dict
| builtins.object
|
| Methods defined here:
|
| __getitem__(self, key: 'str') -> "Union[Parameter, 'list[Parameter]']"
| Return self[key].
|
| __init__(self, init_parameters: "Optional[Union[Parameter, 'list[Parameter]']]" = None)
| Initialize self. See help(type(self)) for accurate signature.
|
| __setitem__(self, key: 'str', value: 'Parameter') -> 'NoReturn'
| Set self[key] to value.
|
| append(self, parameters: "Union['list[Parameter]', Parameter]") -> 'None'
| Appends a ``Parameter`` or list of ``Parameter``s to the dict,
| with the parameter name as its key.
|
| Parameters
| ----------
| parameters: ``Parameter`` or `list` of ``Parameter``s
| The parameter(s) to be added to the dict.
|
| filter(self, predicate: 'Callable[[Parameter], bool]') -> 'Parameters'
| Filters using a predicate
|
| Parameters
| ----------
| predicate : function
| A function that returns a `bool` which takes a ``Parameter`` as an
| argument.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which meet the condition of the predicate
|
| filter_atom_attribute(self, attribute: 'str', value: 'Union[str, float]') -> 'Parameters'
| Filters based on the attribute of ``Atom`` objects which have each
| ``Parameter`` applied to them
|
|
| Parameters
| ----------
| attribute : str
| An attribute of an ``Atom``. Attributes to match to must be either
| `float` or str.
| value : str, float
| The value of the ``Atom`` ``attribute``.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which are applied to an ``Atom`` object
| which has the specified ``value`` of the specified ``attribute``
|
| filter_function(self, function_name: 'str') -> 'Parameters'
| Filters based on the name of the ``InteractionFunction`` of each
| ``Parameter``
|
| Parameters
| ----------
| function_name : str
| The name of the ``InteractionFunction`` of ``Parameter`` objects to
| return, for example ``'LennardJones'`` or ``'HarmonicPotential'``.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which have a ``function`` with the
| specified ``function_name``
|
| filter_interaction(self, interaction_name: 'str') -> 'Parameters'
| Filters based on the name of the ``Interaction`` of each ``Parameter``
|
| Parameters
| ----------
| interaction_name : str
| The name of the ``Interaction`` of ``Parameter`` objects to return,
| for example ``'Bond'``.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which have an ``Interaction`` with the
| specified ``interaction_name``
|
| filter_name(self, name: 'str') -> 'Parameters'
| Filters by ``name``
|
| Parameters
| ----------
| name : str
| The ``name`` of the ``Parameter`` objects to return.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects with ``name``
|
| filter_structure(self, structure_name: 'str') -> 'Parameters'
| Filters based on the name of the ``Structure`` to which each
| ``Parameter`` applies
|
| Parameters
| ----------
| structure_name : str
| The name of a ``Structure``.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which are applied to a ``Structure``
| which has the specified ``structure_name``
|
| filter_value(self, comparison: 'str', value: 'float') -> 'Parameters'
| Filters by ``value``
|
| Parameters
| ----------
| comparison : str
| A `str` representing a comparison operator, ``'>'``, ``'<'``,
| ``'>='``, ``'<='``, ``'=='``, ``'!='``.
| value : float
| A `float` with which ``Parameter`` values are compared, using the
| ``comparison`` operator.
|
| Returns
| -------
| Parameters
| The ``Parameter`` objects which return a `True` when their values
| are compared with ``value`` using the ``comparison`` operator
|
| log_parameters(self) -> 'None'
| Logs all Parameters by ID
|
| ----------------------------------------------------------------------
| Readonly properties defined here:
|
| as_array
| The parameters in the object as a sorted numpy array.
|
| Returns
| -------
| np.ndarray
| An alphabetically-sorted array of parameter values in the object.
|
| ----------------------------------------------------------------------
| Data descriptors defined here:
|
| __dict__
| dictionary for instance variables
|
| __weakref__
| list of weak references to the object
|
| ----------------------------------------------------------------------
| Methods inherited from builtins.dict:
|
| __contains__(self, key, /)
| True if the dictionary has the specified key, else False.
|
| __delitem__(self, key, /)
| Delete self[key].
|
| __eq__(self, value, /)
| Return self==value.
|
| __ge__(self, value, /)
| Return self>=value.
|
| __getattribute__(self, name, /)
| Return getattr(self, name).
|
| __gt__(self, value, /)
| Return self>value.
|
| __ior__(self, value, /)
| Return self|=value.
|
| __iter__(self, /)
| Implement iter(self).
|
| __le__(self, value, /)
| Return self<=value.
|
| __len__(self, /)
| Return len(self).
|
| __lt__(self, value, /)
| Return self<value.
|
| __ne__(self, value, /)
| Return self!=value.
|
| __or__(self, value, /)
| Return self|value.
|
| __repr__(self, /)
| Return repr(self).
|
| __reversed__(self, /)
| Return a reverse iterator over the dict keys.
|
| __ror__(self, value, /)
| Return value|self.
|
| __sizeof__(...)
| D.__sizeof__() -> size of D in memory, in bytes
|
| clear(...)
| D.clear() -> None. Remove all items from D.
|
| copy(...)
| D.copy() -> a shallow copy of D
|
| get(self, key, default=None, /)
| Return the value for key if key is in the dictionary, else default.
|
| items(...)
| D.items() -> a set-like object providing a view on D's items
|
| keys(...)
| D.keys() -> a set-like object providing a view on D's keys
|
| pop(...)
| D.pop(k[,d]) -> v, remove specified key and return the corresponding value.
|
| If the key is not found, return the default if given; otherwise,
| raise a KeyError.
|
| popitem(self, /)
| Remove and return a (key, value) pair as a 2-tuple.
|
| Pairs are returned in LIFO (last-in, first-out) order.
| Raises KeyError if the dict is empty.
|
| setdefault(self, key, default=None, /)
| Insert key with a value of default if key is not in the dictionary.
|
| Return the value for key if key is in the dictionary, else default.
|
| update(...)
| D.update([E, ]**F) -> None. Update D from dict/iterable E and F.
| If E is present and has a .keys() method, then does: for k in E: D[k] = E[k]
| If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v
| In either case, this is followed by: for k in F: D[k] = F[k]
|
| values(...)
| D.values() -> an object providing a view on D's values
|
| ----------------------------------------------------------------------
| Class methods inherited from builtins.dict:
|
| __class_getitem__(...)
| See PEP 585
|
| fromkeys(iterable, value=None, /)
| Create a new dictionary with keys from iterable and values set to value.
|
| ----------------------------------------------------------------------
| Static methods inherited from builtins.dict:
|
| __new__(*args, **kwargs) class method of builtins.dict
| Create and return a new object. See help(type) for accurate signature.
|
| ----------------------------------------------------------------------
| Data and other attributes inherited from builtins.dict:
|
| __hash__ = None
For example, if only the charge parameters should be refined, parameters could be filtered by name:
[7]:
charges = parameters.filter_name('charge')
print(charges)
# charges is a Parameters object
print('\nThe class of charges is: {}'.format(type(charges)))
{'charge (#4)': <Parameter
{ID: 4,
type: 'charge',
value: 0.4238 e,
unit: 'e',
fixed: False,
constraints: None,
interactions_name: 'Coulombic',
functions_name: 'Coulomb',
tied: False}>, 'charge (#27)': <Parameter
{ID: 27,
type: 'charge',
value: -0.8476 e,
unit: 'e',
fixed: False,
constraints: None,
interactions_name: 'Coulombic',
functions_name: 'Coulomb',
tied: False}>}
The class of charges is: <class 'MDMC.MD.parameters.Parameters'>
As each filter returns a Parameters object, filters can be chained together. For example, to find the potential strengths of all bonds:
[8]:
bond_potential_strengths = parameters.filter_name('potential_strength').filter_interaction('Bond')
print(bond_potential_strengths)
# These operations are commutative
print('\nThe order these methods are'
' applied does not matter: {}'.format(bond_potential_strengths
== parameters.filter_interaction('Bond').filter_name('potential_strength')))
{'potential_strength (#24)': <Parameter
{ID: 24,
type: 'potential_strength',
value: 4637.0 kJ / mol Ang ^ 2,
unit: 'kJ / mol Ang ^ 2',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>}
The order these methods are applied does not matter: True
It is also possible to filter parameters based on the properties of the atoms to which they apply. For instance, we can filter the SPCE parameters so that only parameters of interactions on H atoms are shown:
[9]:
H_parameters = parameters.filter_atom_attribute('name', 'H')
print(H_parameters)
{'equilibrium_state (#23)': <Parameter
{ID: 23,
type: 'equilibrium_state',
value: 1.0 Ang,
unit: 'Ang',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>, 'potential_strength (#24)': <Parameter
{ID: 24,
type: 'potential_strength',
value: 4637.0 kJ / mol Ang ^ 2,
unit: 'kJ / mol Ang ^ 2',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>, 'equilibrium_state (#17)': <Parameter
{ID: 17,
type: 'equilibrium_state',
value: 109.47 deg,
unit: 'deg',
fixed: False,
constraints: None,
interactions_name: 'BondAngle',
functions_name: 'HarmonicPotential',
tied: False}>, 'potential_strength (#18)': <Parameter
{ID: 18,
type: 'potential_strength',
value: 383.0 kJ / mol rad ^ 2,
unit: 'kJ / mol rad ^ 2',
fixed: False,
constraints: None,
interactions_name: 'BondAngle',
functions_name: 'HarmonicPotential',
tied: False}>, 'charge (#4)': <Parameter
{ID: 4,
type: 'charge',
value: 0.4238 e,
unit: 'e',
fixed: False,
constraints: None,
interactions_name: 'Coulombic',
functions_name: 'Coulomb',
tied: False}>}
Finally there is also a more flexible method (Parameters.filter) which can be used in conjunction with any function to filter the parameters by some predicate function (a function which takes a parameter and returns True or False according to what is ‘filtered out’). This is most easily done with a lambda function but a regular function can also be used.
[10]:
def is_length(parameter):
"""A function which tells you whether the parameter measures length or not."""
return parameter.unit == 'Ang'
length_parameters = parameters.filter(is_length)
print(length_parameters)
# For those more familiar with Python, this can also be done using a lambda
lambda_length_parameters = parameters.filter(lambda p: p.unit == 'Ang')
print('\nThe same filter can be achieved using lambdas: {}'.format(lambda_length_parameters == length_parameters))
{'equilibrium_state (#23)': <Parameter
{ID: 23,
type: 'equilibrium_state',
value: 1.0 Ang,
unit: 'Ang',
fixed: False,
constraints: None,
interactions_name: 'Bond',
functions_name: 'HarmonicPotential',
tied: False}>, 'sigma (#38)': <Parameter
{ID: 38,
type: 'sigma',
value: 3.166 Ang,
unit: 'Ang',
fixed: False,
constraints: None,
interactions_name: 'Dispersion',
functions_name: 'LennardJones',
tied: False}>}
The same filter can be achieved using lambdas: True