What do random walk polymers look like?¶
In chemistry, and in particular in macromolecular chemistry, there are a bunch of increasingly complex models for how polymer chains in solution look like, based on the random walk model.
They start real easy, the freely jointed chain just assumes a fixed bond length, that's all. So two consecutive bonds may have arbitrary angles to each other. That's not very physical, since it would allow the chain to grow back on itself, but it's easy, and somewhat helpful. There are all kinds of fomulas explaining how such a chain would behave.
When chemists want to know it a bit better, they have to accept that the angle between two bonds will usually be nearly fixed, for example in the tetrahedral angle of $109.47°$, resulting in the freely rotating chain, meaning that the $360°$ rotation around the axis of the first bond is still free. This model is more realistic, but also more complex to describe.
The next step is to acknowledge that atoms/molecules attached to the chain will also interact, making some rotational angles more favorable than others. This results in the model of chains with restricted rotation, which often has zones in which there is a very regular arrangement of bonds, before the chain bends and twists on in a different direction. This is one of the most accurate models.
The exact details of these models can be found in section 2.2 of the SPFP_summary_clean.pdf document supplied with this notebook.
From a squiggly mess to a mathematically correct representation¶
Unfortunately, it is common practice to elaborate all these mathematical models and then to give a visual representation (a sketch) of them that looks like spaghetti at best. Turns out, that's pretty far from the truth! I was interested in what this actually looks like and loved the challenge to put it into code. This notebook is the result of it.
Requirements to run this notebook¶
To run this notebook you obviously need Jupyter Notebook and Python3 installed. You also require some common Python packages. More importantly, gnuplot, a plotting program, and ffmpeg a "command line video editor" need to be accessible from the terminal, in order to bring the 3D models to life in a short animation.
In [1]:
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import random
Let's start with a 2D model of a freely jointed chain¶
In [2]:
bond_length = 3 # Arbitrary, let's assume Angstrom
# Make 6 plots, since we're dealing with statistics, here
for j in range(6):
x = [0]
y = [0]
n = 50 # number of bonds
for i in range(n):
r = random.random()
x.append(x[i] + bond_length*np.sin(2*np.pi*r))
y.append(y[i] + bond_length*np.cos(2*np.pi*r))
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
# Do the same, with longer chains
for j in range(6):
x = [0]
y = [0]
n = 500 # number of bonds
for i in range(n):
r = random.random()
x.append(x[i] + bond_length*np.sin(2*np.pi*r))
y.append(y[i] + bond_length*np.cos(2*np.pi*r))
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
# And even longer chains!
for j in range(6):
x = [0]
y = [0]
n = 50000 # number of bonds
for i in range(n):
r = random.random()
x.append(x[i] + bond_length*np.sin(2*np.pi*r))
y.append(y[i] + bond_length*np.cos(2*np.pi*r))
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
Going to 3D¶
Getting a random point on the circumference of a circle was easy. To get a random point on a sphere is much more challenging, if these points are supposed to be randomly distributet! So, I got an algorithm from here, that works well.
The following cell tests if it reliably gives a distance of $1$, a first indicator, that it works.
In [3]:
for i in range(50):
rad = 2 * np.pi
r1 = random.random()
r2 = random.random()
x = np.sin(r1*rad)*np.sin(r2*rad)
y = np.cos(r1*rad)*np.sin(r2*rad)
z = np.cos(r2*rad)
print(x**2+y**2+z**2)
A better test would be to actually visualize the distribution, to validate that points don't agglomerate somewhere.
In [4]:
%matplotlib inline
x = []
y = []
z = []
# math taken from https://mathworld.wolfram.com/SpherePointPicking.html
for i in range(2000):
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(np.sqrt(1-r1**2)*np.cos(r2))
y.append(np.sqrt(1-r1**2)*np.sin(r2))
z.append(r1)
print('Sum of squared "distances" in the three spatial dimension. They should be near equal.')
print(np.sum(np.array(x)**2))
print(np.sum(np.array(y)**2))
print(np.sum(np.array(z)**2))
print('\n\nVisualization')
fig = plt.figure()
ax = fig.gca(projection='3d')
#ax.plot(np.arange(0, 10*math.pi, 0.1), arr_solution_x, arr_solution_v)
ax.plot(x, y, z, '.')
#plt.gca().set_aspect('equal')
plt.xlabel('x')
plt.ylabel('y')
ax.set_zlabel('z')
plt.show()
Since that looks good, we are ready for a
Freely jointed chain in 3D¶
I had some trouble making this look good using matplotlib, which is why I am resorting to gnuplot.
So first we generate the data and put it into external files
In [ ]:
# pythons 3d plots aren't good enough. Use gnuplot instead,
# but generate data here.
bond_length = 3 # A
for j in range(6):
x = [0]
y = [0]
z = [0]
n = 500 # number of bonds
# math taken from https://mathworld.wolfram.com/SpherePointPicking.html
for i in range(n):
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(x[i] + np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(y[i] + np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(z[i] + r1*bond_length)
f = open('freelyJointedChain_3D_500_{}.txt'.format(j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
for j in range(6):
x = [0]
y = [0]
z = [0]
n = 50000
# math taken from https://mathworld.wolfram.com/SpherePointPicking.html
for i in range(n):
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(x[i] + np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(y[i] + np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(z[i] + r1*bond_length)
f = open('freelyJointedChain_3D_50000_{}.txt'.format(j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
And then we make a freakin' movie out of this (careful, this will take a while...)
In [ ]:
import os
# gnuplot for the images
os.system('cmd /c gnuplot "freelyJointedChain_movie.gp"')
# Ffmpeg for the videos
# large file size for reasonable quality, but works with "Films & TV"
# mpeg4: -c:v libxvid -q:v 4
# default codec, but "pimped" here.
# H.264: -c:v libx264 -preset veryslow -tune animation -crf 32
## dumb players in H.264: -vf format=yuv420p # -> now works with "Films & TV" and probably QuickTime
for i in range(6):
os.system('cmd /c ffmpeg -y -i movie/freelyJointedChain_3D_500_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p -b:a 0k movie/afreelyJointedChain_3D_500_{}.mp4'.format(i, i))
for i in range(6):
os.system('cmd /c ffmpeg -y -i movie/freelyJointedChain_3D_50000_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p movie/afreelyJointedChain_3D_50000_{}.mp4'.format(i, i))
The movie¶
If you don't have the required tools installed, or don't want to spend the time, this what the animations look like:
In [5]:
%%HTML
<div style='display:flex; flex-direction:row; width: 100%; justify-content: space-between;'>
<div>
<h3>500 bonds</h3>
<video width="400px" controls loop>
<source src="afreelyJointedChain_3D_500_5.mp4" type="video/mp4">
</video>
</div>
<div>
<h3>50 000 bonds</h3>
<video width="400px" controls loop>
<source src="afreelyJointedChain_3D_50000_4.mp4" type="video/mp4">
</video>
</div>
</div>
Freely rotating chain in 2D¶
We repeat the procedure for the freely rotating chain. However, this model becomes really boring in 2D, since there cannot be much of a rotation in a flat molecule. This is the main reason I bothered with the 3D stuff.
In [6]:
%matplotlib inline
bond_length = 3 # A
bond_angle = 180 - 109.5 # °
bond_angle *= 2*np.pi/360 # rad
# Use a rotation matrix to generate to new bond out of the old one
rot1 = np.array([[np.cos(bond_angle), -np.sin(bond_angle)], [np.sin(bond_angle), np.cos(bond_angle)]])
rot2 = np.linalg.inv(rot1)
for j in range(6):
x = [0]
y = [0]
n = 50
# We must initialize the first bond in this model
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
# All consecutive bonds are build automatically
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r-.5 < 0:
flip = rot1
else:
flip = rot2
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_freelyRotating_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
for j in range(6):
x = [0]
y = [0]
n = 500
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r-.5 < 0:
flip = rot1
else:
flip = rot2
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_freelyRotating_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
for j in range(6):
x = [0]
y = [0]
n = 50000
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r-.5 < 0:
flip = rot1
else:
flip = rot2
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_freelyRotating_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
3D freely rotating chain¶
In 3D, this is way more interesting!
In [ ]:
## pythons 3d plots aren't good enough. Use gnuplot instead,
# but generate data here.
bond_length = 3 # A
bond_angle = 180 - 109.5 # °
bond_angle *= 2*np.pi/360 # rad
def rotation_matrix(n, alpha):
'''Creates a rotation matrix to rotate by a given angle around a given axis
Parameters:
n (array-like): A vector of length three that defines the rotation axis
alpha (float): The angle by which to rotate
Returns:
matrix (numpy.ndarray): The rotation matrix
'''
a = alpha
# renormalize n to avert numerical issues in long chains
n = np.array(n) / np.sqrt(n[0]**2+n[1]**2+n[2]**2)
# The rotational matrix becomes somewhat more intricate
matrix = np.array([[n[0]**2*(1-np.cos(a))+np.cos(a), n[0]*n[1]*(1-np.cos(a))-n[2]*np.sin(a), n[0]*n[2]*(1-np.cos(a))+n[1]*np.sin(a)],
[n[1]*n[0]*(1-np.cos(a))+n[2]*np.sin(a), n[1]**2*(1-np.cos(a))+np.cos(a), n[1]*n[2]*(1-np.cos(a))-n[0]*np.sin(a)],
[n[2]*n[0]*(1-np.cos(a))-n[1]*np.sin(a), n[2]*n[1]*(1-np.cos(a))+n[0]*np.sin(a), n[2]**2*(1-np.cos(a))+np.cos(a)]])
return matrix
def last_bond():
return np.array([x[-1]-x[-2], y[-1]-y[-2], z[-1]-z[-2]])
for j in range(6):
x = [0]
y = [0]
z = [0]
# Again, we must initialize the first bond
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(r1*bond_length)
n = 500
for i in range(n):
lb = last_bond()
lb /= np.linalg.norm(lb) / bond_length # renormalization to avert numerical issues
orth = np.array([1, 1, -(lb[0]+lb[1])/lb[2]]) # vector orthogonal to last bond
rot = rotation_matrix(orth, bond_angle)
# Tild next bond by defined bond angle
nb = rot @ lb
# The second rotation will be around the axis of the old bond
# by a random angle
rot2 = rotation_matrix(lb, random.random()*rad)
nb = rot2 @ nb
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
f = open('freelyRotating_3D_{}_{}.txt'.format(n, j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
for j in range(6):
x = [0]
y = [0]
z = [0]
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(r1*bond_length)
n = 50000
for i in range(n):
lb = last_bond()
lb /= np.linalg.norm(lb) / bond_length # buggy. Disable this and check
orth = np.array([1, 1, -(lb[0]+lb[1])/lb[2]])
rot = rotation_matrix(orth, bond_angle)
nb = rot @ lb
rot2 = rotation_matrix(lb, random.random()*rad)
nb = rot2 @ nb
#print('norm nb 2', np.linalg.norm(nb), '\n')
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
f = open('freelyRotating_3D_{}_{}.txt'.format(n, j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
Make a movie:
In [ ]:
import os
# gnuplot for the images
os.system('cmd /c gnuplot "freelyRotating_movie.gp"')
# Ffmpeg for the videos
# large file size for reasonable quality, but works with "Films & TV"
# mpeg4: -c:v libxvid -q:v 4
# default codec, but "pimped" here.
# H.264: -c:v libx264 -preset veryslow -tune animation -crf 32
for i in range(6):
pass
os.system('cmd /c ffmpeg -y -i movie/freelyRotating_3D_500_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p movie/afreelyRotating_3D_500_{}.mp4'.format(i, i))
for i in range(6):
os.system('cmd /c ffmpeg -y -i movie/freelyRotating_3D_50000_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p movie/afreelyRotating_3D_50000_{}.mp4'.format(i, i))
The movie¶
If you don't have the required tools installed, or don't want to spend the time, this what the animations look like:
In [7]:
%%HTML
<div style='display:flex; flex-direction:row; width: 100%; justify-content: space-between;'>
<div>
<h3>500 bonds</h3>
<video width="400px" controls loop>
<source src="afreelyRotating_3D_500_0.mp4" type="video/mp4">
</video>
</div>
<div>
<h3>50 000 bonds</h3>
<video width="400px" controls loop>
<source src="afreelyRotating_3D_50000_3.mp4" type="video/mp4">
</video>
</div>
</div>
Chain with restricted rotation in 2D¶
Also very boring
In [8]:
%matplotlib inline
bond_length = 3 # A
bond_angle = 180 - 109.5 # °
bond_angle *= 2*np.pi/360 # rad
# Probability of cis bond
P_trans = .8 # 0 <= P_trans <= 1
# Rotational matrices
rot1 = np.array([[np.cos(bond_angle), -np.sin(bond_angle)], [np.sin(bond_angle), np.cos(bond_angle)]])
rot2 = np.linalg.inv(rot1)
for j in range(6):
x = [0]
y = [0]
n = 50 # number of bonds
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
rot_old = 1 # "turned left"
# We are using monte carlo sampling for the statistics
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r <= P_trans:
if rot_old == 1:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
else:
if rot_old == 2:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_restrictedRotation_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
for j in range(6):
x = [0]
y = [0]
n = 500 # number of bonds
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
rot_old = 1 # "turned left"
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r <= P_trans:
if rot_old == 1:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
else:
if rot_old == 2:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_restrictedRotation_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
for j in range(6):
x = [0]
y = [0]
n = 50000 # number of bonds
r = random.random()
x.append(bond_length*np.sin(2*np.pi*r))
y.append(bond_length*np.cos(2*np.pi*r))
rot_old = 1 # "turned left"
for i in range(1, n):
last_bond = (x[i]-x[i-1], y[i]-y[i-1])
r = random.random()
if r <= P_trans:
if rot_old == 1:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
else:
if rot_old == 2:
flip = rot2
rot_old = 2
else:
flip = rot1
rot_old = 1
new_bond = flip @ last_bond
x.append(x[i] + new_bond[0])
y.append(y[i] + new_bond[1])
plt.plot(x, y)
plt.plot(x[0], y[0], 'o', markersize=12)
plt.plot(x[-1], y[-1], 'mo', markersize=12)
plt.xlabel('width $[\mathsf{\AA}]$')
plt.ylabel('height $[\mathsf{\AA}]$')
plt.gca().set_aspect('equal')
plt.savefig('RandomWalk_restrictedRotation_{}_{}.pdf'.format(n, j), bbox_inches='tight')
plt.show()
3D restricted rotation¶
This is where it gets really interesting! In 2D, we just had to define probabilities for flips in two directions. Discrete, easy. Here, the bond angle is continuous, yet there must be an underlying probability distribution. How to do that, as physically accurately as possible?
I've decided to spare no pains and firstly define a rotational potential. Energy as a function of angle of torsion. Then, the probability of finding a certain angle is given by the Boltzmann distribution. This has to be translated into a function which gives a torsional angle from a random variable, so that it reproduces the statistics correctly! Again, I resorted to monte carlo...
Boltzmann distribution¶
This is the formula that must be used: $$ p(\alpha)\propto \exp\left(-\frac{E(\alpha)}{RT}\right) $$
This is how it's done:
In [9]:
pi = np.pi
# Create the energy polynomial (fit)
# Interpolation points
Xs = [-pi, -2/3*pi, -1/3*pi, 0, 1/3*pi, 2/3*pi, pi]
Ys = np.array([0, 16, 3.8, 19, 3.8, 16, 0]) *1000
mesh = 1000 # number of mesh points
# build the matrix (system of linear equations)
# The interpolation will be defined by the positions
# of extrema (f(x)=y, f'(x)=0)
def row1(x): # f(x)=y
row = np.zeros(2*len(Xs))
for i in range(2*len(Xs)):
row[i] = x**i
return row
def row2(x): # f'(x)=0
row = np.zeros(2*len(Xs))
for i in range(1, 2*len(Xs)):
row[i] = i*x**(i-1)
return row
matrix = np.zeros((2*len(Xs), 2*len(Xs)))
for i in range(2*len(Xs)):
if i % 2 == 0:
matrix[i:] = row1(Xs[int(i/2)]) # f(x)=y
else:
matrix[i:] = row2(Xs[int((i-1)/2)]) # f'(x)=0
# fill up Ys to include extremum conditions
Y_full = np.zeros(2*len(Ys))
for i in range(len(Ys)):
Y_full[2*i] = Ys[i] # f(x)=y
Y_full[2*i+1] = 0 # f'(x)=0
# Solve the system of linear equations
coefs = np.linalg.solve(matrix, Y_full)
print('coefs', coefs)
# A list of possible alpha values, discretized
alphas = np.linspace(-pi, pi, mesh, endpoint=False)
# This defines the potential which we just created
def energy(coefs, alpha):
ret = 0
for i in range(len(coefs)):
ret += coefs[i] * alpha**i
return ret
plt.plot(alphas, [energy(coefs, x) for x in alphas])
plt.title('Our potential')
plt.xlabel('$\\alpha$ [rad]')
plt.ylabel('$E(\\alpha)$')
plt.show()
print('Not exactly beautiful, but going to work')
# Probability according to Boltzmann
T = 298 # temperature
R = 8.314 # universal gas constant
def Boltzmann(coefs, alpha, T):
return np.exp(-energy(coefs, alpha) / (R*T))
# The total probability must be 1...
norm = np.sum([Boltzmann(coefs, x, T) for x in alphas]) * 2*pi/mesh
plt.plot(alphas, [Boltzmann(coefs, x, T) for x in alphas]/norm)
plt.title('Probability distribution')
plt.xlabel('$\\alpha$ [rad]')
plt.ylabel('probability')
plt.show()
# This is our torsional angle as a function of random variable
p_alpha = np.cumsum([Boltzmann(coefs, x, T) for x in alphas]/norm)*2*pi/mesh
plt.plot(p_alpha, alphas)
plt.title('Torsional angle as function of random variable')
plt.xlabel('Random variable')
plt.ylabel('$\\alpha$ [rad]')
plt.show()
And efficient way to find the index for our discrete approximation from a continuous random variable¶
Well still use a float random variable, but the values we can work with are discrete. How do we find the appropriate index to access our discrete data? Divide and conquer!
In [ ]:
# make an efficient algorithm to find discrete x value to continuous x value
def search(x, Xs):
i = int(len(Xs)/2)
left = 0
right = len(Xs)-1
while True:
# border cases
if i == 0:
return 0
if i == len(Xs) - 1:
return i
# found it?!
if Xs[i] <= x and x < Xs[i+1]:
return i
# divide and conquer!
elif Xs[i] < x:
left = i
i += int((right-i)/2)
else:
right = i
i -= math.ceil((i-left)/2)
Creating the data¶
In [ ]:
## pythons 3d plots aren't good enough. Use gnuplot instead,
# but generate data here.
bond_length = 3 # A
bond_angle = 180 - 109.5 # °
bond_angle *= 2*np.pi/360 # rad
def rotation_matrix(n, alpha):
a = alpha
n = np.array(n) / np.sqrt(n[0]**2+n[1]**2+n[2]**2)
matrix = np.array([[n[0]**2*(1-np.cos(a))+np.cos(a), n[0]*n[1]*(1-np.cos(a))-n[2]*np.sin(a), n[0]*n[2]*(1-np.cos(a))+n[1]*np.sin(a)],
[n[1]*n[0]*(1-np.cos(a))+n[2]*np.sin(a), n[1]**2*(1-np.cos(a))+np.cos(a), n[1]*n[2]*(1-np.cos(a))-n[0]*np.sin(a)],
[n[2]*n[0]*(1-np.cos(a))-n[1]*np.sin(a), n[2]*n[1]*(1-np.cos(a))+n[0]*np.sin(a), n[2]**2*(1-np.cos(a))+np.cos(a)]])
return matrix
for j in range(0):
x = [0]
y = [0]
z = [0]
rad = 2 * np.pi
# In this model, we need to initialize the first two bonds
# and always keep track of the last two bonds
# bond 1
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(r1*bond_length)
# bond 2
last_bond = (x[-1], y[-1], z[-1])
lb = np.array(last_bond)
lb /= np.linalg.norm(lb) / bond_length
orth = np.array([1, 1, -(lb[0]+lb[1])/lb[2]])
rot = rotation_matrix(orth, bond_angle)
nb = rot @ lb
rot2 = rotation_matrix(lb, random.random()*rad)
nb = rot2 @ nb
last_last_bond = last_bond
last_bond = nb
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
n = 500
for i in range(n):
r = random.random()
nb = last_last_bond # new bond before rotation
index = search(r, p_alpha)
alpha = alphas[index] + pi
rot = rotation_matrix(last_bond, alpha)
nb = rot @ nb
last_last_bond = last_bond
last_bond = nb
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
# plotting
if True:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
ax.set_zlabel('$x_3$')
plt.show()
# to files
if True:
f = open('restrictedRotation_3D_{}_{}.txt'.format(n, j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
for j in range(6):
x = [0]
y = [0]
z = [0]
rad = 2 * np.pi
r1 = (random.random()-.5)*2 # [-1, 1]
r2 = random.random()*rad # [0, 2pi]
x.append(np.sqrt(1-r1**2)*np.cos(r2)*bond_length)
y.append(np.sqrt(1-r1**2)*np.sin(r2)*bond_length)
z.append(r1*bond_length)
last_bond = (x[-1], y[-1], z[-1])
lb = np.array(last_bond)
lb /= np.linalg.norm(lb) / bond_length # buggy. Disable this and check
orth = np.array([1, 1, -(lb[0]+lb[1])/lb[2]])
rot = rotation_matrix(orth, bond_angle)
nb = rot @ lb
rot2 = rotation_matrix(lb, random.random()*rad)
nb = rot2 @ nb
last_last_bond = last_bond
last_bond = nb
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
n = 50000
for i in range(n):
r = random.random()
nb = last_last_bond
index = search(r, p_alpha)
alpha = alphas[index] + pi
rot = rotation_matrix(last_bond, alpha)
nb = rot @ nb
last_last_bond = last_bond
last_bond = nb
x.append(x[-1] + nb[0])
y.append(y[-1] + nb[1])
z.append(z[-1] + nb[2])
# plotting
if True:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
ax.set_zlabel('$x_3$')
plt.show()
# to files
if True:
f = open('restrictedRotation_3D_{}_{}.txt'.format(n, j), 'w')
f.write('#x\ty\tz\n')
for i in range(len(x)):
f.write('{}\t{}\t{}\n'.format(x[i], y[i], z[i]))
f.close()
Make a movie:
In [ ]:
import os
# gnuplot for the images
os.system('cmd /c gnuplot "restrictedRotation_movie.gp"')
# Ffmpeg for the videos
# large file size for reasonable quality, but works with "Films & TV"
# mpeg4: -c:v libxvid -q:v 4
# default codec, but "pimped" here.
# H.264: -c:v libx264 -preset veryslow -tune animation -crf 32
for i in range(0):
os.system('cmd /c ffmpeg -y -i movie/restrictedRotation_3D_500_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p movie/arestrictedRotation_3D_500_{}.mp4'.format(i, i))
for i in range(6):
os.system('cmd /c ffmpeg -y -i movie/restrictedRotation_3D_50000_{}_%d.png -c:v libx264 -preset veryslow -tune animation -crf 32 -vf format=yuv420p movie/arestrictedRotation_3D_50000_{}.mp4'.format(i, i))
The movie¶
If you don't have the required tools installed, or don't want to spend the time, this what the animations look like:
In [10]:
%%HTML
<div style='display:flex; flex-direction:row; width: 100%; justify-content: space-between;'>
<div>
<h3>500 bonds</h3>
<video width="400px" controls loop>
<source src="arestrictedRotation_3D_500_1.mp4" type="video/mp4">
</video>
</div>
<div>
<h3>50 000 bonds</h3>
<video width="400px" controls loop>
<source src="arestrictedRotation_3D_50000_4.mp4" type="video/mp4">
</video>
</div>
</div>