Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Optimization and normal modes

In this section we will see how to visualize geometry optimization trajectories, as well as molecular vibration modes using py3Dmol.

Geometry optimization

import veloxchem as vlx
import py3Dmol as p3d
from matplotlib import pyplot as plt
import numpy as np
import ipywidgets
import h5py

Let’s assume we have run the geometry optimization and saved all the steps in xyz format in a file. To animate the optimization trajectory we first need a routine which reads the configurations from file.

def read_xyz_file(file_name):
    """Reads all the configurations from an xyz geometry optimization file."""
    xyz_file = open(file_name, "r")
    data = xyz_file.read()
    xyz_file.close()
    
    xyz_file = open(file_name, "r")
    i = 0
    steps = []
    energies = []
    for line in xyz_file:
        if "Energy" in line:
            steps.append(i)
            i += 1
            parts = line.split()
            energy = float(parts[-1])
            energies.append(energy)
    xyz_file.close()
    
    return data, steps, energies

Using this routine, we can plot the energies during optimization, as well as animate the trajectory.

xyz_file_name = '../../data/md/kahweol_optim.xyz'
xyz_data, steps, energies = read_xyz_file(xyz_file_name)

# Plot the energies
plt.figure(figsize=(7,4))
plt.plot(steps, energies,'o--')
plt.xlabel('Iteration')
plt.ylabel('Energy (H)')
plt.title("Kahweol XTB optimization")
plt.tight_layout(); plt.show()

# and animate the optimization
viewer = p3d.view(width=600, height=300)
viewer.addModelsAsFrames(xyz_data)
viewer.animate({"loop": "forward"})

viewer.setViewStyle({"style": "outline", "width": 0.05})
viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})

viewer.show()
<Figure size 700x400 with 1 Axes>
Loading...

If instead we would like to see each configuration individually, we can create an interactive widget with a slider. As in the previous section, we need a routine which reads an individual configuration at a time and returns the corresponding py3Dmol viewer.

def read_xyz_index(file_name, index=0):
    """Reads one configuration defined by its index from the xyz optimization file."""
    xyz_file = open(file_name, "r")
    current_index = 0
    data = ""

    # Read the number of atoms from the first line
    line = xyz_file.readline()
    natm = line.split()[0]
    data += line
    for line in xyz_file:
        if index == current_index:
            data += line
            if 'Energy' in line:
                parts = line.split()
                energy = float(parts[-1])
        if natm in line:
            parts = line.split()
            if len(parts) == 1:
                if index == current_index:
                    break
                current_index += 1
    xyz_file.close()
    return data, energy


def return_viewer(file_name, step=0):
    xyz_data_i, energy_i = read_xyz_index(file_name, step)
    
    # Uncomment if you would also like to see the energy plot
    # or comment out if you would like to skip this.
    plt.figure(figsize=(7,4))
    plt.plot(steps, energies, 'o--')
    plt.plot(step, energy_i, 'o', markersize=15)
    plt.xlabel('Iteration')
    plt.ylabel('Energy (H)')
    plt.title("Kahweol XTB optimization")
    plt.tight_layout(); plt.show()
    
    viewer = p3d.view(width=600, height=300)
    viewer.addModel(xyz_data_i)
    viewer.setViewStyle({"style": "outline", "width": 0.05})
    viewer.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
    viewer.show()
total_steps = steps[-1] # number of configurations in the trajectory file

# Note that the slider only works in a Jupyter notebook.
ipywidgets.interact(return_viewer, file_name=xyz_file_name,
                    step=ipywidgets.IntSlider(min=0, max=total_steps, step=1, value=3))

Vibrational normal modes

Now let’s see how to animate and visualize vibrational normal modes. For this, we will need to determine the molecular Hessian and perform a vibrational analysis.

xtb_drv = vlx.XtbDriver()
xtb_vibanalysis_drv = vlx.VibrationalAnalysis(xtb_drv)

In order to perform a vibrational analysis we need a Molecule object which stores the optimized geometry. To create the Molecule object we will read the last geometry from the xyz file we used in the previous subsection.

# read the last configuration from file
kahweol_xyz, opt_energy = read_xyz_index(file_name=xyz_file_name, index=total_steps)
# create the Mlecul object
xtb_opt_kahweol = vlx.Molecule.read_xyz_string(kahweol_xyz)
# perform vibrational analysis -- diagonalize Hessian, extract frequencies and normal modes
# and calculate IR intensities.
xtb_vibanalysis_drv.compute(xtb_opt_kahweol)

Now we have the frequencies, normal modes and IR intensities. In order to animate a specific normal mode, we need to write a routine which selects the mode of interest and returns it in the format required by py3Dmol. We would also like to plot the IR spectrum and, for this, we also need a routine which adds a Gaussian or Lorentzian broadening to the calculated IR spectrum.

