Brief look at MDAnalysis

Intro

I was going to do a brief tour of VLSI and some of the design steps, but the following blog by Anshul Thakur does a much better job than I can currently. So, instead, we are going to tear apart an example offered MDAnalysis. The goal is to explain like I'm five(ELI5) what they are doing. Here is the example:

import numpy as np
improt MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis


u = mda.Universe(PSF, DCD)


hbonds = HydrogenBondAnalysis(universe=u)
hbonds.run(verbose=True)


hbonds = HydrogenBondAnalysis(universe=u)
hbonds.hydrogens_sel = hbonds.guess_hydrogens("resname ARG HIS LYS")
hbonds.acceptors_sel = hbonds.guess_acceptors("resname ARG HIS LYS")
hbonds.run(verbose=True)
print(f"hydrogen_sel = {hbonds.hydrogens_sel}")
print(f"acceptors_sel = {hbonds.acceptors_sel}")


hbonds = HydrogenBondAnalysis(universe=u)

protein_hydrogens_sel = hbonds.guess_hydrogens("protein")
protein_acceptors_sel = hbonds.guess_acceptors("protein")

water_hydrogens_sel = "resname TIP3 and name H1 H2"
water_acceptors_sel = "resname TIP3 and name OH2"

hbonds.hydrogens_sel = f"({protein_hydrogens_sel}) or ({water_hydrogens_sel} and around 10 not resname TIP3)"
hbonds.acceptors_sel = f"({protein_acceptors_sel}) or ({water_acceptors_sel} and around 10 not resname TIP3)"


hbonds = HydrogenBondAnalysis(
    universe=u,
    between=[
        ["protein", "resname TIP3"], # for protein-water hbonds
        ["protein", "protein"]       # for protein-protein hbonds
    ]
)
hbonds.run(verbose=True)

Questions

What is mda.Universe?

It is a class that describes molecular dynamics systems by largely using atoms and trajectory attributes.^1

What is atoms

It is an AtomGroup, which in turn is a system of atoms. I think this also is what is called a "topology" file.^1

What is trajectory

Is a file that contains dynamic data like atom coordinates, box size, and velocities/forces.^2

What is HydrogenBondAnalysis(universe=u)?

It is a function that determines the hydrogen bond lifetimes. This is also done through the equation shown here. Briefly going through the equation, I believe that this equation is showing the time autocorrelation between time origin t0 and tau + t0. Autocorrelation is simply the coreelation between two observations in a time series.

What is hbonds.run do?

Simply computes the hydrogen bond lifetime using the above equation.^3

What does hbonds.guess_hydrogens do?

It guesses which hydrogen atoms to use.^3 Also, the first parameter is a selection string for the atom group. It seems that in order to do this the program simply looks at the mass and charge to see if it is within an acceptable range.

What does hbonds.guess_acceptors do?

Almost the exact same as above. This method only uses only charge, however.

What does the 'between' parameter do?

The between parameter at the very end is simply looking between two groups and computing the hydrogen bonds. in this case between protein-water hbonds and protein-protein hbonds.

Summary

In short, the only thing that feels like a magical black box is the HydrogenBondAnalysis. The only reason is because they reference a paper, which references a book, and I can't get ahold of the book so I'll have to take their word for it. I'm sure that this is a well established formula, and with a sufficient background this is very straightforward.