Analyzing MD trajectory: molecule orientation

Question

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?

Answer 1

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))