Controlling Tumbling Effects in a Molecular Dynamics Simulation.

Sulstice
4 min readDec 11, 2021

--

Tumbling is the rotational direction and speed of an object as it’s spinning.

Well, why is this useful? By controlling the spin you can control entropic states and perhaps reveal a stable surface area.

Usually, a big thing will spin more slowly and something that is more compact will spin faster. So in the context of proteins, a large protein will spin slower than a smaller protein. Sometimes it can make for an easier entry point for a drug if the system was more rigid. Another example is in polymers, if you would like the polymer to start aggregating with each other you would want the tumbling effects to be minimal so that they can create an accessible surface to connect to each other.

In Molecular Dynamic Simulations they often run their stuff in a solvent like water. Water movement can often control the tumbling effects of the system that lies within it. So if you can control the water you can control the tumbling effect.

Let’s initialize our system:

import MDAnalysis as mdatopology = os.path.join(current_directory, 'topology.psf')
trajectory = os.path.join(current_directory, 'trajectory.xtc')
universe = mda.Universe(topology, trajectory)

If you recall, here is one previous blog to get the end-to-end distance vector (Green), and another blog to get the dipole moment vector (Red). Using both these vectors we will use it to get the angle between them. The more the angle approaches 0 the less tumbling effect we will see.

By applying an electric field in the simulation the water molecules should align and orient the system in a fashion as well to the electric field vector (x, y, or z-axis). The only variable factor would be the strength of the field i.e 0.4 millivolts etc.

Set up standard variables:

starting_frame = 0
ending_frame = 50000
skip = 10
data = universe.trajectory[starting_frame:ending_frame:skip]
atom_group = universe.select_atoms('segid CHAIN_ID')

Get the dipole vectors of each timestep:

dipole_vectors = []for timestep in data: 

# Sum all the atoms individual dipole moments
total_atoms_dipole_vectors = []

for atom in atom_group:

atom_to_atom_dipole_vector = [
atom.position[0] * atom.charge, # X
atom.position[1] * atom.charge, # Y
atom.position[2] * atom.charge # Z
]
total_atoms_dipole_vectors.append(atom_to_atom_dipole_vector)

dipole_vectors.append(add_vectors(total_atoms_dipole_vectors)

Now that we have that, let’s get our end-to-end distance vectors, you’ll have to figure it out for your system:

Note: Select atoms first before looping in this one, otherwise it’s painfully slow.

domains = {}
domains['start_atom'] = universe.select_atoms('ATOM ID')
domains['end_atom'] = universe.select_atoms('ATOM ID')

Get the magnitude vectors:

magnitudes = []      for timestep in data:   
x0 = domain['start_atom'].positions[0][0]
x1 = domain['end_atom'].positions[0][0]
y0 = domain['start_atom'].positions[0][1]
y1 = domain['end_atom'].positions[0][1]
z0 = domain['start_atom'].positions[0][2]
z1 = domain['end_atom'].positions[0][2]

x = x1 - x0
y = y1 - y0
z = z1 - z0
magnitudes.append([x, y, z])

Alrighty, now let’s take the angle between the two.

from numpy import rad2degangles = []     for i in range(0, len(vectors)):         
angle = mda.lib.mdamath.angle(dipole_vectors, magnitude_vectors)
angles.append(rad2deg(angle))

Plot the distribution of the angles with different electric strengths:

# Distribution Plot of the Dipole Vector vs the End to End Distance Vectorimport plotly.figure_factory as ff colors = ["#A800FF", "#0079FF", "#00F11D", "#FF7F00", "#FF0900"] bin_size = .5
bin_width = ''group_labels = [
'High',
'Mid',
'Low',
'No Field',
]
# Create distplot with custom bin_size
fig = ff.create_distplot(
[
dipole_angles,
# dipole_electric_field_high_angles,
# dipole_electric_field_mid_angles,
# dipole_electric_field_low_angles,
],
group_labels,
bin_size=bin_size,
show_hist=False,
show_rug=False,
colors=colors
)
fig.update_layout(legend=dict(itemsizing='constant')) fig.update_layout(
legend=dict(
orientation="h",
yanchor="bottom",
y=1.02,
xanchor="right",
x=1.2,
font = dict(family = "Arial", size = 45),
bordercolor="LightSteelBlue",
borderwidth=2,
),
legend_title = dict(
font = dict(family = "Arial", size = 45))
)
fig.update_traces(line=dict(width=10)) fig.update_xaxes(
ticks="outside",
tickwidth=3,
tickcolor='black',
tickfont=dict(family='Arial', color='black', size=50),
title_font=dict(size=46, family='Arial'),
title_text='Angle (θ)',
ticklen=15,
range=[0, 95],
tickvals=[0,30, 60, 90, 120]
)
fig.update_yaxes(
ticks="outside",
tickwidth=3,
tickcolor='black',
title_text='Probability Density',
tickfont=dict(family='Arial', color='black', size=50),
title_font=dict(size=46, family='Arial'),
ticklen=15,
tickvals=[0.012, 0.024]
)
fig.update_layout(
title_text="Dipole Vector wrt E-Field Distribution",
title_font=dict(size=40, family='Arial'),
template='simple_white',
xaxis_tickformat = 'i',
bargap=0.2,
height=600,
width=1400
)

Here we go:

Interesting that we can start to align systems with the electric field and reduce their tumbling. If you watch your simulations closely, perhaps you will see less tumbling and more stability in your system.

No Field
High Field

--

--

Sulstice
Sulstice

No responses yet