Source code for dcma.biasing

"""
High-level interface to specify biasing potentials and forces that
were applied to the simulation.
"""

from __future__ import print_function

import numpy as np


class Bias(object):
    """
    A class to process command line options that handle biased simulations.
    """

    def __init__(self, pulling_force=None, bias_file=None, bias_function=None,
                 temperature=None):
        """
        Args:
            pulling_force: Pulling force in Piconewton.
            bias_file: A file containing a vector of biases (length: n_bins + 2).
            bias_function: A string containing a lambda function.
            temperature: in Kelvin, required to convert pulling_force into the
                right units

        """
        self._pulling_force = pulling_force
        self._bias_file = bias_file
        self._bias_function = bias_function
        self._temperature = temperature
        assert self._not_none_count() <= 1, \
            ("Error: the options --pulled, --bias_file, and " 
             "--bias_function are mutually exclusive")
        if pulling_force is not None:
            assert temperature is not None, \
                "Option --pulled requires temperature to be specified."
        if (temperature is not None) and (pulling_force is None):
            print("WARNING: --temperature argument is ignored (it is only "
                  "used for pulling forces).")

    def _not_none_count(self):
        return ((self._pulling_force is not None)
                + (self._bias_file is not None)
                + (self._bias_function is not None)
                )

    def get_biasing_type(self):
        if self._not_none_count() == 0:
            return "no"
        else:
            return "potential"

    def get_potential(self, unit_converter):
        """
        Args:
            unit_converter: An instance of :class:`~dcma.tools.UnitConverter`

        Returns:
            np.array: Biasing potential in units of kT.

        """
        if self._not_none_count() == 0:
            return None

        elif self._pulling_force is not None:
            assert self._temperature > 0
            generalized_centers = unit_converter.get_generalized_centers()
            # unit conversion pN -> kBT/Angstrom
            # U = - Fz  (in pN * A)
            # pn * A = 10^-12 J*A/m = 10^-22 J
            # 1 kB T = 1.38 *10^-23*T J / K => 1 J = 10^23 / 1.38 / (T/K) kB
            # pn * A = 10/1.38 kB K
            result = (- 10/1.3806485279 / self._temperature *
                      self._pulling_force * generalized_centers)
            return result

        elif self._bias_file is not None:
            result = np.loadtxt(self._bias_file)
            assert result.shape == (unit_converter.n_bins+2, )
            return result.transpose()

        elif self._bias_function is not None:
            # get centers, including "virtual bin centers"
            # beyond the left and the right edges
            generalized_centers = unit_converter.get_generalized_centers()

            # get objective function
            func = eval(self._bias_function)
            result = np.array([func(z) for z in generalized_centers])
            # no effective free energy difference over boundaries
            result[-1] = result[-2]
            result[0] = result[1]
            assert len(result) == unit_converter.n_bins + 2
            return result

        assert False, "Invalid position reached in biasing.py."