Merge pull request #344 from marrink-lab/bending

Implementation of MC bending sampling

polyply/src/build_file_parser.py

@@ -175,6 +175,15 @@ def _volume(self, line, lineno=0):
         resname, volume = line.split()
         self.topology.volumes[resname] = float(volume)
 
+    @SectionLineParser.section_parser('bending')
+    def _bending(self, line, lineno=0):
+        """
+        Parses the lines in the '[bending]'
+        directive and stores it.
+        """
+        resA, resB, resC, bending_const = line.split()
+        self.topology.bending[(resA, resB, resC)] = float(bending_const)
+
     def finalize_section(self, previous_section, ended_section):
         """
         Called once a section has finished. Here we perform all

polyply/src/nonbond_engine.py

@@ -11,12 +11,13 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+import random
 import itertools
 import numpy as np
 import scipy.spatial
 from polyply import jit
 from .topology import lorentz_berthelot_rule
-from .linalg_functions import not_exceeds_max_dimensions
+from .linalg_functions import angle, not_exceeds_max_dimensions
 
 def _lennard_jones_force(dist, point, ref, params):
     """
@@ -54,7 +55,6 @@ def _lennard_jones_force(dist, point, ref, params):
 
 POTENTIAL_FUNC = {"LJ": lennard_jones_force}
 
-
 def _n_particles(molecules):
     """
     Count the number of meta_molecule nodes
@@ -77,6 +77,8 @@ def __init__(self,
                  nodes_to_idx,
                  atom_types,
                  interaction_matrix,
+                 bending_matrix,
+                 torsion_matrix,
                  cut_off,
                  boxsize):
         """
