Freely rotating chain

Usage examples and full code reference for FreelyRotatingChain. For the underlying theory, see Freely rotating chain.

Quick start

from afrc.polymer_models.frc import FreelyRotatingChain

# defaults: bond length b = 3.8 A, characteristic ratio c_inf = 2.0
model = FreelyRotatingChain('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG')

model.get_mean_end_to_end_distance()
model.get_root_mean_squared_end_to_end_distance()
model.get_mean_radius_of_gyration()

r, p = model.get_end_to_end_distribution()
samples = model.sample_end_to_end_distribution(n=1000)

Tune the stiffness via the characteristic ratio (c_inf = 1 recovers the freely jointed chain; larger values give a stiffer, more extended ideal chain):

flexible = FreelyRotatingChain('MASNDYTQQATQSYG', c_inf=1.0)
stiff    = FreelyRotatingChain('MASNDYTQQATQSYG', c_inf=4.0)

flexible.get_root_mean_squared_end_to_end_distance()
stiff.get_root_mean_squared_end_to_end_distance()

See also the demo/demo_FreelyRotatingChain.ipynb notebook for a worked, plotted example.

Generating 3D ensembles

Generate explicit conformations (one bead per residue), write them to disk, and check how well they reproduce the model:

model = FreelyRotatingChain('MASNDYTQQATQSYGAYPTQPGQGYSQQSSQPYG', c_inf=2.0)

xyz = model.sample_conformations(n=5000, seed=1)                 # [n x N x 3] array, Angstroms
model.save_ensemble('ens', n=5000, seed=1)                       # ens.pdb + ens.xtc (needs mdtraj)
print(model.check_ensemble(xyz).format())                  # the model-check report

Or from the command line: afrc-ensemble -s <sequence> -m frc (see Command-line tools).

Code reference

class afrc.polymer_models.frc.FreelyRotatingChain(seq: str, p_of_r_resolution: float = 0.05, b: float = 3.8, c_inf: float = 2.0)

Freely rotating chain of N bonds of length b with a fixed bond angle and free torsions.

The freely rotating chain is an ideal chain: like the AFRC it has Gaussian end-to-end statistics with a scaling exponent of 0.5, but its size is set by a single stiffness parameter, the characteristic ratio \(C_\infty\). The mean-squared end-to-end distance is the exact finite-N result

\[\langle R^2 \rangle = C_\infty N b^2 - 2 b^2 \frac{\alpha (1 - \alpha^N)}{(1 - \alpha)^2}, \qquad \alpha = \frac{C_\infty - 1}{C_\infty + 1}\]

where \(\alpha\) is the cosine of the angle between successive bonds. c_inf = 1 recovers the freely jointed chain, and c_inf = 2 is a tetrahedral backbone. A freely rotating chain cannot reach the much larger characteristic ratio of a real polypeptide (\(C_\infty \approx 9\)), which comes from hindered rotation - use the AFRC for that.

This is a composition-independent reference model: the sequence is only used to set the number of bonds.

References

[1] Flory, P. J. (1969). Statistical Mechanics of Chain Molecules. Wiley-Interscience.

[2] Rubinstein, M., & Colby, R. H. (2003). Polymer Physics. Oxford University Press.

__init__(seq: str, p_of_r_resolution: float = 0.05, b: float = 3.8, c_inf: float = 2.0) → None

Create a FreelyRotatingChain object.

Parameters:

seq (str) – Amino acid sequence. Only its length (the number of bonds) is used, and it is not validated.
p_of_r_resolution (float) – Grid spacing (in Angstroms) used for the distribution. Default is 0.05 A.
b (float) – Bond length, in Angstroms. The default of 3.8 A is the Cα-Cα distance, i.e. one virtual bond per residue. Must be > 0.
c_inf (float) – Characteristic ratio \(C_\infty = (1 + \alpha)/(1 - \alpha)\), where \(\alpha\) is the cosine of the angle between successive bonds. c_inf = 1 recovers the freely jointed chain; c_inf = 2 (the default) is a tetrahedral bond angle. Must be > 0.

Raises:

FRCException – If b or c_inf is not positive.

check_ensemble(conformations: NDArray[float64]) → EnsembleReport

Check how well an ensemble reproduces the freely rotating chain.

The exact checks are the mean-squared distance of every residue pair (in bands of separation), the root-mean-square radius of gyration and first-to-last distance that follow from them, and the chain’s geometry: every bond exactly b long, every bond angle fixed by c_inf, and no pair further apart than its contour length. The model’s analytical (Gaussian) end-to-end distribution and its radius of gyration are reported as context.

