Orientation as collective variables

Overview

pysages.colvars.orientation.RMSD(indices, ...)

Use a reference to calculate the RMSD of a set of atoms.

pysages.colvars.orientation.ERMSD(indices, ...)

Use a reference to calculate the eRMSD of a set of RNA structures.

pysages.colvars.orientation.ERMSDCG(indices, ...)

Use a reference to calculate the eRMSD of a set of Coarse-grained RNA structures.

Details

Collective variable for orientations describe the orientation measures of particles in the simulations with respect to a reference

RMSD describes distances between three-dimensional structures.

eRMSD describes distances between three-dimensional RNA structures. The standard RMSD is not accurate enough for distinguish RNA 3D structures. eRMSD provides a way to measuring the difference in the base interacting network. eRMSD also only requires the knowledge of C2, C4 and C6 for each base.

pysages.colvars.orientation.kabsch(P, Q)

Using the Kabsch algorithm with two sets of paired point P and Q, centered around the centroid. Each vector set is represented as an N x D matrix, where D is the the dimension of the space.

The algorithm works in three steps:

1. a centroid translation of P and Q (assumed done before this function call) 2. the computation of a covariance matrix C 3. computation of the optimal rotation matrix U

For more info see http://en.wikipedia.org/wiki/Kabsch_algorithm

Parameters:

  • P (array) – (N,D) matrix, where N is points and D is dimension.

  • Q (array) – (N,D) matrix, where N is points and D is dimension.

Returns:

U – Rotation matrix (D,D)

Return type:

matrix

pysages.colvars.orientation.rmsd(r1: Array, Q: Array, w_0: Array | None, optimal_rotation: Callable)

Calculate the RMSD respect to a reference using rotation matrix.

Parameters:

  • r1 (np.ndarray) – Atomic positions.

  • Q (np.ndarray) – Cartesian coordinates of the reference position of the atoms.

  • w_0 (np.ndarray) – Weights of the selected atom positions

  • optimal_rotation (Callable) – functions for the optimal rotation for aligning two positions

class pysages.colvars.orientation.RMSD(indices, references, weights=None, optimal_rotation: ~typing.Callable = <function kabsch>)

Use a reference to calculate the RMSD of a set of atoms. For more info see http://en.wikipedia.org/wiki/Kabsch_algorithm.

Parameters:

  • indices (list[int], list[tuple(int)]) – Select atom groups via indices.

  • references (list[tuple(float)]) – Cartesian coordinates of the reference position of the atoms in indices. The coordinates must match the ones of the atoms used in indices.

class pysages.colvars.orientation.ERMSD(indices, references, cutoff=2.4, a=0.5, b=0.3)

