PyGRPY Documentation

PyGRPY is a python rewrite of Generalized Rotne Prager Yakamava hydrodynamic tensors package avaliable on github: GRPY. You can use it to obtain grand mobility matricies and to compute hydrodynamic sizes of particles in bead approximation.

New release also supports jax and jax.grad.

There is an online, interactive, variant (written in javascript) that you play with on my website ~rwaszkiewicz/rigidmolecule

How to install

python3 -m pip install pygrpy

and you’ll be good to go.

How to cite

Pychastic: Precise Brownian dynamics using Taylor-Itō integrators in Python
Radost Waszkiewicz, Maciej Bartczak, Kamil Kolasa, Maciej Lisicki
SciPost Phys. Codebases 11 (2023)
doi.org/10.21468/SciPostPhysCodeb.11.

Package contents

pygrpy.grpy.ensembleAveragedStokesRadius(ensemble_centres, sizes)

Returns stokes radius estimate given ensemble of configurations assuming no rigid bonds. See: https://doi.org/10.1017/jfm.2019.652 for details.

Parameters

  • ensemble_centres (np.array) – An M by N by 3 array describing M conformers each with N beads

  • sizes (np.array) – An array of length N describing radii of N beads.

Returns

Diffusive apparent hydrodynamic size of ensemble of conformers.

Return type

float

Example

>>> import numpy as np
>>> ensemble_centres = 0.1*np.arange(10*5*3).reshape(10,5,3)
>>> sizes = 0.5*np.arange(5) + 1
>>> pygrpy.grpy.ensembleAveragedStokesRadius(ensemble_centres,sizes)
3.0000008218278578
pygrpy.grpy.stokesRadius(centres, sizes)

Returns stokes radius of given conglomerate of beads treating it as a rigid body. See: https://doi.org/10.1016/j.bpj.2018.07.015 for details.

Parameters

  • centers (np.array) – An N by 3 array describing locations of centres of N beads.

  • sizes (np.array) – An array of length N describing radii of N beads.

Returns

Diffusive hydrodynamic size of the specified particle.

Return type

float

Example

>>> import numpy as np
>>> centres = np.array([[0,0,0],[0,0,1]])
>>> sizes = np.array([1,1])
>>> pygrpy.grpy.stokesRadius(centres,sizes)
1.1860816944024204
pygrpy.grpy_tensors.conglomerateMobilityMatrix(centres, radii)

Returns mobility matrix centred at [0,0,0] of bead conglomerate specified by centres and radii. Treats conglomerate as rigid body.

Parameters

  • centers (np.array) – An N by 3 array describing locations of centres of N beads.

  • radii (np.array) – An array of length N describing sizes of N beads.

Returns

A 6 by 6 array specifying mobility matrix centred at [0,0,0]. Indicies are ordered ux,uy,uz,wx,wy,wz.

Return type

np.array

pygrpy.grpy_tensors.mobilityCentre(mobility_matrix)

Returns mobility cetre realtive to zero centered mobility matrix mobility_matrix.

Parameters

mobility_matrix (np.array) – A 6 by 6 array specifying mobility matrix centred at [0,0,0]

Returns

A length 3 vector. Mobility centres location relative to [0,0,0].

Return type

np.array

pygrpy.grpy_tensors.mu(centres, radii, blockmatrix=False)

Returns grand mobility matrix in RPY approximation.

Parameters

  • centers (np.array) – An N by 3 array describing locations of centres of N beads.

  • radii (np.array) – An array of length N describing sizes of N beads.

  • blockmatrix ({True,False}) – Whether to retun rank 4 tensor instead of rank 2 tensor.

Returns

A 6N by 6N array containing translational mobility coefficients of the beads. Indicies are ordered: ux1,uy1,uz1,  ux2,uy2,uz2, ..., wx1,wy1,wz1, ...`. All translations before rotations, then by bead index, then by coordinate.

Unless blockmatrix is turned to True, then np.array with size=(2,2,N,N,3,3) is returned

Return type

np.array

pygrpy.grpy_tensors.muTT(centres, radii)

Returns grand mobility matrix in RPY approximation.

Parameters

  • centers (np.array) – An N by 3 array describing locations of centres of N beads.

  • radii (np.array) – An array of length N describing sizes of N beads.

Returns

A 3N by 3N array containing translational mobility coefficients of the beads. Indicies are ordered: ux1,uy1,uz1,  ux2,uy2,uz2, ...,`. by bead index, then by coordinate.

