# Analyzing MD trajectory: molecule orientation

I'm running a MD simulation (in Gromacs) of dye molecules incorporated into polymer to determine how that system responds to external electric field.

There are two order parameters that are used to describe such systems: $$<\cos\theta>$$ and $$<\cos^3\theta>$$, where $$\theta$$ being the angle between dipole moment of dye molecule and the direction of electric field, and averaging is over the ensemble of dye molecules.

I'm thinking of writing a script, that will find 3 specific atoms in each dye molecule, and with a help of linear algebra will determine orientation of that molecule. That way, assuming the dipole moment isn't changing during MD run, I'll be able to calculate that parameters.

The problem is I'm not very good in linear algebra and programming, so I'm wondering is there any software that can calculate such parameters?

• I have never used it, but gromacs has a rotmat module (gmx rotmat) which is supposed to determines the orientation of the frames from trajectory relative to a reference frame. If you use gmx trjconv combined with an index file, you could only keep the dye molecule in the trajectory and remove everything else. Then comparing it with the initial orientation of dye molecule should give you what you need. – Shoubhik R Maiti Mar 23 at 15:23
• @ShoubhikRMaiti gmx rotmat gives rotational matrix for only one molecule, and I need to average that angle over all molecules. Probably gmx gangle will work, but I don't quite understand how it works. – roma ichenko Mar 24 at 8:23
• If you calculate the rotational matrix R, that brings the two molecules in alignment, then you can apply the inverse R to your dipole moment you've pre-calculated. Then you can calculate the angle using arccos(E*dipole) – Cody Aldaz Mar 24 at 8:30
• I guess if gromacs does not have good inbuilt features for this, you will have to write your own script. There are many python modules like mdtraj (mdtraj.org/1.9.4/index.html) that help in analysis of trajectories, so you don't need to write code from scratch. As for the linear algebra, numpy and scipy will have almost everything you need. – Shoubhik R Maiti Mar 24 at 15:31
• @ShoubhikRMaiti thank you for the library, it would be a good occasion to revise my linear algebra knowledge. – roma ichenko Mar 25 at 3:59

## 1 Answer

So, with a lot of tries and errors I've found out how to perform such analysis. I'm using a library for Python called MDAnalysis. It can load Gromacs topologies and extract charges from them. Using them dipole moment for each molecule can be easily calculated and ordering parameters can be calculated from it.

Here is the code. It may not be perfect, but it gives $$<\cos \theta>$$ values consistent with the literature.

import MDAnalysis as md
import numpy as np

u=md.Universe('topol.top', '50ns/dr1.pdb', topology_format='ITP')

mean_cos_theta = []
mean_cos3_theta = []

z_unit = np.array([0,0,1])

for ts in u.trajectory:
mean_cos_theta_frame = []
mean_cos3_theta_frame = []
for res in u.residues:
dipole_moment = np.array([0,0,0])

for atom in res.atoms:
dipole_moment = dipole_moment + atom.charge*atom.position

cos_theta = np.dot(z_unit, dipole_moment)/np.linalg.norm(dipole_moment)
mean_cos_theta_frame.append(cos_theta)
mean_cos3_theta_frame.append(cos_theta**3)

mean_cos_theta.append(np.mean(mean_cos_theta_frame))
mean_cos3_theta.append(np.mean(mean_cos3_theta_frame))

• Nice solution! And you said you weren't very good with programming... – WaterMolecule Apr 1 at 16:40