Make animations of phonon displacement patterns

[mstekiel]

I'm working on the crystal system where I was interested in a particular phonon, and its displacement pattern. I adapted the functionality of euphonics to make gifs that show atomic vibrations of that phonon.

I introduced the calculate_displacement function that creates the displaced structure in a form of a Crystal object. It follows the definition of Martin Dove's book, that you cite for neutron structure factor. Then I also made a plot_crystal function that plots any Crystal on 3D axes.

From here, I was able to make such animation, with the code below.

I think such functionality might be interesting to you and other authors. I had a lot of fun making and investigating it :)

import numpy as np
import matplotlib.pyplot as plt
import sys

sys.path.append('C:/Users/Stekiel/Documents/GitHub/Euphonic')
from euphonic import ForceConstants
from euphonic.plot import plot_crystal_to_axis

# Input for euphonics
yaml_desciptor = 'phonopy_disp_handReduced.yaml'
fc_data = 'FORCE_CONSTANTS'


# Define the mode
Qpoints = np.array([[0.5, -0.5, 0.5]])
qpt_index = 0
mode = 2
ampl = 1

# More options for animation
frames = 9
SC = [3,3,3]

# Plotting styles
atom_styles = {
    'Ce':dict(label='Ce', s=250, color='black'),
    'Au':dict(label='Au', s=125, color='orange'),
    'Al':dict(label='Al', s=100, color='gray')
}

# Box and lattice
a, c = 4.3, 10.84
A = 0.3
bbox = dict(xmin=-A, xmax=1*a+A,
            ymin=-A, ymax=1*a+A,
            zmin=-A, zmax=2*c+A)


########################################
# Read force constants
fc = ForceConstants.from_phonopy(path='.',
                                summary_name=yaml_desciptor,
                                fc_name='FORCE_CONSTANTS')

# Calculate frequencies/eigenvectors to choose proper mode
phonons = fc.calculate_qpoint_phonon_modes(Qpoints, asr='reciprocal')

# Calculate displacements
crystal_snapshots = [phonons.calculate_displacement(qpt_index, mode, disp_amplitude=ampl*np.exp(2j*np.pi*n/frames), supercell=SC) for n in range(frames)]

# Generate images
for n,crystal_snapshot in enumerate(crystal_snapshots):
    print(f'Preparing frame {n}')

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(projection='3d')
    ax.grid(visible=None)
    ax.set_proj_type('ortho')

    ax.set_xlim(bbox['xmin']-A, bbox['xmax']+A)
    ax.set_ylim(bbox['ymin']-A, bbox['ymax']+A)
    ax.set_zlim(bbox['zmin']-A, bbox['zmax']+A)
    ax.set_box_aspect((4.3,4.3,10))

    # Plot displaced atoms
    plot_crystal_to_axis(crystal_snapshot, ax=ax, atom_styles=atom_styles, bbox=bbox)

    ax.legend()

    plt.savefig(f"animation/CeAuAl3-phonon-{n}.png", dpi=100)
    plt.close()

# Use pillow to save all frames as an animation in a gif file
from PIL import Image

images = [Image.open(f"animation/CeAuAl3-phonon-{n}.png") for n in range(frames)]

out_filename =  f'CeAuAl3-phonon.gif'
print(f'Preparing: {out_filename}')
images[0].save(out_filename, save_all=True, append_images=images, duration=60, loop=0)

[ajjackson]

This file (and the LICENSE) seem to have been mistakenly copied from the project directory into the Python module folder?

[mstekiel]

Yes

[ajjackson]

This takes me back to the now rather under-maintained ascii-phonons program; it even uses Pillow! 😄 I did find that in practice it can be challenging to choose a good viewing direction; specifying as a Miller index works quite well but often needs a few interative adjustments. https://ascii-phonons.readthedocs.io/en/latest/GUI.html