Return type

np.array

pygrpy.grpy_tensors.muTT_trace(centres, radii)

Returns beadwise trace of grand mobility matrix in RPY approximation.

Parameters

  • centers (np.array) – An N by 3 array describing locations of centres of N beads.

  • radii (np.array) – An array of length N describing sizes of N beads.

Returns

A N by N array containing traces of translational mobility coefficients of the beads.

Return type

np.array

pygrpy.grpy_tensors.muTT_trace_vectorised(list_of_centres, radii)

Returns beadwise trace of grand mobility matrix in RPY approximation. This version uses broadcasting for faster computation.

Parameters

  • list_of_centres (np.array) – A k by n by 3 array, where each n by 3 slice represents the locations of centres of n beads for a specific time step.

  • radii (np.array) – An array of length n describing sizes of n beads.

Returns

A k by n by n array containing traces of translational mobility coefficients of the beads for each time step in list_of_centres.

Return type

np.array

pygrpy.grpy_tensors.rigidProjectionMatrix(centres)

Returns rigid projection matrix. When this matrix is multiplied with friction matrix it gives friction matrix of rigid arrangement of beads.

Parameters

centers (np.array) – An N by 3 array describing locations of centres of N beads.

Returns

A 6 by 6N matrix containing projection coefficients.

Return type

np.array

pygrpy.pdb_loader.centres_and_radii(pdb_string)

Extracts the coordinates and radii of C-alpha atoms from a PDB file string.

Parameters

pdb_string (str) – The PDB file content as a string.

Returns

A tuple containing the C-alpha atom coordinates (np.array) and their corresponding radii (np.array).

Return type

tuple

Examples

>>> # Lysozyme C structure
>>> pdb_content = get_pdb_from_alphafold("P61626")
>>> coordinates, radii = centres_and_radii(pdb_content)
>>> print(stokesRadius(coordinates, radii))

>>> # Lysozyme C structure
>>> pdb_content = get_pdb_from_alphafold("P61626")
>>> coordinates, radii = centres_and_radii(pdb_content)
>>> print(stokesRadius(coordinates, radii))
pygrpy.pdb_loader.get_pdb_from_alphafold(uniprot, alphafold_download_server='https://alphafold.ebi.ac.uk/api/prediction/')

Fetches the PDB file from AlphaFold for a given UniProt ID.

Parameters

  • uniprot (str) – The UniProt ID of the protein.

  • alphafold_download_server (str [optional]) – The address of the alphafold server.

Returns

The PDB file content as a string if successful, None otherwise.

Return type

str or None

pygrpy.pdb_loader.get_pdb_from_pdb(pdb_id, pdb_download_server='https://files.rcsb.org/download/')

Fetches the PDB file of a given PDB ID.

Parameters

  • pdb_id (str) – The PDB ID of the protein.

  • pdb_download_server (str [optional]) – The address of the pdb server

Returns

The PDB file content as a string if successful, None otherwise.

Return type

str or None

Experimental features – jax support

pygrpy.jax_grpy_tensors.mu(centres, radii)

Returns grand mobility matrix in RPY approximation.

Parameters

  • centers (jnp.array) – An N by 3 array describing locations of centres of N beads.

  • radii (jnp.array) – An array of length N describing sizes of N beads.

Returns

A 6N by 6N array containing translational and rotational mobility coefficients of the beads. Indicies are ordered: ux1,uy1,uz1,  ux2,uy2,uz2, ..., wx1,wy1,wz1, ...`. All translations before rotations, then by bead index, then by coordinate.

Return type

jnp.array

Example

>>> import jax
>>> import jax.numpy as jnp
>>> import pygrpy.jax_grpy_tensors
>>> centres = jnp.array([[0,0,0],[0,0,1]])
>>> radii = jnp.array([1,1])
>>> pygrpy.jax_grpy_tensors.mu(centres,radii).shape()
(12,12)
>>> fast_mu = jax.jit(pygrpy.jax_grpy_tensors.mu)
>>> fast_mu(centres,radii).shape #compiled variant of mu
(12,12)
pygrpy.jax_grpy_tensors.muTT(centres, radii)

Returns grand mobility matrix in RPY approximation.

Parameters

  • centers (jnp.array) – An N by 3 array describing locations of centres of N beads.

  • radii (jnp.array) – An array of length N describing sizes of N beads.

Returns

A 3N by 3N array containing translational mobility coefficients of the beads. Indicies are ordered: ux1,uy1,uz1,  ux2,uy2,uz2, ...,`. by bead index, then by coordinate.

