Molecular orientations and rotations.¶
This notebook shows how to use GEMDAT to compute information related to molecular orientations and rotations.
First, it computes the bonds between central and satellite atoms, following the user's instructions. Then it computes their trajectory in cartesian, spherical and conventional form. If a symmetry group is defined, the statistics can be enhanced by exploiting the symmetry. Now the user can compute and plot the relevant information.
As input you will need:
- trajectory
- Orientational indication including:
- central atom species
- satellite atom species
- number of expected neighbors
- Optionally, a list of transformations that you want to apply to the unit vector trajectory. The user can combine the following operations:
normalize
: the trajectory of unit vectors is scaled to $[0,1]$transform
: convert the trajectory to conventional coordinate setting or otherwise. The conventional form may be simpler to visualize and compare.symmetrize
: applies symmetry operations specified via 4.
- Optionally, the point group, either as:
- Hermann–Mauguin notation
- Symmetry operations as 3x3 rotation matrices
The resulting orientation information can be visualized on a rectilinear plot, while the user can also access the statistics of the bond length. Finally, it is possible to plot the autocorrelation obtaining information about expected time for rotations.
Load the trajectory¶
We start by loading the trajectory from MD simulations.
from gemdat import Trajectory
# To use your own data:
# trajectory = Trajectory.from_vasprun('path/to/your/vasprun.xml')
from gemdat.utils import VASPCACHE_ORIENTATIONS as VASPCACHE
trajectory = Trajectory.from_cache(VASPCACHE)
trajectory
/home/stef/python/gemdat/.venv/lib64/python3.12/site-packages/MDAnalysis/topology/TPRParser.py:161: DeprecationWarning: 'xdrlib' is deprecated and slated for removal in Python 3.13 import xdrlib
Full Formula (Li16 S8 O32) Reduced Formula: Li2SO4 abc : 9.891300 9.891300 9.891300 angles: 60.000000 60.000000 60.000000 pbc : True True True Constant lattice (True) Sites (56) Time steps (6666)
Extract the trajectories of the unit vectors¶
The entry point to orientation finding is via an Orientations object. It computes the orientations and contains the orientation vectors.
from gemdat import Orientations
oris = Orientations(
trajectory,
center_type='S',
satellite_type='O',
)
oris
Orientations(trajectory=Full Formula (Li16 S8 O32) Reduced Formula: Li2SO4 abc : 9.891300 9.891300 9.891300 angles: 60.000000 60.000000 60.000000 pbc : True True True Constant lattice (True) Sites (56) Time steps (6666), center_type='S', satellite_type='O', vectors=array([[[ 0.47113781, 0.67949494, -1.33022098], [-1.58341442, 0.07208362, 0.07481367], [ 0.23692508, -1.47092223, 0.07815367], ..., [ 1.0762701 , 0.98581556, -0.20011628], [-0.57727755, -0.34050784, -1.40376095], [-1.11144327, 0.50560628, 0.86053202]], [[ 0.5133705 , 0.65825029, -1.3283309 ], [-1.584008 , 0.05856863, 0.05963306], [ 0.26752253, -1.46328622, 0.06763651], ..., [ 1.04147903, 1.03050253, -0.2146953 ], [-0.55834332, -0.30509991, -1.42183269], [-1.13202602, 0.51663661, 0.86796779]], [[ 0.55669039, 0.63464872, -1.31883796], [-1.57243621, 0.04706539, 0.05083524], [ 0.30466535, -1.45481033, 0.06230588], ..., [ 1.00815292, 1.07543505, -0.22493198], [-0.5393451 , -0.27013102, -1.42946754], [-1.1466608 , 0.52496499, 0.87855127]], ..., [[-0.69591073, -1.20342737, -0.75200444], [ 1.37375155, -0.04188055, -0.3618892 ], [-0.78866744, 1.24437287, -0.36940896], ..., [-0.40402687, -0.11054431, 1.40876756], [ 1.46986185, 0.44588816, 0.02082548], [-0.92918464, 1.05362794, -0.62993111]], [[-0.70978956, -1.18708479, -0.72122131], [ 1.36897242, -0.0544473 , -0.32346549], [-0.7943129 , 1.23325816, -0.36538725], ..., [-0.40091779, -0.10346154, 1.41105369], [ 1.4797438 , 0.45118816, 0.02199153], [-0.97405246, 1.05113743, -0.63068608]], [[-0.71984071, -1.16900686, -0.68794885], [ 1.3742471 , -0.06196398, -0.28769166], [-0.79134383, 1.22075506, -0.36058756], ..., [-0.39693313, -0.09806034, 1.41540952], [ 1.48652073, 0.45604486, 0.02042022], [-1.01219539, 1.04630722, -0.63173938]]]))
You can easily access the unit vectors in direct cartesian coordinates.
vectors = oris.vectors
vectors.shape
(6666, 32, 3)
transform, normalize, or symmetrize¶
For example, you can get the normalized conventional coordinates like this by chaining Orientations.normalize()
and Orientations.transform()
. The transformation matrix converts from primitive to conventional setting in this case.
import numpy as np
matrix = np.array(
[
[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
[1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5],
[0, 2 / 6**0.5, 1 / 3**0.5],
]
)
norm_conv = oris.normalize().transform(matrix)
It is also possible to set different transformation.
eye = np.eye(3)
# and then you can apply them
oris_conv = oris.transform(eye)
It is also possible to use symmetry transformations, which are very useful to improve the statistics of the trajectory. To do so, you first have to specify the point group of your molecule or cluster.
GEMDAT does this via pymatgen, following the Hemann-Mauguin notation. For example, here we select the Oh point group which is denoted as m-3m
oris_sym = oris.symmetrize(sym_group='m-3m')
You can also manually supply the rotation matrices.
sym_ops = np.asarray(
[
[[0, 1, 0], [1, 0, 0], [1, 0, 0]],
]
)
oris_manual_sym = oris.symmetrize(sym_ops=sym_ops)
It is also possible to convert any trajectory in spherical coordinates
oris_sym.vectors_spherical
Plotting¶
We compute the rectilinear plot which is a 2d map of the azimutal and elevation angles.
By default, it uses the transformed_trajectory
, so remember to specify which transformation you want to use.
Here we show the effect of different transformations their ordering.
for title, obj in (
('-', oris),
('norm', oris.normalize()),
('norm->symm', oris.normalize().symmetrize(sym_group='m-3m')),
('norm->conv', oris.normalize().transform(matrix)),
('norm->symm->conv', oris.normalize().symmetrize(sym_group='m-3m').transform(matrix)),
('norm->conv->symm', oris.normalize().transform(matrix).symmetrize(sym_group='m-3m')),
):
fig = obj.plot_rectilinear(normalize_histo=False)
fig.suptitle(title)
Plot the bond length distribution. A skewed gaussian is automatically fitted to the data in order to provide useful information.
fig1 = oris.plot_bond_length_distribution(bins=1000)
Now we compute and plot the autocorrelation of the unit vectors:
fig1 = oris.plot_autocorrelation()
Autocorrelation¶
Note that the user can calculate the autocorrelation function like this:
autocorrelation = oris.autocorrelation()
In this notebook, we want to compare the autocorrelation implemented in GEMDAT
using FFT, with a direct implementation that was developed for a previous matlab code. For more details look here.
For this comparison we define in the next block the matlab functions that we want to compare:
def simple_autocorr(direct_cart, Npt=100):
nts = len(direct_cart)
ep = nts - 1
dts = np.round(np.logspace(1, np.log10(ep), Npt)).astype(int)
mean_autocorr = np.zeros(len(dts))
for k, dt in enumerate(dts):
autocorr = np.sum(direct_cart[:-dt] * direct_cart[dt:], axis=-1)
mean_autocorr[k] = np.mean(autocorr)
return dts * 0.002, mean_autocorr # time in ps
from concurrent.futures import ThreadPoolExecutor
def autocorr_par(direct_cart, Npt=100):
def compute_autocorr(dt):
autocorr = np.sum(direct_cart[:-dt] * direct_cart[dt:], axis=-1)
mean_autocorr = np.mean(autocorr)
std_autocorr = np.std(autocorr)
return mean_autocorr, std_autocorr
nts = len(direct_cart)
ep = nts - 1
dts = np.round(np.logspace(1, np.log10(ep), Npt)).astype(int)
mean_autocorr = np.zeros(len(dts))
std_autocorr = np.zeros(len(dts))
with ThreadPoolExecutor() as executor:
results = list(executor.map(compute_autocorr, dts))
for i, (mean, std) in enumerate(results):
mean_autocorr[i] = mean
std_autocorr[i] = std
return dts * 0.002, mean_autocorr, std_autocorr
And we compare .autocorrelation()
with the matlab implementation.
import matplotlib.pyplot as plt
oris2 = oris.normalize().transform(matrix).symmetrize(sym_group='m-3m')
fig1 = oris2.plot_autocorrelation()
ax1 = plt.gca()
x, y = simple_autocorr(oris2.vectors)
ax1.scatter(x, y, marker='.', color='tab:red', label='Simple Autocorrelation')
x, y, y_std = autocorr_par(oris2.vectors)
ax1.errorbar(
x, y, zorder=-1, yerr=y_std, fmt='o', alpha=0.5, label='Mean Autocorrelation with Std Dev'
)
ax1.legend()
<matplotlib.legend.Legend at 0x7fe8f60b7320>