Use a reference to calculate the eRMSD of a set of RNA structures. Mathematical details can be found in [S. Bottaro, F. Di Palma, and G. Bussi, NAR, 2014](https://doi.org/10.1093/nar/gku972)

To use this method, the model has to have access to the coordinates of C2, C4, and C6 of each base. This is usually the case for all-atom simulations or coarse-grained models with SPQR-type of mappings.

Parameters:

  • indices (list[int].) – Must be a list of indices with dimension (n_nucleotides*3,). The order is tricky, for Purines (A, G) the order in the tuple should be C2, C6, C4; for Pyrimidines (U, C) the order in the tuple should be C2, C4, C6. This weird order is simply due to the numbering of the base. Please make sure the order is correct, otherwise the results will be wrong.

  • references (list[tuple(float)]) – Must be a list of tuple of floats, with dimension (n_nucleotides*3, 3) The last dimension is the x,y,z of the cartesian coordinates of the reference position of the C2/C4/C6 atoms, with the same order as the indices.

  • cutoff (float) – unitless cutoff. Default is 2.4

  • a (float) – re-scaling scale in x and y direction. Default is 0.5 (nm).

  • b (float) – re-scaling scale in z direction. Default is 0.3 (nm).

pysages.colvars.orientation.calc_local_reference_systems(rs)

calculate the local coordinate system for eRMSD calculation: origin and orthogonal x/y/z axis in the six membered rings. The x-axis is given by center of geometry (cog) -> C2, the z-axis is given by cross product between x-axis and cog -> second atom (C2 for U/C, C4 for A/G), and y-axis is given by cross product of z-axis and x-axis.

Parameters:

rs ((n_nucleotides, 3, 3) array) – positions of C2/4/6 atoms of every bases. The last dimension is x/y/z. The order is C2, C6, C4 for A/G and C2, C4, C6 for U/C.

Returns:

  • local_reference_systems ((n_nucleotides, 3, 3) array) – x, y, z unit vectors of each base

  • origins ((n_nucleotides, 3) array) – origins (center of geometry) of the six membered rings of each base

pysages.colvars.orientation.g_vector(local_reference_systems, origins, cutoff, a, b)

Calculate the smoothing function \(G(\tilde{r})\) based on local coordinates and origins for eRMSD calculations

Mathematical details can be found in [S. Bottaro, F. Di Palma, and G. Bussi, NAR, 2014](https://doi.org/10.1093/nar/gku972)

\(\tilde{\mathbf{r}}=(\frac{r_x}{a}, \frac{r_y}{a}, \frac{r_z}{a})\)

\(\gamma = \pi/\tilde{r}_\mathrm{cutoff}\)

\(\mathbf{G}(\tilde{\mathbf{r}})= \begin{pmatrix} \sin(\gamma |\tilde{r}|)\frac{\tilde{r}_x}{|\tilde{r}|}\\ \sin(\gamma |\tilde{r}|)\frac{\tilde{r}_y}{|\tilde{r}|}\\ \sin(\gamma |\tilde{r}|)\frac{\tilde{r}_z}{|\tilde{r}|}\\ 1+\cos(\gamma |\tilde{r}|) \\ \end{pmatrix} \frac{\Theta(\tilde{r}_\mathrm{cutoff}-|\tilde{r}|)}{\gamma}\)

Parameters:

  • local_reference_systems ((n_nucleotides, 3, 3) array) – arrays of coordinates of the xyz unit vectors

  • origins ((n_nucleotides, 3) array) – arrays of origins

  • cutoff (float) – unitless cutoff for eRMSD

  • a (float) – re-scaling factor in x, y direction, usually it’s 0.5 nm

  • b (float) – re-scaling factor in z direction, usually it’s 0.3 nm

Returns:

G

Return type:

(n_nucleotides, n_nucleotides, 4) array

pysages.colvars.orientation.reshape_coordinates(rs, reference)

reshape the coordinates so that they have easier dimensions to work with. In pysages, we have been concatenating the nucleotides with sites in each base. E.g [ C2, C4, C6, C2, C6, C4] for a 2nt RNA with sequence: UA. we want to reshape the coordinates so that it’s [[C2, C4, C6], [C2, C6, C4]]

Parameters:

  • rs ((n_nucleotides*3, 3) array) –

  • reference ((n_nucleotides*3, 3) array) –

Returns:

  • rs ((n_nucleotides, 3, 3) array)

  • reference ((n_nucleotides, 3, 3) array)

pysages.colvars.orientation.ermsd_core(rs, reference, cutoff, a, b)

compute the eRMSD given the current snapshots and the reference coordinates Mathematical details can be found in [S. Bottaro, F. Di Palma, and G. Bussi, NAR, 2014](https://doi.org/10.1093/nar/gku972)

First calculate the local coordinates based on C2/4/6 positions, and then calculate the G vectors. At last we sum up the square differences in G vectors.

\(\varepsilon RMSD = \sqrt{\frac{1}{N} \sum_{j, k} |\mathbf{G}(\tilde{\mathbf{r}}_{jk}^\alpha) - \mathbf{G}(\tilde{\mathbf{r}}_{jk}^\beta)|^2}\)

Parameters:

  • rs – (n_nucleotides, 3, 3) array, positions of the selected C2/4/6 atoms

  • reference – (n_nucleotides, 3, 3) array, reference positions

Returns:

eRMSD – value of the eRMSD

Return type:

float

pysages.colvars.orientation.ermsd(rs, reference, cutoff, a, b)

compute the eRMSD between the current snapshot and the reference

Parameters:

  • rs – (n_nucleotides*3, 3) array, positions of the selected C2/4/6 atoms

  • reference – (n_nucleotides*3, 3) array, reference positions

Returns:

eRMSD – value of the eRMSD

Return type:

float

class pysages.colvars.orientation.ERMSDCG(indices, reference, sequence, local_coordinates=[[[-0.1333021, -0.1762476, -1e-07], [-0.0857174, 0.0499809, -3.3e-06], [0.0795015, -0.1210861, -0.0001387]], [[-0.0253117, -0.122488, 0.0003341], [-0.0236469, 0.1239022, -0.0006548], [0.181128, 0.0, -0.0]], [[-0.116757, -0.1372227, -0.0002264], [-0.0409667, 0.0959927, -0.0012938], [0.1013132, -0.0983306, 0.0005419]], [[-0.0267224, -0.1188717, 0.0004939], [-0.025468, 0.1138277, 0.0008671], [0.1815304, -0.0, -0.0]]], cutoff=2.4, a=0.5, b=0.3)

Use a reference to calculate the eRMSD of a set of Coarse-grained RNA structures. Mathematical details can be found in [S. Bottaro, F. Di Palma, and G. Bussi, NAR, 2014](https://doi.org/10.1093/nar/gku972)

To use this method, we assume the model has access to 3 sites in the base. The method will use the ideal base geometry to reconstruct C2/C4/C6 that’s needed for eRMSD. The known base sites are assumed to be RACER-like mapping: A: C8, N6, C2 U: C6, O4, O2 G: C8, O6, N2 C: C6, N4, O2 Will refer to these mapped sites as B1/B2/B3

Parameters:

  • indices (list[int].) – Must be a list of indices with dimension (n_nucleotides*3,). The order is tricky. A: C8, N6, C2 U: C6, O4, O2 G: C8, O6, N2 C: C6, N4, O2 Please make sure the order is correct, otherwise the results will be wrong.

  • reference (list[tuple(float)]) – Must be a list of tuple of floats, with dimension (n_nucleotides*3, 3) The last dimension is the x,y,z of the cartesian coordinates of the reference position of the selected sites with the same order as the indices.

  • sequence (list[int]) – a list of int with length=n_nucleotides, each element is 0 (A) or 1 (U) or 2 (G) or 3 (C)

  • local_coordinates ((4,3,3) array) – Must be a list of tuple of floats, with dimension (4, 3, 3) Local coordinates of [C2, C6, C4] (A/G) or [C2, C4, C6] (U/C). The reference system is given by this: x-axis is unit vector, pointing from the center of B1/B2/B3 to B1 z-axis is cross-product between B1->B2, B1->B3, normalized y-axis is cross-product between z-axis and x-axis. The axis 0 of the local_coordinates is A/U/G/C. The axis 1 is [C2, C6, C4] (A/G) or [C2, C4, C6] (U/C). The axis 2 is coefficient for x-axis, y-axis and z-axis (Coefficient for z-axis is 0 because sites in the base share the same plane) The default value is local coordinates from RACER like mapping B1/B2/B3 to SPQR type of mapping. Notice the unit is in nm.

  • cutoff (float) – unitless cutoff. Default is 2.4

  • a (float) – re-scaling scale in x and y direction. Default is 0.5 (nm).

  • b (float) – re-scaling scale in z direction. Default is 0.3 (nm).

pysages.colvars.orientation.infer_base_coordinates(local_reference_systems, origins, sequence, local_coordinates)

Infer the coordinates of C2/C4/C6 based on the local reference system These are needed for eRMSD calculation.

Parameters:

  • rs ((n_nucleotides, 3, 3) array) – current coordinates of selected sites of the base

  • sequence ((n_nucleotides, ) array) – an array of int with length=n_nucleotides, each element is 0 (A) or 1 (U) or 2 (G) or 3 (C)

  • local_coordinates ((4,3,2) array) – Must be a list of tuple of floats, with dimension (4, 3, 2) Local coordinates of [C2, C6, C4] (A/G) or [C2, C4, C6] (U/C). The reference system is given by this: x-axis is unit vector, pointing from the center of B1/B2/B3 to B1 z-axis is cross-product between B1->B2, B1->B3, normalized y-axis is cross-product between z-axis and x-axis. The axis 0 of the local_coordinates is A/U/G/C The axis 1 is [C2, C6, C4] (A/G) or [C2, C4, C6] (U/C). The axis 2 is coefficient for x-axis, y-axis and z-axis (coefficient for z-axis is close to 0 because sites in the base are almost in the same plane) The default value is local coordinates from RACER like mapping B1/B2/B3 to SPQR type of mapping.

Returns:

C246_coords – coordinates of C2/4/6 in correct orders, ready for eRMSD calculation

Return type:

(n_nucleotides, 3, 3) array

pysages.colvars.orientation.ermsd_cg(rs, reference, sequence, local_coordinates, cutoff, a, b)

coarse-grained version of eRMSD which pass in different base sites other than the default C2/4/6 used by eRMSD. Thus we first need to infer the positions of C2/4/6 using rigid ideal base references and then do standard eRMSD calculations.

Parameters:

  • rs ((n_nucleotides*3, 3) array) – current snapshots of selected sites

  • reference ((n_nucleotides*3, 3) array) – reference of the selected sites The last dimension is the x,y,z of the cartesian coordinates of the reference position of the C2/C4/C6 atoms, with the same order as the indices.

  • sequence (list[int]) – a list of int with length=n_nucleotides, each element is 0 (A) or 1 (U) or 2 (G) or 3 (C)

  • local_coordinates ((4,3,2) array) – Must be a list of tuple of floats, with dimension (4, 3, 2) Local coordinates of [C2, C6, C4] (A/G) or [C2, C4, C6] (U/C).

  • cutoff (float) – unitless cutoff. Default is 2.4

  • a (float) – re-scaling scale in x and y direction. Default is 0.5 (nm).

  • b (float) – re-scaling scale in z direction. Default is 0.3 (nm).