13

How can I visualize the normal mode vibrations of a molecule in Mathematica? Something like the following GIF:

vibrational modes of benzene (from here)

I have the XYZ file of the benzene molecule and the normal mode coordinates.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574

1 Answers1

29

Using a few undocumented functions that have been discussed here on the stack exchange before,

{atoms, coords} = Import[
    "https://www.dropbox.com/s/q3rwpedngv5wbqi/ben.xyz?dl=1",
    {"XYZ", {"VertexTypes", "VertexCoordinates"}}
];
mode = Most @ Import[
    "https://www.dropbox.com/s/3pjxams7ndr0l7h/19.mode?dl=1", "Table"
];
mode = 100 * mode;
bonds = Graphics`MoleculePlotDump`InferBonds[atoms, coords, 40, 25];
range = Range[-0.5, 0.5, 0.025];
range = Join[range, Reverse @ range];
plot[v_] := ImportExport`MoleculePlot3D[
    {
        "VertexCoordinates" -> (coords + v * mode),
        "VertexTypes" -> atoms, "EdgeRules" -> bonds
    }
];
ListAnimate @ Map[plot, range] 

enter image description here

When importing the mode file, it returns an empty list as the last element so I need to use Most to make it a matrix. Then I multiplied it by 100 to convert to picometers. Finally I had to precompute the bonds, because as the molecule moves along the normal coordinate, the atoms could get far enough apart that the bonds wouldn't be detected by MoleculePlot3D.

Sources:

  1. Connectivity in a molecule and permutations
  2. How can I access the internal function that plots a molecule from a formatted XYZ file?
  3. Visualizing .xyz file in Mathematica
Jason B.
  • 68,381
  • 3
  • 139
  • 286