This notebook outlines how to visualize the dipolar structure of 10% Sr-doped (Ba,Sr)TiO3 from bond valence molecular dynamics (BVMD) simulations. This recipe could be used to visualize the dynamics of any atom-centered vector property, e.g. magnetic moment, molecular orientation, etc. Please let me know if you have any questions, comments, or suggestions.

# Import modules I need
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
import glob
from scipy.interpolate import griddata
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

numpy offers efficient array operations, pandas provides efficient dataframe operations, matplotlib is used to make sophicated plots (e.g. animations), os gives command line control, glob is used to find files with particular names, and scipy has built-in functions for interpolation, which we’ll need to plot contours.

def count_lines(filename) :
    """function that counts the number of lines in a file named filename"""
    count = 0
    with open(filename, "r") as file :
        for line in file :
            count += 1
    return count

# Define user parameters

# Displacements
filename = "filesIwant/files.LocalP.dat100"
ntis = 1000
nstrs = int(count_lines(filename) / ntis)

# A-sites
datafilename = "data.BST"
data = np.array(os.popen("grep Atoms -A 1001 {} | tail -1000".format(datafilename)).read().split()).reshape((1000, 7)) # extracts Ti positions
subset = data[:, [2, -3, -2, -1]] # selects only the element id and cartesian coordinates (indices 2, -3, -2, -1 where negative indices count from the end of the list)
subset_df = pd.DataFrame(subset, columns = ["element", "x", "y", "z"]) # convert to a dataframe (which are easy to subset)
subset_df.element = subset_df.element.astype(float) # converts columns to floats
subset_df.x = subset_df.x.astype(float)
subset_df.y = subset_df.y.astype(float)
subset_df.z = subset_df.z.astype(float)
aplane = subset_df[subset_df.x < 2.5] # extracts a 10x10 slice of the structures for x < 2.5 A from the origin

filename is the name of a file containing the positions (x, y, z) and components (u, v, w) of the displacement (i.e. local dipole). You can download an example file here. The columns of this file correspond to x y z u v w. The rows correspond to a particular Ti’s displacement. There are nstrs structures and ntis Ti’s per structure; therefore, there are nstrs * ntis lines in the file.

datafilename is the name of a file containing the positions of the atoms. The structure is a 10x10x10 supercell of BST containing 10% Sr, and 5000 atoms total (hence ntis = 1000). You can download an example data file here.

# Read the Ti displacements into an array of dimension nstrs * ntis * 6, where 6 corresponds to x, y, z, u, v, w
# This may take a few seconds to complete
dip_traj = np.zeros((nstrs, ntis, 6))
with open(filename, "r") as file :
    for istr in range(nstrs) :
        for jti in range(ntis) :
            for line in file :
                dip_traj[istr, jti, :] = np.array(line.split()).astype(float)
                break

# Visualize dipolar structure dynamics

# Step 1: set up figure and axis
fig, ax = plt.subplots(figsize=(5, 5))
ax.set(xlim=(-5, 45), ylim=(-5, 45))

# Step 2: curate data
traj0 = dip_traj[0, :, :] # extract structure at time 0
df_traj0 = pd.DataFrame(traj0, columns = ["X", "Y", "Z", "U", "V", "W"]) # convert to a dataframe
plane0 = df_traj0[df_traj0.X < 2.5] # extract a 10x10 slice of the structures for x < 2.5 A from the origin
y = plane0.Y
z = plane0.Z
v = plane0.V
w = plane0.W
        
# Step 3: plot the Ti displacements and positions for time 0
qax = ax.quiver(y, z, v, w, scale = 4, zorder = 2, pivot = "middle")
ax.scatter(subset_df.y, subset_df.z, c = subset_df.element)

# Step 4: create a function that updates the displacements
def animate(i):
    """update the Ti displacements
    Note: the positions of the A-site atoms, Ba and Sr, are not updated (for simplicity)"""
    traji = dip_traj[i, :, :] # extract structure at time i
    df_traji = pd.DataFrame(traji, columns = ["X", "Y", "Z", "U", "V", "W"]) # convert to a dataframe
    planei = df_traji[df_traji.X < 2.5] # extract a 10x10 slice of the structures for x < 2.5 A from the origin
    v = planei.V
    w = planei.W
    qax.set_UVC(v, w) # update quiver plot
        
# Step 5: animate and show
anim = FuncAnimation(fig, animate, interval=100, frames=201)
#anim.save('dipolar_structure_dynamics.mp4')
HTML(anim.to_html5_video())

Purple and yellow dots correspond to the positions of Ba and Sr, respectively, at time 0. Black arrows corresponds to the dynamic Ti displacements. For 100 K, you can see (partially) that the polarization is pointing along [111], which is consistent with low-temperature, ferroelectric, rhombohedral structure of 10% Sr BST. For more information on quiver plots, please see matplotlib’s documentation.

References:

  1. Matplotlib animations the easy way by Ken Hughes
  2. Embedding Matplotlib Animations in Jupyter Notebooks by Louis Tiao