@@ -93,6 +95,12 @@ def __init__(self,
              Dict mapping the atom_types to LJ type interaction parameters,
              that is sigma, epsilon or C6, C12 depending on the potential
              used. Currently only the sigma epsilon form is implemented.
+        bending_matrix: dict[frozenset(str, str), float]
+             Dict mapping the residue types to a bending constant for the
+             bending probability
+        torsion_matrix: dict[frozenset(str, str), float]
+             Dict mapping the residue types to a torsion constant for the
+             torsion probability
         cut_off: float
              cut-off for which to compute the interaction in nm
         boxsize: np.ndarray
@@ -103,6 +111,8 @@ def __init__(self,
         self.nodes_to_gndx = nodes_to_idx
         self.atypes = np.asarray(atom_types, dtype=str)
         self.interaction_matrix = interaction_matrix
+        self.bending_matrix = bending_matrix
+        self.torsion_matrix = torsion_matrix
         self.cut_off = cut_off
         self.boxsize = boxsize
 
@@ -292,6 +302,42 @@ def compute_force_point(self, point, mol_idx, node, exclude=[], potential="LJ"):
                     force += POTENTIAL_FUNC[potential](dist, point, self.positions[gndx_pair], params)
         return force
 
+    def compute_bending_probability(self, lp, point, mol_idx, node_b, node_c):
+        """
+        Compute probability of an angle between three points (`point`, `node_b`,
+        `node_c`) according to a simple decay function modulated by a bending
+        constant `lp`. If lp is large the distribution increases with the highest
+        probability being at 180 degrees (i.e. more straight angles), whilst if lp
+        is low the angles occur with almost equal probability.
+
+        Parameters
+        ----------
+        lp: float
+            bending constant
+        point: np.ndarray(1, 3)
+            coordinates of the new point
+        mol_idx: int
+            index of the molecule at hand
+        node_b: abc.hashable
+            first predecessor of the node corresponding to `point`
+        node_c: abc.hashable
+            second predecessor of the node corresponding to `point`
+
+        Returns
+        -------
+        float
+            the probability
+        """
+        # get the missing positions
+        B_pos = self.get_point(mol_idx, node_b)
+        C_pos = self.get_point(mol_idx, node_c)
+        # compute angle
+        ang_val = angle(point, B_pos, C_pos)
+        # compute probability
+        # denominator is normalization
+        prob = np.exp(lp*ang_val/180)/(180*(np.exp(lp)-1)/lp)
+        return prob
+
     @classmethod
     def from_topology(cls, molecules, topology, box):
         """
@@ -360,5 +406,7 @@ def from_topology(cls, molecules, topology, box):
         # dynamically set the cut-off as twice the largest vdw-radius
         cut_off = max(list(inter_matrix.values()))[0] * 2.
         nonbond_matrix = cls(positions, nodes_to_gndx,
-                             atom_types, inter_matrix, cut_off=cut_off, boxsize=box)
+                             atom_types, inter_matrix,
+                             cut_off=cut_off, boxsize=box,
+                             torsion_matrix=None, bending_matrix=topology.bending)
         return nonbond_matrix

polyply/src/random_walk.py

@@ -17,7 +17,7 @@
 from numpy.linalg import norm
 import networkx as nx
 from .processor import Processor
-from .linalg_functions import _vector_angle_degrees, not_exceeds_max_dimensions, norm_sphere, pbc_complete
+from .linalg_functions import angle, _vector_angle_degrees, not_exceeds_max_dimensions, norm_sphere, pbc_complete
 from .graph_utils import neighborhood
 from .meta_molecule import _find_starting_node
 """
@@ -246,6 +246,7 @@ def __init__(self,
         self.start_node = start_node
         self.nrewind = nrewind
         self.placed_nodes = []
+        self.prev_prob = 1
 
     def _rewind(self, current_step):
         nodes = [node for _, node in self.placed_nodes[-self.nrewind:-1]]
@@ -279,6 +280,53 @@ def checks_milestones(self, current_node, current_position, fudge=0.7):
 
         return True
 
+    def bendiness(self, point, node):
+        """
+        Perform Monte-Carlo like sampling of bending probability.
+
+        Parameters
+        ----------
+        point: np.ndarray
+            new coordinate
+        node: `abc.hashable`
+            the current node for which to get a new point
+        """
+        b_node = list(self.molecule.search_tree.predecessors(node))[0]
+        c_nodes = list(self.molecule.search_tree.predecessors(b_node))
+
+        if len(c_nodes) != 1:
+            return True
+
+        c_node = c_nodes[0]
+        # get the atom types aka resnames
+        typea = self.molecule.nodes[node]['resname']
+        typeb = self.molecule.nodes[b_node]['resname']
+        typec = self.molecule.nodes[c_node]['resname']
+        # get the bending constant
+        lp = self.nonbond_matrix.bending_matrix.get((typea, typeb, typec), None)
+        if not lp:
+            return True
+
+        prob = self.nonbond_matrix.compute_bending_probability(lp,
+                                                               point,
+                                                               self.mol_idx,
+                                                               b_node,
+                                                               c_node)
+        if self.prev_prob < prob:
+            self.prev_prob = prob
+            return True
+
+        # compute test probability from uniform sampling of prob function
+        norm = 180/lp * (np.exp(lp)-1)
+        test_prob = random.uniform(np.exp(lp*1/180)/norm,
+                                   np.exp(lp*179/180)/norm)
+
+        if test_prob < prob:
+            self.prev_prob = prob
+            return True
+
+        return False
+
     def update_positions(self, vector_bundle, current_node, prev_node):
         """
         Take an array of unit vectors `vector_bundle` and generate the coordinates
@@ -308,7 +356,8 @@ def update_positions(self, vector_bundle, current_node, prev_node):
             if fulfill_geometrical_constraints(new_point, self.molecule.nodes[current_node])\
                 and self.checks_milestones(current_node, new_point, step_length)\
                 and is_restricted(new_point, last_point, self.molecule.nodes[current_node])\
-                and not self._is_overlap(new_point, current_node):
+                and not self._is_overlap(new_point, current_node)\
+                and self.bendiness(new_point, current_node):
 
                     self.nonbond_matrix.add_positions(new_point,
                                                       self.mol_idx,
@@ -356,10 +405,8 @@ def _random_walk(self, meta_molecule):
         count = 0
         path = list(meta_molecule.search_tree.edges)
         step_count = 0
-
         while step_count < len(path):
             prev_node, current_node = path[step_count]
-
             if not meta_molecule.nodes[current_node]["build"]:
                 step_count += 1
                 continue

polyply/src/topology.py

@@ -215,6 +215,10 @@ class Topology(System):
         A dictionary of all typed parameter
     defines: list
         A list of everything that is defined
+    volumes: dict
+        Volume constants for meta-models
+    bending: dict
+        Sequence dependent bending constants for meta-model
     box: np.array(3,1)
         Box vectors as a, b, c in nanometers
     """
@@ -232,6 +236,7 @@ def __init__(self, force_field, name=None):
         self.persistences = []
         self.distance_restraints = defaultdict(dict)
         self.volumes = {}
+        self.bending = {}
         self.box = None
 
     def preprocess(self):

polyply/tests/test_nb_engine.py

@@ -172,6 +172,31 @@ def test_update_positions_in_molecules(topology):
          for node in mol.nodes:
             assert all(mol.nodes[node]["position"] == np.array([1., 1., 1.]))
 
+@pytest.mark.parametrize('lp, point, ref_prob',(
+                        # low lp, 180 degree,
+                        (0.001, np.array([0., 0., 2.]), 0.00555833379629605),
+                        # low lp, 90 degree,
+                        (0.001, np.array([0., 1., 1.]), 0.005555555324073842),
+                        # low lp, close to zero angle,
+                        (0.001, np.array([0., 0., 0.1]), 0.00555833379629605),
+                        # high lp, 180 degree,
+                        (10, np.array([0., 0., 2.]), 0.05555807788838943),
+                        # high lp, 90 degree,
+                        (10, np.array([0., 1., 1.]), 0.00037434738418303014),
+                        # high lp, close to 0 degree,
+                        (10, np.array([0., 0., 0.1]), 3.3299656173078084e-06),))
+def test_compute_bending_probability(topology, lp, point, ref_prob):
+    # we add two fixed positions to the first molecule
+    positions = [np.array([0., 0., 0,]), np.array([0., 0., 1.])]
+    for node, position in zip([0, 1], positions):
+        topology.molecules[0].nodes[node]["position"] = position
+
+    nb_engine =  NonBondEngine.from_topology(topology.molecules,
+                                             topology,
+                                             box=np.array([10., 10., 10.]))
+    prob = nb_engine.compute_bending_probability(lp, point, 0, 1, 0)
+    assert prob == pytest.approx(ref_prob, abs=10**-5)
+
 @pytest.mark.parametrize('mol_idx_a, mol_idx_b, node_a, node_b, expected',
                         ((0, 0, 0, 1, (0.53+0.67)/2.0),
                          (0, 2, 0, 0, (0.43+0.53)/2.0),

polyply/tests/test_random_walk.py

@@ -16,6 +16,7 @@
 """
 import math
 import pytest
+import random
 import numpy as np
 from numpy.linalg import norm
 import networkx as nx
@@ -188,6 +189,33 @@ def test_is_overlap(nonbond_matrix, molecule, new_point, result):
     # node 4 is already placed and hence is skipped over
     assert proccessor._is_overlap(new_point, 7, nrexcl=1) == result
 
+@pytest.mark.parametrize('n_coords, new_point, prev_prob, lp, result', (
+    # only 1 previous node -> always True
+    (1, None, 1, 10, True,),
+    # no lp is given -> always True
+    (1, None, 1, None, True,),
+    # 180 degrees should be fine even with high prev prob
+    (2, np.array([1., 1., 1.11]), 1, 10, True),
+    # small angle and prev prob high
+    (2, np.array([1., 1.47, 0.1]), 1., 10, False),
+    # small angle and prev prob low
+    (2, np.array([1., 1.47, 0.1]), 0., 10, True),
+    # medium angle, prev prob high, unifrom prob low
+    (2, np.array([1., 1.57, 1.74]), 1., 10, True),
+))
+def test_bendiness(nonbond_matrix, molecule, n_coords, new_point, prev_prob, lp, result):
+    # set random seed for reproducability
+    random.seed(1)
+    nb_matrix = add_positions(nonbond_matrix, n_coords)
+    # the bending constant for a series of PEO monomers
+    # the value is the same as in the high test case from
+    # test nb matrix
+    nb_matrix.bending_matrix = {("PEO", "PEO", "PEO"): lp}
+    processor = RandomWalk(mol_idx=0, nonbond_matrix=nb_matrix)
+    processor.molecule = molecule
+    processor.prev_prob = prev_prob
+    assert processor.bendiness(new_point, n_coords) == result
+
 @pytest.mark.parametrize('new_point, restraint, result', (
    # distance restraint true upper_bound
    # ref_node, upper_bound, lower_bound