|
| 1 | +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*- |
| 2 | +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 |
| 3 | +""" Module: topological_order_parameter |
| 4 | + =================================== |
| 5 | +
|
| 6 | + Implements the topological order parameter Psi as introduced in |
| 7 | + "Correlated Fluctuations of Structural Indicators Close to the |
| 8 | + Liquid−Liquid Transition in Supercooled Water", |
| 9 | + R. Foffi and F. Sciortino, J. Phys. Chem. B 2023, 127, 378−386 |
| 10 | + DOI: 10.1021/acs.jpcb.2c07169 |
| 11 | +
|
| 12 | +""" |
| 13 | +from __future__ import print_function |
| 14 | +from abc import ABCMeta, abstractmethod |
| 15 | +import numpy as np |
| 16 | +import warnings |
| 17 | +from scipy import stats, sparse |
| 18 | +from scipy.sparse.csgraph import dijkstra |
| 19 | +from scipy.spatial import KDTree |
| 20 | +from MDAnalysis.core.groups import Atom, AtomGroup, Residue, ResidueGroup |
| 21 | +from MDAnalysis.lib import distances |
| 22 | +from ..utilities_geometry import minimum_image |
| 23 | +from .observable import Observable |
| 24 | + |
| 25 | +class TopologicalOrderParameter(Observable): |
| 26 | + r"""Bond Topological order parameter. This observable provides |
| 27 | + also several helper functions and stores additional information that |
| 28 | + is useful for analyzing the topological properties of a set of molecules. |
| 29 | + """ |
| 30 | + |
| 31 | + def __init__(self, bonding_criterion, topological_distance=4): |
| 32 | + r""" |
| 33 | + :param BondingCriterion bonding_criterion: instance of a class to determine whether |
| 34 | + atoms in two groups are bonded or not. |
| 35 | + :param int topological_distance: which topological distance to consider when |
| 36 | + computing the order parameter |
| 37 | +
|
| 38 | + """ |
| 39 | + |
| 40 | + Observable.__init__(self, None) |
| 41 | + self.bonding_criterion = bonding_criterion |
| 42 | + self.topological_distance = topological_distance |
| 43 | + self.M = None |
| 44 | + |
| 45 | + def compute(self, inp ,kargs=None): |
| 46 | + r""" |
| 47 | + Compute the order parameters and Euclidean distances |
| 48 | + for atoms in group separated by the chosen topologigal |
| 49 | + distance |
| 50 | +
|
| 51 | + :param AtomGroup inp: an AtomGroup of including donors, |
| 52 | + acceptors and hydrogen |
| 53 | +
|
| 54 | + Returns: |
| 55 | +
|
| 56 | + psi : an ndarray containing the order parameters |
| 57 | + of each atom in group and an array with the indices |
| 58 | + of nearest neighbors at the chosen topological |
| 59 | + distance. |
| 60 | +
|
| 61 | + Note that the observable stores also other information |
| 62 | + in the class: |
| 63 | +
|
| 64 | + - neighbors: the index of the nearest neighbor at topological |
| 65 | + distance N. The distances to these neighbors are the values |
| 66 | + of psi. |
| 67 | +
|
| 68 | + - dists: the Euclidean distance to all neighbors at the |
| 69 | + chosen topological distance N |
| 70 | +
|
| 71 | + - predecessors: the matrix of predecessors. This can be |
| 72 | + used to obtain the sequence of molecules at topological |
| 73 | + distance 1,2, ..., N. |
| 74 | +
|
| 75 | + - pairs: the lists of pairs at topological distance N |
| 76 | +
|
| 77 | + - dist: the topological distances |
| 78 | +
|
| 79 | + Example: |
| 80 | +
|
| 81 | + >>> import numpy as np |
| 82 | + >>> import MDAnalysis as mda |
| 83 | + >>> import pytim |
| 84 | + >>> from pytim.observables import TopologicalOrderParameter,WaterHydrogenBondingAdjacencyMatrix |
| 85 | + >>> from pytim.datafiles import WATER_GRO |
| 86 | + >>> |
| 87 | + >>> u = mda.Universe(WATER_GRO) |
| 88 | + >>> pos = u.atoms.positions[:] |
| 89 | + >>> # First, we need to define the criterion |
| 90 | + >>> # used to define the bonding |
| 91 | + >>> HB = WaterHydrogenBondingAdjacencyMatrix() |
| 92 | + >>> # we create a molecule in the vapor phase, |
| 93 | + >>> # its value of the order parameter must be inf |
| 94 | + >>> pos[u.atoms.resindices==0]+=75 |
| 95 | + >>> u.atoms.positions=pos |
| 96 | + >>> g = u.select_atoms('name OW') |
| 97 | + >>> psi = TopologicalOrderParameter(HB,4) |
| 98 | + >>> with np.printoptions(precision=4): |
| 99 | + ... print(psi.compute(g)[:4]) |
| 100 | + [ inf 3.2736 3.9129 3.5797] |
| 101 | +
|
| 102 | + """ |
| 103 | + |
| 104 | + M = self.bonding_criterion.compute(inp) |
| 105 | + # we might use directed=True to speedup because we computed the matrix symmetrically |
| 106 | + res = dijkstra(M, |
| 107 | + unweighted=True, |
| 108 | + limit=self.topological_distance, |
| 109 | + directed=True, |
| 110 | + return_predecessors=True) |
| 111 | + self.dist,self.predecessors = res |
| 112 | + |
| 113 | + pairs = np.argwhere(self.dist==self.topological_distance) |
| 114 | + g = self.bonding_criterion.base_group |
| 115 | + dists = np.linalg.norm(minimum_image(g[pairs[:,0]].positions - g[pairs[:,1]].positions, |
| 116 | + self.bonding_criterion.group.universe.dimensions[:3]),axis=1) |
| 117 | + self.dists = dists |
| 118 | + # in some case (gas) a molecule might have no 4th neighbor, so we initialize |
| 119 | + # psi with inf. |
| 120 | + neighbors,psi=-np.ones(len(g),dtype=int),np.inf*np.ones(len(g),dtype=float) |
| 121 | + |
| 122 | + for i in np.arange(len(g)): |
| 123 | + # pairs has symmetric entries (eg if [1,2] is present, also [2,1] will) |
| 124 | + # because of how M has been built. These are the pairs that are connected |
| 125 | + # by a given topological_distance. Now we take the first column, of pairs |
| 126 | + # and we are sure, because of symmetry, that these are all the indices of |
| 127 | + # molecules that have a linked neighbour. We check for each of the donors |
| 128 | + # and acceptors if they are in this list. |
| 129 | + cond = np.where(pairs[:,0]==i)[0] |
| 130 | + if len(cond) > 0: # cond gives the list of donors or acceptors connected to |
| 131 | + # the i-th molecule of the donor+acceptor group. |
| 132 | + # select the specific pair with minimum distance |
| 133 | + amin = np.argmin(dists[cond]) |
| 134 | + # the neighbor is the second column in the pair where the distance is |
| 135 | + # minimum, and the associated value of psi is the minimum distance |
| 136 | + # itself. |
| 137 | + neighbors[i],psi[i] = pairs[cond,1][amin], np.min(dists[cond]) |
| 138 | + # otherwise, there is no 4th neighbor, we keep the default value (inf) |
| 139 | + self.psi = psi |
| 140 | + # so far neighbors and pairs are indices of the local array, now we store |
| 141 | + # them as MDAnalysis indices. To be able to search for paths, we store also |
| 142 | + # the indices of the local arrau |
| 143 | + self.neighbors = g.indices[neighbors] |
| 144 | + self._neighbors_localindex = neighbors |
| 145 | + self.pairs = self.bonding_criterion.base_group.indices[pairs] |
| 146 | + return self.psi |
| 147 | + |
| 148 | + def path(self, start, end): |
| 149 | + r""" |
| 150 | + Given the start and end of a path, returns the |
| 151 | + molecules within the connected path, if it exists. |
| 152 | + Accepts indices (zero-based) of the input group |
| 153 | + used to compute the observable, or two Atom objects, |
| 154 | + in which case it identifies the molecules to which they |
| 155 | + belong and returns an AtomGroup. |
| 156 | +
|
| 157 | + :param int|Atom : start, the first atom in the chain |
| 158 | + :param int|Atom : end, the last atom in the chain |
| 159 | + """ |
| 160 | + try: self.psi |
| 161 | + except (AttributeError): warnings.warn("The order parameter has not yet been computed",category=UserWarning) |
| 162 | + ret_type = list |
| 163 | + if isinstance(start,Atom) and isinstance(end,Atom): |
| 164 | + ret_type=AtomGroup |
| 165 | + if not start in self.group: raise ValueError("the start Atom does not belong to the input group") |
| 166 | + if not end in self.group: raise ValueError("the end Atom does not belong to input group") |
| 167 | + |
| 168 | + if ret_type == AtomGroup: |
| 169 | + start_id, end_id = start.resindex, end.resindex |
| 170 | + else: |
| 171 | + start_id, end_id = start, end |
| 172 | + p = [end_id] |
| 173 | + try: |
| 174 | + while p[-1] != start_id: |
| 175 | + p.append(self.predecessors[start_id,p[-1]] ) |
| 176 | + except: |
| 177 | + return [] |
| 178 | + p = np.array(p).tolist() |
| 179 | + if ret_type == AtomGroup: return self.bonding_criterion.base_group[p] |
| 180 | + else: return p |
| 181 | + |
| 182 | + def path_to_nearest_topological_neighbor(self,start): |
| 183 | + r""" |
| 184 | + Given a starting atom belonging to the input group |
| 185 | + of the observable, return the set of connected neighbors |
| 186 | + leading to the nearest neighbor at the chosen |
| 187 | + topological distance. |
| 188 | +
|
| 189 | + :param int|Atom : start, the first atom in the chain. |
| 190 | + If an Atom is passed, this must be part |
| 191 | + of either the acceptor or donor groups. |
| 192 | + If an index is passed, this is interpreted |
| 193 | + as the residue index of one of the acceptor |
| 194 | + or donor group members. |
| 195 | +
|
| 196 | + """ |
| 197 | + try: self.psi |
| 198 | + except (AttributeError): warnings.warn("The order parameter has not yet been computed",category=UserWarning) |
| 199 | + if isinstance(start, Atom): |
| 200 | + if not start in self.bonding_criterion.base_group: raise ValueError("the start Atom does not belong to either of the acceptor or donor groups") |
| 201 | + start_id = start.resindex |
| 202 | + else: |
| 203 | + start_id = start |
| 204 | + end_id = self._neighbors_localindex[start_id] |
| 205 | + path = self.path(start_id,end_id) |
| 206 | + if isinstance(start, Atom): return self.bonding_criterion.base_group[path] |
| 207 | + else: return path |
| 208 | + |
| 209 | + |
0 commit comments