Parameters:: conformations (np.ndarray) – Array of shape [n_conformations x N x 3] (in Angstroms).
Returns:: The checks; report.passed says whether the ensemble is model-like and report.format() gives a printable report.
Return type:: EnsembleReport
Raises:: AFRCException – If the conformations do not match this chain.

get_end_to_end_distribution()

Return the end-to-end distance distribution.

The distribution is computed on first use and then cached.

Returns:: (distances, probabilities), where distances are in Angstroms and the probabilities are a normalized probability mass function (they sum to 1).
Return type:: tuple of np.ndarray

get_mean_end_to_end_distance()

Return the mean end-to-end distance, \(\langle R_e \rangle\).

This is the expectation over the end-to-end distribution, \(\sum r P(r)\).

Returns:: The mean end-to-end distance (in Angstroms).
Return type:: float

get_mean_radius_of_gyration()

Return the root-mean-square radius of gyration, \(\sqrt{\langle R_g^2 \rangle}\).

This uses the ideal-chain relation \(\langle R_g^2 \rangle = \langle R_e^2 \rangle / 6\). Note that despite the method name this is the root-mean-square radius of gyration, not \(\langle R_g \rangle\). The name is kept for consistency with the other models.

Returns:: The root-mean-square radius of gyration (in Angstroms).
Return type:: float

get_mean_squared_distance_map() → NDArray[float64]

Return the exact mean-squared distance between every pair of residues.

For beads k bonds apart this is the exact finite-k freely rotating chain result (see the class docstring) - the same expression the model uses for the whole chain.

Returns:: Symmetric [N x N] matrix of mean-squared distances (in Angstroms squared), with zeros on the diagonal.
Return type:: np.ndarray

get_root_mean_squared_end_to_end_distance()

Return the root-mean-square end-to-end distance, \(\sqrt{\langle R_e^2 \rangle}\).

This is the square root of \(\sum r^2 P(r)\) over the end-to-end distribution. For long chains it approaches \(\sqrt{C_\infty} b \sqrt{N}\).

Returns:: The root-mean-square end-to-end distance (in Angstroms).
Return type:: float

sample_conformations(n: int = 1000, seed: int | Generator | None = None) → NDArray[float64]

Generate 3D conformations (one bead per residue) of the freely rotating chain.

Consecutive beads are joined by bonds of exactly b, every pair of consecutive bonds meets at the same bond angle (set by c_inf), and each torsion is uniformly random - the freely rotating chain itself. Its mean-squared distances are exact; its distance distributions are not Gaussian for short separations, unlike the model’s analytical P(r), which is a Gaussian with the right mean-squared size. Note that c_inf = 1 gives a fixed 90 degree bond angle, which has the same mean-squared size as the freely jointed chain but is a different chain.

Parameters:

n (int) – Number of conformations. Default is 1000.
seed (int, np.random.Generator or None) – Seed (or generator) for reproducible ensembles. If None (default) a fresh, unpredictable seed is used.

Returns:

Array of shape [n x N x 3] with the bead coordinates, in Angstroms, each conformation centred on the origin.

Return type:

np.ndarray

Raises:

AFRCException – If n is not a positive integer.

sample_end_to_end_distribution(n=1000)

Draw random end-to-end distances from the distribution.

Useful for building a size-matched sample to compare against simulation data.

Parameters:: n (int) – Number of values to draw. Default is 1000.
Returns:: n independent end-to-end distances (in Angstroms). For a zero-length chain every value is 0.
Return type:: np.ndarray

save_ensemble(filename: str, n: int = 1000, seed: int | Generator | None = None, pdb_only: bool = False) → NDArray[float64]

Generate a freely rotating chain ensemble and write it to disk.

Writes <filename>.pdb (topology and first conformation) and <filename>.xtc (every conformation; needs mdtraj), or with pdb_only=True a single multi-model <filename>.pdb. See sample_conformations() for how the conformations are generated.

Parameters:

filename (str) – Output path without an extension (a trailing .pdb or .xtc is dropped).
n (int) – Number of conformations. Default is 1000.
seed (int, np.random.Generator or None) – Seed (or generator) for reproducible ensembles.
pdb_only (bool) – Write a single multi-model PDB instead of a PDB/XTC pair. Default is False.

Returns:

The conformations that were written, shape [n x N x 3] (Angstroms).

Return type:

np.ndarray

Raises:

AFRCException – If n is invalid, the sequence has non-standard amino acids, mdtraj is missing (for a PDB/XTC pair), or the PDB format limits are exceeded.