Return type

jnp.array

Example use

# Copyright (C) Radost Waszkiewicz 2022
# This software is distributed under MIT license
# Test if line of four identical beads has correct hydrodynamic size

import pygrpy
import numpy as np
import json

centres_four = np.array([[0,0,0],[0,0,1],[0,0,2],[0,0,3]])
sizes_four = np.array([1,1,1,1])

def test_hydrosize():
    testsize = pygrpy.grpy.stokesRadius(centres_four,sizes_four)
    assert np.allclose(testsize, 1.5409546371938094)

if __name__ == "__main__":
    test_hydrosize()

# Copyright (C) Radost Waszkiewicz 2024
# This software is distributed under MIT license
# Load shape of Lysozyme-C from different databases. Compare hydrodynamic size

import pygrpy.pdb_loader
import pygrpy.grpy

pdb_content = pygrpy.pdb_loader.get_pdb_from_alphafold("P61626")
coordinates, radii = pygrpy.pdb_loader.centres_and_radii(pdb_content)
alphafold_size = pygrpy.grpy.stokesRadius(coordinates, radii)

pdb_content = pygrpy.pdb_loader.get_pdb_from_pdb("253L")
coordinates, radii = pygrpy.pdb_loader.centres_and_radii(pdb_content)
pdb_size = pygrpy.grpy.stokesRadius(coordinates, radii)

print("Alphafold size [Ang]:")
print(alphafold_size)
print("Protein Data Bank size [Ang]:")
print(pdb_size)

# Copyright (C) Radost Waszkiewicz 2022
# This software is distributed under MIT license
# Load an ensemble from a .pdb file and compute R_h using locations of C_alpha atoms

import argparse
import numpy as np
import pygrpy
from tqdm import tqdm

# Console arguments.
parser = argparse.ArgumentParser()
parser.add_argument("-i", "--input", default="test.pdb", help="specify input file")
parser.add_argument(
    "-s",
    "--sigmas",
    help="compute standard deviations using bootstrap",
    action="store_true",
)

args = parser.parse_args()

with open(args.input, encoding="utf-8") as f:
    contents = f.read()

lines = contents.splitlines()

ensemble = list()
residues = list()
for line in lines:
    if "ATOM" in line:
        if "CA" in line:
            x = float(line[30:38])
            y = float(line[38:46])
            z = float(line[46:54])
            residues.append([x, y, z])
    elif "END" in line:
        ensemble.append(residues)
        residues = list()

ensemble = np.array(ensemble)
(ensemble_size, molecule_size, _) = ensemble.shape

hydrodynamic_size = pygrpy.grpy.ensembleAveragedStokesRadius(
    ensemble, 4.2 * np.ones(molecule_size)
)  # sizes in angstroms

centre_of_mass = np.mean(ensemble, axis=1)  # shape = (conformer,3)
gyration_radius = np.sqrt(3) * np.sqrt(np.mean((ensemble - centre_of_mass.reshape(-1, 1, 3)) ** 2))

bootstrap_rounds = 5
if args.sigmas:
    hydrodynamic_sizes_stats = np.zeros(bootstrap_rounds)
    for i in tqdm(range(bootstrap_rounds)):
        chosen = np.random.choice(np.arange(ensemble_size), ensemble_size)
        hydrodynamic_sizes_stats[i] = pygrpy.grpy.ensembleAveragedStokesRadius(
            ensemble[chosen], 4.2 * np.ones(molecule_size)
        )

    print(
        f"Hydrodynamic radius [Ang] = {hydrodynamic_size:.4f} +/- {np.std(hydrodynamic_sizes_stats):.4f}"
    )
    print(f"Gyration radius     [Ang] = {gyration_radius:.4f}")

else:
    print(f"Hydrodynamic radius [Ang] = {hydrodynamic_size:.2f}")
    print(f"Gyration radius     [Ang] = {gyration_radius:.2f}")