I have found Henrique Miranda's phonon website to be a superb tool for this stuff: https://henriquemiranda.github.io/phononwebsite/phonon.html Somewhere on the to-do list it seems a very good idea to provide a QpointPhononModes.to_web_json() method or similar, that produces the required data format for this website. euphonic-dispersion could then also provide an option to call this so we can easily get a Seekpath band structure, etc.

Another relevant tool for directly rendering images from crystal structure is ASE, which can create PNG, .gif or povray outputs as well as write to multi-frame file formats supported by viewers such as Jmol or Ovito. I'm not sure it's a good idea to duplicate that functionality here when it is quite simple to create an Atoms from a Crystal.

The calculate_displacement() method does seem very useful, though, and its applications go beyond visualisation. Perhaps we could keep that part, and include a tutorial that interfaces with ASE to write trajectory files and animated images?

[ajjackson]

Crystal has a reciprocal_cell() method for this purpose

[ajjackson]

This double-duty feels a bit awkward. Can you use complex numbers here as well? Is it useful to set disp_amplitude if also passing amplitudes here? Is it permitted to combine modes from multiple Q vectors?
It's helpful to minimise the number of paths through a function if possible; it makes the overall system easier to learn, document and maintain.

One approach might be to structure this as two public methods:
get_mode_displacement(qpt_index, mode, amplitude, supercell) and
get_combined_displacements(amplitude, supercell) where the latter amplitude array runs over all q-points and indices. These would each generate an appropriate (scaled) eigenvector then call a common private method _displaced_crystal(crystal, eigenvector, supercell).

Then again, the code to generate the eigenvectors is very simple. Perhaps a better API would be to add nothing at all to QpointPhononModes, and instead create a Crystal.apply_eigenvector(eigenvector, supercell) -> Crystal method on Crystal. That is more general and would support other (e.g. random) schemes of generating mass-scaled displacements. And to recreate these functions the user does

modes.crystal.apply_eigenvector(modes.eigenvectors[qpt_index, mode] * disp_amplitude, supercell)

or

weighted_evec = np.sum(modes.eigenvectors * weights, axis=(0, 1))
modes.crystal.apply_eigenvector(weighted_evec, supercell)

[mstekiel]

The backstory is that I was first implemented get_mode_displacement(qpt_index, mode, amplitude, supercell)-like function, then I wanted to look at a double degenerated phonon. The eigenvectors of the degenerated subspace were not to my liking, and I wanted to take their linear combination with get_combined_displacements(amplitude, supercell)-like function. In the end I managed to get what I want, but left the function in this awkward state.

With that experience, I feel like the apply_eigenvector function on Crystal is indeed a better choice. Feels more approachable and, with your snippets with the choice of the eigenvector, it is very clean.

[ajjackson]

@mstekiel any thoughts on this alternative approach?

[mstekiel]

@mstekiel any thoughts on this alternative approach?

I like it very much!

I think every person investigating lattice dynamics develops their own set of tools in the end :) Here, I wanted to be able to see the displacement to get an idea how a particular soft phonon looks like, within the packages provided by euphonics, Publication quality picture would require better tools, following those you listed in the first reply.

Another utility/Pandora box, would be a Crystal.to_cif() method, that allows to save the displaced structure in a standard cif format. My current approach is to export the supercell into a big cell with P1 space group and find the symmetry with FINDSYM of the ISODISTORT package.

[ajjackson]

Again, as a very active contributor to ASE my biased recommendation would be to toss the attributes of Crystal into ASE and write to one's favourite format from there. One would still need to use findsym or similar to get the symops back, it's true.

Even though it is only a few lines of code I suppose it would be helpful if that euphonic.Crystal <-> ase.Atoms interface existed somewhere. We don't really want to make ASE a dependency of euphonic. But I do intend that it is eventually possible to get ASE vibrations/phonons into Euphonic. Maybe a small ase_euphonic package should be responsible, or maybe it becomes another optional dependency of Euphonic and that stuff is tucked into its own module for safe-keeping 🤔

In any case it would be good for the Euphonic docs to have some nice tutorial examples that show people how to do this stuff, and can be used as the basis of user scripts.