# To animate the normal mode we will need both the geometry and the displacements 
def get_normal_mode(molecule, normal_mode):
    elements = molecule.get_labels()
    coords = molecule.get_coordinates_in_angstrom()
    natm = molecule.number_of_atoms()
    vib_xyz = "%d\n\n" % natm
    nm = normal_mode.reshape(natm, 3)
    for i in range(natm):
        # add coordinates:
        vib_xyz += elements[i] + " %15.7f %15.7f %15.7f " % (coords[i,0], coords[i,1], coords[i,2])
        # add displacements:
        vib_xyz += "%15.7f %15.7f %15.7f\n" % (nm[i,0], nm[i,1], nm[i,2])
    return vib_xyz

# Broadening function
def add_broadening(list_ex_energy, list_osci_strength, line_profile='Lorentzian', line_param=10, step=10):
        x_min = np.amin(list_ex_energy) - 50
        x_max = np.amax(list_ex_energy) + 50
        x = np.arange(x_min, x_max, step)
        y = np.zeros((len(x)))
        #print(x)
        #print(y)

        # go through the frames and calculate the spectrum for each frame
        for xp in range(len(x)):
            for e, f in zip(list_ex_energy, list_osci_strength):
                if line_profile == 'Gaussian':
                    y[xp] += f * np.exp(-(
                        (e - x[xp]) / line_param)**2)
                elif line_profile == 'Lorentzian':
                    y[xp] += 0.5 * line_param * f / (np.pi * (
                        (x[xp] - e)**2 + 0.25 * line_param**2))
        return x, y

After adding the broadening, we can plot the spectrum and then animate a specific normal mode, selected based on its index.

wvn, ir = xtb_vibanalysis_drv.vib_frequencies, xtb_vibanalysis_drv.ir_intensities
wvng, irg = add_broadening(wvn, ir, line_profile='Gaussian', line_param=10, step=2)
# Plot the IR spectra
plt.figure(figsize=(7,4))

# Plot the IR spectrum
plt.plot(wvng, irg, label='IR spectrum')
plt.xlabel('Wavenumber (cm**-1)')
#plt.axis(xmin=3200, xmax=3500)
#plt.axis(ymin=-0.2, ymax=1.5)
plt.ylabel('IR intensity (km/mol)')
plt.title("Calculated IR sepctrum, Kahweol")
plt.legend()
plt.tight_layout(); plt.show()

# Visualize the vibrational mode
index = 64
print("\n\n\n       Normal mode %d: %.2f cm-1, %.2f km/mol." % (index, 
                                                                xtb_vibanalysis_drv.vib_frequencies[index-1],
                                                                xtb_vibanalysis_drv.ir_intensities[index-1]))
normal_mode = get_normal_mode(xtb_opt_kahweol, xtb_vibanalysis_drv.normal_modes[index-1])
view = p3d.view(viewergrid=(1,1), width=600, height=300)
view.addModel(normal_mode, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})

view.setViewStyle({"style": "outline", "width": 0.05})
view.setStyle({"stick": {}, "sphere": {"scale": 0.25}})

view.animate({'loop': 'backAndForth'})
view.zoomTo()
view.show()
<Figure size 700x400 with 1 Axes>



       Normal mode 64: 1078.97 cm-1, 68.44 km/mol.
Loading...

To make a more interactive normal mode selection, we can create a slider using the ipywidgets Python module for Jupyter notebooks. To make use of the slider widget, we first need a routine which can plot the IR spectrum and animate a specific normal mode selected based on its index.

def vibration_viewer(index):
    freq = xtb_vibanalysis_drv.vib_frequencies[index-1]
    ir_intens = xtb_vibanalysis_drv.ir_intensities[index-1]
    # Plot the IR spectrum
    plt.figure(figsize=(7,4))
    plt.plot(wvng, irg)
    plt.plot(freq, ir_intens, 'o',  markersize=10)
    plt.xlabel('Wavenumber (cm**-1)')
    plt.ylabel('IR intensity (km/mol)')
    plt.title("Calculated IR sepctrum, Kahweol")
    plt.tight_layout(); plt.show()
    
    normal_mode = get_normal_mode(xtb_opt_kahweol, xtb_vibanalysis_drv.normal_modes[index-1])
    view = p3d.view(viewergrid=(1,1), width=600, height=300)
    view.addModel(normal_mode, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
    view.setViewStyle({"style": "outline", "width": 0.05})
    view.setStyle({"stick": {}, "sphere": {"scale": 0.25}})
    view.animate({'loop': 'backAndForth'})
    view.zoomTo()
    view.show()
no_norm_modes = xtb_vibanalysis_drv.vib_frequencies.shape[0] # number of vibrational modes

# Note that the slider only works in a Jupyter notebook.
ipywidgets.interact(vibration_viewer,
                    index=ipywidgets.IntSlider(min=1, max=no_norm_modes, step=1, value=92))
Visualization
MD trajectories
Visualization
Electronic structure