A simple model of a polymer in solution treats its monomer units as totally uncorrelated in position (each monomer unit adopts a random orientation with respect to all the others): this is the *random flight model*.

In the following code, a random flight polymer with $N$ segments, each of length $a$, is constructed and displayed. The inertia matrix (actually a tensor) of this polymer can be calculated and diagonalized to find the principal axes. The moment of inertia can be represented as an inertia ellipsoid and this is plotted as a blue wireframe in the code below.

The inertia ellipsoid displayed in blue is that for the polymer at an instant in time. Considering all configurations that the polymer could adopt relative to its centre of mass, one can quantify its size as the *radius of gyration*, $R_g$: this is the root mean squared distance between all the segments and is also equal to the radius of a hollow sphere with the same moment of inertia as a theoretical random flight polymer. Since for a hollow sphere, $I = \frac{2}{3}mr^2$, the inertia ellipsoid corresponding to the radius of gyration (which is a sphere) can also be displayed (here, in red).

from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import numpy as np import numpy.linalg as linalg fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.set_aspect("equal") # Random flight segment length a = 1. # Number of random flight segments in the polymer N = 10000 # pos olds the locations of each segment (Cartesian coordinates) pos = np.empty((N,3)) pos[0] = (0,0,0) for i in range(1, N): # Pick a new random direction for the next segment uniformly on # the unit sphere theta = np.arccos(2 * np.random.random() - 1) phi = np.random.random() * 2. * np.pi # Add on the corresponding displacement vector for this segment. dpos = (a * np.sin(theta) * np.cos(phi), a * np.sin(theta) * np.sin(phi), a * np.cos(theta)) pos[i] = pos[i-1] + dpos # Shift the coordinate system to the centre of mass of the polymer cofm = np.mean(pos, axis=0) pos -= cofm # Calculate the moment of interia tensor mofi = np.empty((3,3)) x, y, z = pos[:,0], pos[:,1], pos[:,2] Ixx = np.sum(y**2 + z**2) Iyy = np.sum(x**2 + z**2) Izz = np.sum(x**2 + y**2) Ixy = -np.sum(x * y) Iyz = -np.sum(y * z) Ixz = -np.sum(x * z) mofi = np.array(((Ixx, Ixy, Ixz), (Ixy, Iyy, Iyz), (Ixz, Iyz, Izz))) # Find the rotation matrix and lengths of the semi-principal axes # by singular value decomposition U, s, rot = linalg.svd(mofi) abc = np.sqrt(s/N) # Rotate the polymer to align with its principal axes for i,(x,y,z) in enumerate(pos[:]): pos[i] = np.dot(np.array((x,y,z)), rot) # A mesh representing the unit sphere u = np.linspace(0.0, 2.0 * np.pi, 100) v = np.linspace(0.0, np.pi, 100) unit_sphere = np.array((np.outer(np.cos(u), np.sin(v)), np.outer(np.sin(u), np.sin(v)), np.outer(np.ones_like(u), np.cos(v)))).T # The ellipsoid of inertia ellipsoid = abc * unit_sphere ellipsoid = np.dot(ellipsoid, rot) # The radius of gyration of a theoretical random flight polymer Rg = np.sqrt(N * a**2 / 6) # Rg is the radius of a hollow sphere with the same moment of inertia as # the (theoretical) polymer; the ellipsoid of inertia of this theoretical # polymer is therefore (2/3) times this radius sphere = 2 / 3 * Rg * unit_sphere ax.plot(*pos.T, color="k") ax.plot_wireframe(*ellipsoid.T, rstride=4, cstride=4, color='b', alpha=0.3) ax.plot_wireframe(*sphere.T, rstride=4, cstride=4, color='r', alpha=0.3) ax.set_axis_off() plt.show()

Two different outputs are shown below for two different polymer configurations.

Share on Twitter Share on Facebook

## Comments

There are currently no comments

## New Comment