8 Auxiliary Tool User Guide

备注

Want to quickly analyze results, plot data, or perform common data processing tasks after completing a DFT calculation with DSPAW?

dspawpy (Python >= 3.9) is such a tool. It can be called programmatically (see example scripts below) and also provides a command-line interactive program.

After following the tutorial installation instructions, you can use the interactive program by typing dspawpy and pressing Enter in the command line:

 ... loading dspawpy cli ...

 This is the dspawpy command-line interactive tool.  Enjoy!
    ( )
   _| |  ___  _ _      _ _  _   _   _  _ _    _   _
 /'_` |/',__)( '_`\  /'_` )( ) ( ) ( )( '_`\ ( ) ( )
( (_| |\__, \| (_) )( (_| || \_/ \_/ || (_) )| (_) |
 \__,_)(____/| ,__/'`\__,_) \___x___/ | ,__/  \__, |
             | |                      | |    ( )_| |
             (_)                      (_)     \___/

 Version: Installation Path

 ======================================
 |  1: Update
 |  2: structure conversion
 |  3: Volumetric data processing
 |  4: Band structure calculation
 |  5: Density of States (DOS) data processing
 |  6: Joint display of band structure and density of states (DOS)
 |  7: Optical properties data processing
 |  8: NEB (Nudged Elastic Band) transition state calculation data processing
 |  9: phonon calculation data processing
 |  10: aimd ab initio molecular dynamics data processing
 |  11: Polarization data processing
 |  12: ZPE zero-point energy data processing
 |  13: Thermal correction energy of TS
 |
 |  q: Quit
 ======================================
 --> Enter a number and press Enter to select a function:

Highlights:

Autocompletion: Works by pressing the Tab key, helping to quickly and correctly enter the required program arguments.
Multithreaded lazy loading: Loads modules in the background while waiting for user input, significantly reducing waiting time; loads only necessary modules, minimizing memory usage.

Note:

When using on a remote server, the startup time may be longer due to poor disk I/O performance, potentially taking up to half a minute in extreme cases (directly related to the server's current disk I/O performance). If this is unacceptable, please install and use dspawpy on your own computer.
After typing dspawpy and pressing Enter, Python will first load built-in modules. Once this is complete, the prompt "... loading dspawpy cli ..." will appear, indicating the second stage (loading third-party dependencies) has begun.
After the second stage is completed, a welcome screen will be displayed, indicating that dspwapy has finished the initial loading and has entered the third stage. Subsequently, it will dynamically load the corresponding dependency libraries based on the selected functional modules, thereby minimizing waiting time.

8.1 Installation and Updates

On the HZW machine, dspawpy has been pre-installed. Activate the virtual environment using the following command to start using it:

source /data/hzwtech/profile/dspawpy.env

On other machines, please install dspawpy yourself (choose one of the following two methods):

Using mamba or conda, you can install the package from https://conda-forge.org/download/.

mamba install dspawpy -c conda-forge
#conda install dspawpy -c conda-forge

Or, use pip3 (some operating systems may not have the executable pip3, in which case try pip)

pip3 install dspawpy

For information on how to configure pip and conda mirror addresses to speed up the installation process, please refer to https://mirrors.tuna.tsinghua.edu.cn/help/pypi/ and https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/.

If the installation still fails, try the mamba/conda installation method above.

警告

On clusters, due to permission issues, the pip in the public path may not support global installation of Python libraries. You must add the --user option after pip install to install them in your home directory under ~/.local/lib/python3.x/site-packages/, where 3.x represents the Python interpreter version, and x can be any integer between 9 and 13.

Python will prioritize loading dspawpy from the home directory, even if the version in the public environment is newer! Therefore, if you have previously installed dspawpy with --user and have forgotten to manually update the old version in your home directory, even after sourcing the public environment, you will not be able to call the dspawpy in the public environment. Instead, the old version will still be used, leading to some bugs.

Therefore, considering that the HZW cluster automatically updates dspawpy weekly, it is recommended not to install it redundantly in your home directory; delete any existing installations. On other clusters, ensure that you manually update dspawpy in your home directory in a timely manner.

If you prefer not to delete and update the dspawpy in your home directory, you can use the -s option when running your Python scripts to prevent importing dspawpy from your home directory: python -s your-script.py.

Update dspawpy

To update dspawpy if it was installed with mamba/conda, use the following command:

mamba update dspawpy
#conda update dspawpy

If dspawpy was installed via pip:

pip install dspawpy -U # -U for upgrading to the latest version

备注

If pip uses a domestic mirror site, it may fail to upgrade smoothly because the mirror site has not yet synchronized the latest version of dspawpy . Please use the following command to tell pip to download and install from the official PyPI site:

pip install dspawpy -i https://pypi.org/simple --user -U # -i specifies the download address, --user installs for the current user only, and -U installs the latest version

If you encounter errors related to dspawpy during runtime, first verify that you have correctly imported the latest version of dspawpy and check the installation path:

$ python3 # or python

>>> import dspawpy
>>> dspawpy.__version__ # will output the version number
>>> dspawpy.__file__ # will output the installation path

8.2 structure structure conversion

To read structure information, use the read function; to write structure information to a file, use the write function; for quick structure conversion, use the convert function:

See the 2conversion.py script for conversion:

# coding:utf-8
from dspawpy.io.structure import convert

convert(
    infile="dspawpy_proj/dspawpy_tests/inputs/2.1/relax.h5",  # Structure to be converted, if in the current path, you can just write the filename
    si=None,  # Select configuration number, if not specified, read all by default
    ele=None,  # Filter element symbol, default reads atomic information for all elements
    ai=None,  # Filter atomic indices, starting from 1, default to read all atomic information
    infmt=None,  # Input structure file type, e.g., 'h5'. If None, it will be matched ambiguously based on the filename rule.
    task="relax",  # Task type, this parameter is only valid when infile is a folder rather than a filename
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/relaxed.xyz",  # Structure file name
    outfmt=None,  # Output structure file type, e.g., 'xyz'. If None, it will be fuzzy matched according to filename rules.
    coords_are_cartesian=True,  # Written in Cartesian coordinates by default
)

The rules for setting several key parameters of the convert function are shown in the table below:

**dspawpy Supported IO format**
infmt (input file format)	infile (Input file name fuzzy match)	outfmt (output file format)	outfile (output filename fuzzy match)	Description
h5	*.h5	X	X	HDF5 files saved after DS-PAW calculations are completed
json	*.json	json	*.json	json files saved after DS-PAW calculations are completed
pdb	*.pdb	pdb	*.pdb	Protein Data Bank
as	*.as	as	*.as	DS-PAW structure file containing atomic coordinates and other information
hzw	*.hzw	hzw	*.hzw	DeviceStudio's default structure file
xyz	*.xyz	xyz	*.xyz	Supports only single conformation of molecular structure when reading, and extended-xyz type trajectory files including unit cell when writing
X	X	dump	*.dump	LAMMPS dump-type trajectory files
X	.cif/.mcif	cif/mcif	.cif/.mcif	Crystallographic Information File
X	POSCAR/CONTCAR/.vasp/CHGCAR/LOCPOT/vasprun.xml*	poscar	POSCAR	VASP files
X	.cssr	cssr	.cssr	Crystal Structure Standard Representation
X	.yaml/.yml	yaml/yml	.yaml/.yml	YAML Ain't Markup Language
X	.xsf	xsf	.xsf	eXtended Structural Format
X	rndstr.in/lat.in/bestsqs	mcsqs	rndstr.in/lat.in/bestsqs	Monte Carlo Special Quasirandom Structure
X	inp.xml/.in/inp_	fleur-inpgen	.in	FLEUR structure file, requires the additional installation of the pymatgen-io-fleur library
X	*.res	res	*.res	ShelX res structure file
X	.config/.pwmat	pwmat	.config/.pwmat	PWmat files
X	X	prismatic	prismatic	A file format used for STEM simulations
X	CTRL*	X	X	Stuttgart LMTO-ASA files

备注

In the table above, * represents any character, and X indicates unsupported formats.
h5, json, pdb, xyz, dump, and CONTCAR formats support trajectory information consisting of multiple structures (common in structure optimization, NEB, or AIMD tasks)
The in(out)fmt parameter has higher priority than filename wildcard matching; for example, specifying in(out)fmt='h5' allows any filename, even a.json.
When writing structural information in json format, only visualization of NEB chain tasks is supported. See Observing the NEB Chain for details.
Structure information from DS-PAW output files such as neb.h5, phonon.h5, phonon.json, neb.json, and wannier.json is currently not readable.

8.3 Volumetric Data Processing

volumetricData Visualization

Differential volumetric data visualization

See 3dvol.py:

# coding:utf-8
from dspawpy.io.write import write_delta_rho_vesta

# Read the data file (h5 or json format), process it, and output it to a cube file, which can be directly opened with Vesta and has a small volume
write_delta_rho_vesta(
    total="dspawpy_proj/dspawpy_tests/inputs/supplement/AB.h5",  # Data file for the system containing all components
    individuals=[
        "dspawpy_proj/dspawpy_tests/inputs/supplement/A.h5",
        "dspawpy_proj/dspawpy_tests/inputs/supplement/B.h5",
    ],  # Data files for the system containing each component
    output="dspawpy_proj/dspawpy_tests/outputs/us/3delta_rho.cube",  # Output file path
)

The above script supports the processing of charge density differences in multi-component systems. As an example using a binary system, it generates the charge density difference file delta_rho.cube from AB.h5, A.h5, and B.h5. This file can be directly opened using VESTA.

Volumetric data plane average

See also 3planar_ave.py:

# coding:utf-8
from dspawpy.plot import average_along_axis

axes = [
    "2"
]  # "0", "1", "2" correspond to the x, y, z axes respectively; select which axes to average along
axes_indices = [int(i) for i in axes]
for ai in axes_indices:
    plt = average_along_axis(
        datafile="dspawpy_proj/dspawpy_tests/inputs/3.3/scf.h5",  # Data file path
        task="potential",  # Task name, can be 'rho', 'potential', 'elf', 'pcharge', 'rhoBound'
        axis=ai,  # Axis along which to plot the potential curve
        smooth=False,  # Whether to smooth
        smooth_frac=0.8,  # Smoothing coefficient
        subtype=None,  # Used to specify the subclass of task data, currently only used for Potential
        label=f"axis{ai}",  # Legend label
    )
if len(axes_indices) > 1:
    plt.legend()

plt.xlabel("Grid Index")
plt.ylabel("TotalElectrostaticPotential (eV)")
plt.savefig("dspawpy_proj/dspawpy_tests/outputs/us/3pot_ave.png", dpi=300)  # Image name

Processing the electrostatic potential file obtained from Section 3.3 of Application Cases yields the following vacuum direction potential function curve:

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's possible that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. You should either switch programs (e.g., VSCode or the system's built-in terminal command line) or add the following code, starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.4 band data processing

备注

The script calls get_band_data() to read data, and setting efermi=XX during data reading can shift the energy zero point to the specified value; setting zero_to_efermi=True can shift the energy zero point to the Fermi level in the read file.
When plotting using pymatgen's BSPlotter.get_plot() in the script, you can set zero_to_efermi=True to shift the energy zero to the Fermi level. Due to a critical update in pymatgen on August 17, 2023, which changed the return object of the plotting function from plt to axes, subsequent scripts may become incompatible. Therefore, a conditional statement has been added to handle this in the relevant parts of the user's script.
For two-step band calculations, obtain the accurate Fermi level from the first-step self-consistent calculation (from the self-consistent system.json). If this fails, users can modify the energy zero point when calling get_band_data to read data, using the efermi parameter. For example: band_data=get_band_data('band.h5',efermi=-1.5)
When plotting, the script calls BSPlotter.get_plot from pymatgen. When the system is determined to be non-metallic, setting zero_to_efermi will consider the VBM as the Fermi level energy, rather than the Fermi level from the data file. Therefore, when the system is non-metallic, setting zero_to_efermi=True during data reading and setting zero_to_efermi=True during plotting will result in different plots.

Running the Python script listed in this section, the program will determine whether the system is metallic. If it is a non-metallic system, you will be prompted to choose whether to shift the Fermi level to the zero energy point; please follow the prompts.

Conventional Band Treatment

See 4bandplot.py:

# coding:utf-8
import os
import matplotlib.pyplot as plt
from pymatgen.electronic_structure.plotter import BSPlotter

from dspawpy.io.read import get_band_data

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # Specifies the data file path
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # path to system.json file, required only when band_dir is a json file
    efermi=None,  # Used for manually correcting the Fermi level
    zero_to_efermi=True,  # For non-metallic systems, the zero point energy should be shifted to the Fermi level
)

bsp = BSPlotter(band_data)
axes_or_plt = bsp.get_plot(
    zero_to_efermi=False,  # The data has already been shifted when read, so this should be turned off
    ylim=[-10, 10],  # Range of the y-axis for the band structure plot
    smooth=False,  # Whether to smooth the band structure plot
    vbm_cbm_marker=False,  # Whether to mark the valence band maximum and conduction band minimum in the band structure plot
    smooth_tol=0,  # Threshold for smoothing
    smooth_k=3,  # Order of the smoothing process
    smooth_np=100,  # Number of points for smoothing
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# Add a reference line for the energy zero point
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot.png"  # Filename for the output band plot
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

备注

For band structure calculations, an accurate Fermi level is required, which is obtained from the self-consistent calculation (from system.json). If the acquisition fails, users can modify the efermi parameter in the get_band_data function.

Executing the code will generate a band structure plot similar to the following:

The band is projected onto each element separately, with the size of the data points representing the element's contribution to the orbital.

See 4bandplot_elt.py:

# coding:utf-8
import os

import matplotlib.pyplot as plt
import numpy as np
from pymatgen.electronic_structure.plotter import BSPlotterProjected

from dspawpy.io.read import get_band_data

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # Specify the data file path
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # path to system.json file, required only when band_dir is a json file
    efermi=None,  # Used to manually adjust the Fermi level
    zero_to_efermi=True,  # For non-metallic systems, shift the zero-point energy to the Fermi level
)

bsp = BSPlotterProjected(bs=band_data)  # Initialize the BSPlotterProjected class
axes_or_plt = bsp.get_elt_projected_plots(
    zero_to_efermi=False,  # The data has already been shifted when read, so this should be disabled
    ylim=[-8, 5],  # Set the energy range
    vbm_cbm_marker=False,  # Whether to mark the conduction band minimum (CBM) and valence band maximum (VBM)
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# Add a reference line for the energy zero point
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot_elt.png"  # The filename for the output band plot
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

备注

To plot projected band structure data, use the BSPlotterProjected module.
Use the get_elt_projected_plots function in the BSPlotterProjected module to plot band diagrams with orbital contributions for each element.

Executing the code will generate band plots similar to the following:

警告

If you execute the script above by connecting to a remote server via SSH, and you encounter QT-related error messages, it's possible that the program you're using (such as MobaXterm) is incompatible with the QT libraries. You can either switch to a different program (e.g., VSCode or the system's built-in terminal command line), or add the following code, starting on the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Band projections onto different elements' different orbitals

Refer to 4bandplot_elt_orbit.py:

# coding:utf-8
import os

import matplotlib.pyplot as plt
import numpy as np
from pymatgen.electronic_structure.plotter import BSPlotterProjected

from dspawpy.io.read import get_band_data

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # Specify the data file path
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # path to system.json file, only required when band_dir is a json file
    efermi=None,  # Used for manually correcting the Fermi level
    zero_to_efermi=True,  # For non-metallic systems, shift the zero point energy to the Fermi level
)

bsp = BSPlotterProjected(bs=band_data)  # Initialize the BSPlotterProjected class
# Select elements and orbitals, create a dictionary
dict_elem_orbit = {"Mo": ["d"], "S": ["s"]}

axes_or_plt = bsp.get_projected_plots_dots(
    dictio=dict_elem_orbit,
    zero_to_efermi=False,  # The data has already been shifted when read, so this should be turned off
    ylim=[-8, 5],  # Set the energy range
    vbm_cbm_marker=False,  # Whether to mark the conduction band minimum and valence band maximum
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# Add a reference line for the energy zero point
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandplot_elt_orbit.png"  # Filename for the output band plot
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

备注

Use the get_projected_plots_dots method in the BSPlotterProjected module, which allows users to customize the band structure plots by specifying elements and orbitals (L) to be plotted.
For example, get_projected_plots_dots({'Mo': ['d'], 'S': ['s']}) plots the d-orbitals of Mo and the s-orbitals of S.

Executing the code will generate a band structure plot similar to the following:

警告

If you execute the above script by connecting to a remote server via SSH, and QT-related error messages appear, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Projecting band structure onto different atomic orbitals

See 4bandplot_patom_porbit.py:

# coding:utf-8
import os

import matplotlib.pyplot as plt
import numpy as np
from pymatgen.electronic_structure.plotter import BSPlotterProjected

from dspawpy.io.read import get_band_data

datafile = "dspawpy_proj/dspawpy_tests/inputs/supplement/pband.h5"  # Specify the data file path
band_data = get_band_data(
    band_dir=datafile,
    syst_dir=None,  # path to system.json file, required only when band_dir is a JSON file
    efermi=None,  # Used to manually adjust the Fermi level
    zero_to_efermi=True,  # For non-metallic systems, shift the zero point energy to the Fermi level
)

bsp = BSPlotterProjected(bs=band_data)
# Specify elements, orbitals, and atomic numbers
dict_elem_orbit = {"Mo": ["px", "py", "pz"]}
dict_elem_index = {"Mo": [1]}

axes_or_plt = bsp.get_projected_plots_dots_patom_pmorb(
    dictio=dict_elem_orbit,  # Specify the element-orbit dictionary
    dictpa=dict_elem_index,  # Specify the element-atomic number dictionary
    sum_atoms=None,  # Whether to sum over atoms
    sum_morbs=None,  # Whether to sum orbitals
    zero_to_efermi=False,  # Data has already been shifted during reading, should be turned off here
    ylim=None,  # Set the energy range
    vbm_cbm_marker=False,  # Whether to mark the conduction band minimum and valence band maximum
    selected_branches=None,  # Specify the energy band branches to be plotted
    w_h_size=(12, 8),  # Set image width and height
    num_column=None,  # Number of images displayed per row
)

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# Add a reference line for the energy zero point
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = " dspawpy_proj/dspawpy_tests/outputs/us/4band_patom_porbit.png"  # Output bandpass figure filename
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

备注

The get_projected_plots_dots_patom_pmorb function in the BSPlotterProjected module offers greater flexibility, allowing users to customize the band diagrams for specific atoms and orbitals.
Use dictpa to specify the atom, and dictio to specify the orbitals of that atom.
To superimpose projected components of some atoms or orbitals, specify the sum_atoms or sum_morbs parameters according to the documentation of the get_projected_plots_dots_patom_pmorb function.

警告

If only a single orbital is selected and the orbital name has more than one letter (e.g., px, dxy, dxz), the get_projected_plots_dots_patom_pmorb function will raise an error. See here for details.

Executing the code will generate band diagrams similar to the following:

警告

If you execute the above script by connecting to a remote server via SSH, and you encounter QT-related error messages, it's possible that the program you're using (e.g., MobaXterm) is incompatible with the QT libraries. Either switch to another program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Band unfolding processing

See 4bandunfolding.py:

# coding:utf-8
import os

from dspawpy.plot import plot_bandunfolding

plt = plot_bandunfolding(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.22.1/band.h5",  # Read data
    ef=None,  # Fermi level, read from the file
    de=0.05,  # Band width, default 0.05
    dele=0.06,  # Band gap, default 0.06
)

plt.ylim(-15, 10)
figname = "dspawpy_proj/dspawpy_tests/outputs/us/4bandunfolding.png"  # Output band structure plot filename
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
plt.savefig(figname, dpi=300)
# plt.show()

Executing the code yields a band diagram similar to the following:

警告

This feature currently does not support setting the Fermi level of non-metallic materials as the zero-energy point (the default is the valence band top as the zero-energy point).

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it might be due to incompatibility between the program used (such as MobaXterm, etc.) and the QT library. You can either switch programs (e.g., VSCode or the system's built-in terminal) or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

band-compare band structure comparison figure processing

Plotting regular band structure and Wannier band structure on the same figure.

Refer to 4bandcompare.py:

# coding:utf-8
import os

from pymatgen.electronic_structure.plotter import BSPlotter

from dspawpy.io.read import get_band_data

band_data = get_band_data(
    band_dir="dspawpy_proj/dspawpy_tests/inputs/2.30/wannier.h5",  # Wannier band file path
    syst_dir=None,  # system.json file path, only needed when band_dir is a json file
    efermi=None,  # Used for manually adjusting the Fermi level
    zero_to_efermi=False,  # Whether to shift zero energy to the Fermi level
)
bsp = BSPlotter(bs=band_data)
band_data = get_band_data(
    band_dir="dspawpy_proj/dspawpy_tests/inputs/2.3/band.h5",  # Read DFT band structure
    syst_dir=None,  # path to system.json file, required only when band_dir is a json file
    efermi=None,  # Used for manually correcting the Fermi level
    zero_to_efermi=False,  # Whether to shift the zero point energy to the Fermi level
)

bsp2 = BSPlotter(bs=band_data)
bsp.add_bs(bsp2._bs)
axes_or_plt = bsp.get_plot(
    zero_to_efermi=True,  # Move the zero energy level to the Fermi level
    ylim=[-10, 10],  # Energy band plot y-axis range
    smooth=False,  # Whether to smooth the band structure plot
    vbm_cbm_marker=False,  # Whether to mark the valence band maximum and conduction band minimum in the band structure plot
    smooth_tol=0,  # Threshold for smoothing
    smooth_k=3,  # Order of the smoothing process
    smooth_np=100,  # Number of points for smoothing
    bs_labels=["wannier interpolated", "DFT"],  # Band structure labels
)

import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

# Add a reference line for the energy zero point
for ax in fig.axes:
    ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/4wanierBand.png"  # File name for the output band structure plot
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)
fig.savefig(figname, dpi=300)

Executing the code will generate band comparison curves similar to the following:

警告

If you are running the above script by connecting to a remote server via SSH and encounter QT-related error messages, it may be due to incompatibility between the program you are using (such as MobaXterm, etc.) and the QT libraries. You can either switch to another program (such as VSCode or the system's built-in terminal), or add the following code, starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.5 DOS Data Processing

Total Density of States

See 5dosplot_total.py:

# coding:utf-8
import os

from pymatgen.electronic_structure.plotter import DosPlotter

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # Read projected density of states data
    return_dos=False,  # If False, always return a CompleteDos object (regardless of whether projection was enabled during calculation)
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # Whether to set the Fermi level as the zero point
    stack=False,  # True indicates drawing an area chart
    sigma=None,  # Gaussian broadening, None indicates no smoothing process
)
dos_plotter.add_dos(
    label="total dos", dos=dos_data
)  # Set the legend for the density of states plot  # Pass the density of states data

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # Set the energy range
    ylim=[-15, 15],  # Set the density of states range
)
ax.axhline(0, lw=2, ls="-.", color="gray")

filename = "dspawpy_proj/dspawpy_tests/outputs/us/5dos_total.png"  # File name for the output density of states plot
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

备注

Use the get_dos_data function to convert the dos.h5 file obtained from DS-PAW calculations into a format supported by pymatgen.
Use the DosPlotter module to obtain the data from the DS-PAW calculated dos.h5 file.
The DosPlotter function can pass parameters: the stack parameter indicates whether to fill the DOS plots, and zero_at_efermi indicates whether to set the Fermi energy to zero in the DOS plot. Here, stack=False and zero_at_efermi=False are set.
Use add_dos in the DosPlotter module to add the DOS data.
Use the get_plot function in the DosPlotter module to plot the DOS.

Executing the code will generate a density of states plot similar to the following:

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely that the program you're using (such as MobaXterm, etc.) is incompatible with the QT libraries. You can either switch programs (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Project Density of States onto different orbitals

See 5dosplot_spd.py:

# coding:utf-8
import os

from pymatgen.electronic_structure.plotter import DosPlotter

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # Read projected DOS data
    return_dos=False,  # If False, always return a CompleteDos object (regardless of whether projection was enabled during calculation)
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # Whether to set the Fermi level as the zero point
    stack=False,  # True indicates drawing an area chart
    sigma=None,  # Gaussian broadening, None indicates no smoothing is applied
)
dos_plotter.add_dos_dict(
    dos_dict=dos_data.get_spd_dos(),
    key_sort_func=None,  # Orbital projection  # Specifies the sorting function
)
ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # Set the energy range
    ylim=None,  # Set the density of states range
)

ax.axhline(0, lw=2, ls="-.", color="gray")

filename = "dspawpy_proj/dspawpy_tests/outputs/us/5dos_spd.png"  # Filename of the output density of states plot
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(filename, dpi=300)

备注

Use the add_dos_dict function in the DosPlotter module to obtain the projected density of states (DOS) data, and then use get_spd_dos to project the information onto spd orbitals.

The code execution will produce a density of states plot similar to the following:

警告

If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it might be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Projecting the density of states onto different elements

Projecting the density of states onto different orbitals of different atoms

See 5dosplot_atom_orbit.py:

# coding:utf-8
import os

from pymatgen.electronic_structure.core import Orbital
from pymatgen.electronic_structure.plotter import DosPlotter

from dspawpy.io.read import get_dos_data
from dspawpy.plot import plot_dos

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/3.2.4/dos.h5",  # Reads projected density of states data
    return_dos=False,  # If False, always return a CompleteDos object (regardless of whether projection was enabled during calculation)
)
dos_plotter = DosPlotter(
    zero_at_efermi=True,  # Whether to set the Fermi level as the zero point
    stack=False,  # True indicates drawing an area plot
    sigma=None,  # Gaussian broadening, None indicates no smoothing treatment
)

# ! Specify atomic number and orbital
dict_index_orbit = {0: ["dxy"], 2: ["s"]}

print("Plotting...")
for index in dict_index_orbit:
    _os = dict_index_orbit[index]
    _e = str(dos_data.structure.sites[index].species)
    for _orb in _os:
        dos_plotter.add_dos(
            f"{_e}(atom-{index}) {_orb}",  # label
            dos_data.get_site_orbital_dos(
                dos_data.structure[index],
                getattr(Orbital, _orb),
            ),
        )

ax = plot_dos(
    dosplotter=dos_plotter,
    xlim=[-10, 5],  # Set the energy range
    ylim=None,  # Set the density of states range
)
ax.axhline(0, lw=2, ls="-.", color="gray")

figname = "dspawpy_proj/dspawpy_tests/outputs/us/5dos_atom_orbit.png"  # Output density of states figure filename
os.makedirs(os.path.dirname(os.path.abspath(figname)), exist_ok=True)

fig = ax.get_figure()
fig.savefig(figname, dpi=300)

备注

Use the get_site_orbital_dos function to extract the contribution of a specific atom and specific orbital from the DOS data. dos_data.structure[0], Orbital(4) represents obtaining the density of states for the dxy orbital of the first atom; the index in the get_site_orbital_dos function starts from 0.
Running this script and selecting the element and orbital as prompted will generate the corresponding density of states (DOS) plot.

Executing the code will produce a density of states plot similar to the following:

警告

If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it's likely due to incompatibility between the program you're using (e.g., MobaXterm) and the QT library. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code to your Python script starting from the second line:

import matplotlib
matplotlib.use('agg')

Projecting the density of states onto the split d-orbitals (t2g, eg) of different atoms

d-centered analysis

Taking the Pb-slab system as an example, a d-band center analysis is performed on Pt atoms:

See 5center_dband.py:

# coding:utf-8
from dspawpy.io.read import get_dos_data
from dspawpy.io.utils import d_band

dos_data = get_dos_data(
    dos_dir="dspawpy_proj/dspawpy_tests/inputs/supplement/dos.h5",  # Read projected density of states data
    return_dos=False,  # If False, always returns a CompleteDos object (regardless of whether projection was enabled during calculation)
)
for spin in dos_data.densities:
    print("spin=", spin)
    c = d_band(spin, dos_data)
    print(c)

Executing the code yields results similar to the following:

spin=1
-1.785319344084034

备注

Currently, only the d-orbital center averaged over all atoms is supported. Element-resolved, atom-projected, or other orbitals are not supported, nor is the selection of spin direction or energy range.

The get_dos_data function is responsible for processing density of states data:

8.6 bandDos: Displaying Band Structure and Density of States Together

Using the Si system from the application tutorial as an example:

Display band structure and density of states in a single figure.

See 6bandDosplot.py:

# coding:utf-8
import os

import numpy as np
from matplotlib.axes import Axes
from pymatgen.electronic_structure.plotter import BSDOSPlotter

from dspawpy.io.read import get_band_data, get_dos_data

bandfile = "dspawpy_proj/dspawpy_tests/inputs/2.3/band.h5"  # Normal band data
band_data = get_band_data(
    band_dir=bandfile,
    syst_dir=None,  # path to system.json file, required only when band_dir is a json file
    efermi=None,  # Used for manually correcting the Fermi level
)
band_efermi = band_data.efermi
dosfile = "dspawpy_proj/dspawpy_tests/inputs/2.5/dos.h5"  # Density of states data
dos_data = get_dos_data(
    dos_dir=dosfile,
    return_dos=False,  # If False, always return a CompleteDos object (regardless of whether projection was enabled during calculation)
)
dos_efermi = dos_data.efermi
bdp = BSDOSPlotter(
    bs_projection=None,  # Band structure projection method, None means no projection
    dos_projection=None,  # Projection method for density of states, None means no projection
    vb_energy_range=4,  # Valence band energy range
    cb_energy_range=4,  # Conduction band energy range
    fixed_cb_energy=False,  # Whether to fix the conduction band energy range
    egrid_interval=1,  # Energy grid interval
    font="DejaVu Sans",  # Default is Times New Roman, change to DejaVu Sans to avoid warnings due to missing font on Linux
    axis_fontsize=20,  # Axis font size
    tick_fontsize=15,  # Tick label font size
    legend_fontsize=14,  # Legend font size
    bs_legend="best",  # Band structure legend position
    dos_legend="best",  # Density of States legend position
    rgb_legend=True,  # Use colored legend
    fig_size=(11, 8.5),  # Figure size
)
if band_efermi != dos_efermi:
    print(f"{band_efermi=:.4f} eV")
    print(f"{dos_efermi=:.4f} eV")
    d_efermi = band_efermi - dos_efermi

    print(
        "! Band and DOS Fermi levels are inconsistent, using DOS Fermi level as reference"
    )
    band_data.bands = {spin: v + d_efermi for spin, v in band_data.bands.items()}

    # ! Band and DOS Fermi levels are inconsistent, using Band level as the reference
    # dos_data.energies -= d_efermi

axes_or_plt = bdp.get_plot(
    bs=band_data, dos=dos_data
)  # Pass band data  # Pass density of states data

if isinstance(axes_or_plt, Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = "dspawpy_proj/dspawpy_tests/outputs/us/6bandDos.png"  # Filename for the band structure - density of states plot output
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)
print("==> Saved", filename)

Executing the code yields a band density of states plot similar to the following:

警告

If you are connecting to a remote server via SSH and running the above script, and you encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT libraries. You should either switch programs (e.g., VSCode or the system's built-in terminal command line) or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Display band structure and projected density of states on a single plot.

See 6bandPdosplot.py:

# coding:utf-8
import os

import numpy as np
from matplotlib.axes import Axes
from pymatgen.electronic_structure.plotter import BSDOSPlotter

from dspawpy.io.read import get_band_data, get_dos_data

bandfile = "dspawpy_proj/dspawpy_tests/inputs/2.4/band.h5"  # Normal band data
band_data = get_band_data(
    band_dir=bandfile,
    syst_dir=None,  # path to system.json file, required only when band_dir is a json file
    efermi=None,  # Used for manually correcting the Fermi level
)
band_efermi = band_data.efermi
dosfile = (
    "dspawpy_proj/dspawpy_tests/inputs/2.6/dos.h5"  # DOS data for projected states
)
dos_data = get_dos_data(
    dos_dir=dosfile,
    return_dos=False,  # If False, always return a CompleteDos object (regardless of whether projection was enabled during calculation)
)
dos_efermi = dos_data.efermi
bdp = BSDOSPlotter(
    bs_projection="elements",  # Projection method for band structure, None means no projection
    dos_projection="elements",  # Project DOS onto elements
    vb_energy_range=4,  # Valence band energy range
    cb_energy_range=4,  # Conduction band energy range
    fixed_cb_energy=False,  # Whether to fix the conduction band energy range
    egrid_interval=1,  # Energy grid interval
    font="DejaVu Sans",  # Default is Times New Roman, can be changed to DejaVu Sans to avoid warnings due to font not being installed on Linux
    axis_fontsize=20,  # Axis font size
    tick_fontsize=15,  # Tick label font size
    legend_fontsize=14,  # Legend font size
    bs_legend="best",  # Band structure legend position
    dos_legend="best",  # Position of the projected density of states legend
    rgb_legend=True,  # Use colored legend
    fig_size=(11, 8.5),  # Figure size
)
if band_efermi != dos_efermi:
    print(f"{band_efermi=:.4f} eV")
    print(f"{dos_efermi=:.4f} eV")
    d_efermi = band_efermi - dos_efermi

    print(
        "! Band and DOS Fermi levels are inconsistent, using DOS Fermi level as reference"
    )
    band_data.bands = {spin: v + d_efermi for spin, v in band_data.bands.items()}

    # ! Band and DOS Fermi levels are inconsistent, using Band level as reference
    # dos_data.energies -= d_efermi

axes_or_plt = bdp.get_plot(
    bs=band_data,
    dos=dos_data,
)  # Pass band structure data  # Pass projected density of states data

if isinstance(axes_or_plt, Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif np.iterable(axes_or_plt):
    fig = np.asarray(axes_or_plt).flatten()[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = "dspawpy_proj/dspawpy_tests/outputs/us/6bandPdos.png"  # filename for the band structure-projected density of states plot
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)
print("==> Saved", filename)

Executing the code yields a band-decomposed density of states plot similar to the following:

警告

Given projected band data, it will be projected along the element by default; given ordinary band data (or if the system contains more than 4 types of elements), it will not be projected and a warning will be output.
Given projected density of states (PDOS) data, projection along elements is also the default. You can switch to projection along orbitals, or no projection at all. For ordinary density of states (DOS) data and without disabling the DOS projection option BSDOSPlotter(dos_projection=None), the pymatgen plotting program will report an error, which is why a 6bandDosplot.py file was specifically prepared, as mentioned above.

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely due to incompatibility between the program you are using (e.g., MobaXterm) and the QT library. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.7 optical data processing

Using the scf.h5 file obtained from a quick start calculation of the optical properties of the Si system as an example (Note: the output file name is the same as the task, task = scf; io.optical = true can calculate optical properties):

Processing the reflectivity data, referring to 7optical.py:

# coding:utf-8
from dspawpy.plot import plot_optical

plot_optical(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.12/scf.h5",
    keys=["ExtinctionCoefficient", "Reflectance"],
    axes=["X"],  # ["X", "Y", "Z", "XY", "YZ", "ZX"]
    prefix="dspawpy_proj/dspawpy_tests/outputs/optical",  # Where to save, if empty, it means the current folder
    save=True,  # Whether to save the image with the tool's name, if False, please refer to the script below to save manually
)

# The above function will plot and save the images of ExtinctionCoefficient and Reflectance separately
# To plot multiple properties on the same figure, uncomment the following code and set the save parameter above to False

# import os
# import matplotlib.pyplot as plt
#
# plt.tick_params(labelsize=16)
# plt.tight_layout()
# filename = "outputs/us/7optical.png"  # Filename for the output optical properties plot
# os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
# plt.savefig(filename, dpi=300)

备注

Reflectance is an optical property, and users can modify this keyword to "AbsorptionCoefficient", "ExtinctionCoefficient", or "RefractiveIndex" based on their needs, corresponding to the absorption coefficient, extinction coefficient, and refractive index, respectively.

Executing the code will generate a curve showing the reflectance as a function of energy, similar to the following:

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it is likely that the program you are using (such as MobaXterm, etc.) is incompatible with the QT libraries. You can either switch to another program (such as VSCode or the system's built-in terminal command line) or add the following code to your Python script, starting from the second line:

import matplotlib
matplotlib.use('agg')

8.8 neb data processing

Let's start with a quick introduction using the H diffusion on Pt(100) surface example:

Generating intermediate configurations for input files

See 8neb_interpolate_structures.py:

# coding:utf-8
from dspawpy.diffusion.neb import NEB, write_neb_structures
from dspawpy.diffusion.nebtools import write_json_chain
from dspawpy.io.structure import read

# Read initial configuration
init_struct = read("dspawpy_proj/dspawpy_tests/inputs/2.15/00/structure00.as")[0]
# Read final state configuration
final_struct = read("dspawpy_proj/dspawpy_tests/inputs/2.15/04/structure04.as")[0]

neb = NEB(
    initial_structure=init_struct,  # Initial structure
    final_structure=final_struct,  # Final state configuration
    nimages=8,  # Total of 8 configurations, including initial and final states
)
structures = neb.linear_interpolate()  # Linear interpolation
# structures = neb.idpp_interpolate()   # IDPP interpolation

# Save as structure file to dest path
write_neb_structures(
    structures=structures,  # Insert interpolated structure chains
    coords_are_cartesian=True,  # Whether to save in Cartesian coordinates
    fmt="as",  # Save format, supported formats: 'json', 'as', 'hzw', 'pdb', 'xyz', 'dump'
    path="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures",  # Save path
    prefix="structure",  # File name prefix
)

# Preview initial structure chain
write_json_chain(
    preview=True,  # whether to enable preview mode
    directory="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures",  # Directory for NEB calculations
    step=-1,  # Default to saving the structure chain of the last ion step (latest)
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # Save path
)
# write_xyz_chain(preview=True, # Whether to run in preview mode
#                  directory="dspawpy_proj/dspawpy_tests/outputs/us/8neb_interpolate_structures", # NEB calculation directory
#                  step=-1, # Default to saving the structure chain of the last ionic step (latest)
#                  dst='dspawpy_proj/dspawpy_tests/outputs/us/8neb' # save path
# )

备注

Users can modify the number of interpolated points as needed. Setting it to 8 will generate a folder containing 8 structure files, with 6 intermediate configurations.
neb.linear_interpolate is a linear interpolation method. The pbc parameter, when set to True, will lock the search for the shortest diffusion path. It defaults to False to increase user control, because
For example, if the initial fractional coordinate of an atom is 0.2 and the final state is 0.8. When pbc = True, the diffusion path will be forced to be 0.2 -> -0.2. When pbc = False, the user can make the program perform interpolation along the diffusion path 0.2 -> 0.8; if the shortest path is desired, manually change 0.8 to -0.2, thereby ensuring the program completes the initial guess of interpolation according to the user's intent.

Plotting the energy barrier diagram

neb.iniFin = true/false

When neb.iniFin = true/false, you can use the path from the NEB calculation for barrier analysis (ensure that the initial and final state calculation files are in the NEB calculation path):

Refer to 8neb_barrier_CubicSpline.py:

# coding:utf-8
from dspawpy.diffusion.nebtools import plot_barrier

directory_of_neb_task = (
    "dspawpy_proj/dspawpy_tests/inputs/2.15"  # <-- Please modify to the actual NEB path
)

# Plotting the energy barrier using CubicSpline interpolation
plot_barrier(
    directory=directory_of_neb_task,  # path of the neb task
    method="CubicSpline",  # Cubic spline interpolation
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_CubicSpline.png",  # Output filename for the energy barrier plot
    show=False,  # Whether to display the energy barrier plot
)

备注

After running the above script, you can obtain a barrier curve similar to the following, with cubic spline interpolation:

For this specific example, the curve will exhibit an undesirable "dip" after cubic spline interpolation, which is inherent to the characteristics of the cubic spline interpolation algorithm.

dspawpy internally calls scipy's interpolation algorithms. Taking the cubic spline interpolation algorithm as an example in the script above, it is defined in the scipy documentation as:

class scipy.interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)

The keyword arguments include axis, bc_type, and extrapolate, whose specific meanings can be found in scipy.interpolate.CubicSpline. We can specify the corresponding keyword arguments (axis, bc_type, extrapolate) in the plot_barrier function and pass them to the scipy.interpolate.CubicSpline class for processing.

Here we use the script 8neb_barrier.py to compare the curves plotted by interpolating with three algorithms:

# coding:utf-8
import os

import matplotlib.pyplot as plt

from dspawpy.diffusion.nebtools import plot_barrier

# Compare the differences in energy barrier curves drawn by different interpolation methods, where show should be set to False
# 1. interp1d
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # path for NEB calculation
    ri=None,  # Reaction coordinate between the initial structure and the second structure, required when the NEB task only calculated intermediate structures
    rf=None,  # Reaction coordinate between the last configuration and the second-to-last configuration, when the NEB task only calculated intermediate configurations
    ei=None,  # Energy of the initial configuration, required when the NEB task only calculated intermediate configurations
    ef=None,  # Energy of the final configuration, required when the NEB task only calculated intermediate configurations
    method="interp1d",  # Interpolation method
    figname=None,  # Name of the output energy barrier plot file
    show=False,  # Whether to display the energy barrier plot
    kind="quadratic",  # Parameter of the interpolation method
)
# 2. CubicSpline
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    method="CubicSpline",
    figname=None,
    show=False,
)
# 3. pchip
plot_barrier(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    method="pchip",
    figname=None,
    show=False,
)

filename = "dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_comparison.png"  # Filename for the energy barrier plot output
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
plt.savefig(filename, dpi=300)
# plt.show()

备注

Choosing the appropriate interpolation algorithm is crucial for optimizing the final curve presentation.
In most cases, selecting the pchip (piecewise cubic Hermite interpolating polynomial) monotonic cubic spline interpolation algorithm will achieve good results, and it is also the default interpolation algorithm called.

neb.iniFin = true

When neb.iniFin = true is set, reading the neb.h5/neb.json files generated by the NEB calculation allows for a quick barrier analysis:

See 8neb_barrier_CubicSpline.py:

# coding:utf-8
from dspawpy.diffusion.nebtools import plot_barrier

# Plot energy barrier using CubicSpline interpolation
plot_barrier(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.15/neb.h5",  # Path to neb.h5
    method="CubicSpline",  # Cubic spline interpolation
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb_barrier_.png",  # Output file name for the energy barrier plot
    show=False,  # Whether to display the energy barrier plot
)

Processing the resulting barrier diagram is consistent with the previously read path.

备注

The energy stored in neb.h5 and neb.json files is TotalEnergy. If you need an accurate barrier value, it is recommended to process it by reading the NEB calculation path (taking TotalEnergy0).

警告

If you are connecting to a remote server via SSH to execute the script above and encounter QT-related error messages, it's possible that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. Either switch to a different program (such as VSCode or the system's built-in terminal command line), or add the following code starting on the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Processing Data for Transition State Calculations

After NEB calculations, it is generally necessary to plot the energy barrier diagram and check the forces on each interpolated structure to ensure they are below a specified threshold. If the results are abnormal, the force and energy changes of each interpolated structure during the structure optimization process should also be checked to determine if they have truly "converged." These operations require at least three cycles. To simplify the process, we provide an all-in-one summary function summary:

Refer to 8neb_check_results.py:

# coding:utf-8
from dspawpy.diffusion.nebtools import summary

# Import the neb calculation directory, a complete folder after neb calculation needs to be provided
summary(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",
    show_converge=False,  # Whether to display the convergence plots of energy and force
    outdir="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # Path to save convergence plots of energy and forces
    figname="dspawpy_proj/dspawpy_tests/outputs/us/8neb/neb_barrier_summary.png",  # Path to save the energy barrier plot
)
# Additional keyword arguments can be set for plotting the barrier diagram, such as:
# summary(directory='dspawpy_proj/dspawpy_tests/inputs/2.15', method='CubicSpline') # Change to CubicSpline for spline interpolation

备注

This script will print the energies and forces of each structure in a table, plot the energy barrier, and also plot the convergence of energy and forces for intermediate structures.
If neb.iniFin = false, the user must copy the results file of the self-consistent calculation, either scf.h5 or system.json, to the corresponding initial and final state subfolders. Otherwise, the program cannot read the energy and force information of the initial and final states and will exit with an error.
By default, the energy barrier plot is stored in the parent directory of the NEB calculation, and the energy and force convergence plots for each intermediate structure are stored in the respective subfolders.

Executing the code will generate a table similar to the following, displaying the energy and force information for each NEB configuration:

Image	Force (eV/Å)	Reaction coordinate (Å)	Energy (eV)	Delta energy (eV)
00	0.1803	0.0000	-39637.0984	0.0000
01	0.0263	0.5428	-39637.0186	0.0798
02	0.0248	1.0868	-39636.8801	0.2183
03	0.2344	1.5884	-39636.9984	0.1000
04	0.0141	2.0892	-39637.0900	0.0084

In addition to the energy barrier diagram, you can also obtain the energy and force convergence curves for each intermediate configuration (taking configuration 02 as an example).

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it might be due to incompatibility between the program you are using (such as MobaXterm) and the QT library. Either switch to another program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Observing the NEB Chain

Here, the NEB chain refers to the geometric relationship between the interpolated structures (structure00.as, structure01.as, ...), rather than the changes of a single structure during the optimization process.

NEB calculations are computationally expensive, and observing the NEB chain helps to judge the convergence speed of the NEB calculation. Furthermore, after generating intermediate structures via interpolation, previewing the NEB chain is often necessary. These needs can be met using the 8neb_visualize.py script:

# coding:utf-8
from dspawpy.diffusion.nebtools import write_json_chain, write_xyz_chain

# Convert the configuration chain under the NEB calculation path to a JSON format file
write_json_chain(
    preview=False,  # If the NEB calculation is already completed, preview mode is not required
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # NEB calculation directory
    step=-1,  # Default to saving the configuration chain of the last ion step (latest)
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # Save path
    ignorels=False,  # Set to True to ignore latestStructureXX.as files
)

# Convert the configuration chain in the NEB calculation path to xyz format files
write_xyz_chain(
    preview=False,  # If the NEB calculation is already completed, preview mode is not required
    directory="dspawpy_proj/dspawpy_tests/inputs/2.15",  # NEB calculation directory
    step=-1,  # Default to saving the configuration chain of the last ionic step (latest)
    dst="dspawpy_proj/dspawpy_tests/outputs/us/8neb",  # Save path
    ignorels=False,  # Set to True to ignore latestStructureXX.as files
)

备注

After this script generates the neb_movie*.json files, you can view them by opening the json file via Device Studio --> Simulator --> DS-PAW --> Analysis Plot.
The directory parameter specifies the main path of the NEB calculation; the complete folder after the NEB calculation is finished must be provided.
This script supports processing ongoing (i.e., incomplete) NEB calculation files, allowing users to monitor the trajectory in real time.
The xyz file can be opened and viewed using OVITO software: Open the visualization interface via Device Studio --> Simulator --> OVITO, and then drag and drop the xyz file.
Structure information reading priority: latestStructureXX.as > h5 > json; When ignorels is set to True, it first attempts to read data from h5, and if it fails, it reads from json.

Calculate the inter-configuration distance

Refer to this script: 8calc_dist.py:

# coding:utf-8
from dspawpy.diffusion.nebtools import get_distance
from dspawpy.io.structure import read

# Please modify the paths of structure01.as and structure02.as structure files according to the actual situation
# First read the fractional coordinates, element list, and cell information of the two configurations
s1 = read("dspawpy_proj/dspawpy_tests/inputs/2.15/01/structure01.as")[0]
s2 = read("dspawpy_proj/dspawpy_tests/inputs/2.15/02/structure02.as")[0]
# Calculate the distance between the two configurations, note that this function only accepts fractional coordinates
dist = get_distance(
    spo1=s1.frac_coords,
    spo2=s2.frac_coords,
    lat1=s1.lattice.matrix,
    lat2=s2.lattice.matrix,
)
print("The distance between the two configurations is:", dist, "Angstrom")

Continued calculation with neb

To restart a NEB calculation, refer to 8neb_restart.py:

# coding:utf-8
import os
from shutil import copytree, rmtree

from dspawpy.diffusion.nebtools import restart

if os.path.isdir("dspawpy_proj/dspawpy_tests/outputs/us/neb4bk"):
    rmtree("dspawpy_proj/dspawpy_tests/outputs/us/neb4bk")

copytree(
    "dspawpy_proj/dspawpy_tests/inputs/2.15",
    "dspawpy_proj/dspawpy_tests/outputs/us/neb4bk",
)
restart(
    directory="dspawpy_proj/dspawpy_tests/outputs/us/neb4bk",  # NEB task path
    output="dspawpy_proj/dspawpy_tests/outputs/us/8neb_restart",  # Backup destination
)

See Continued calculation with neb for details.

Energy and maximum atomic force variation trend during NEB calculation

To view plots showing the energy and maximum atomic force trends during the NEB calculation, refer to 8neb_energy_force_curves.py:

# coding:utf-8
from dspawpy.diffusion.nebtools import monitor_force_energy

# Specify the path to the NEB calculation folder; after running, Energies.png and MaxForce.png images will be generated in the specified directory
unfinished_neb_folder = "dspawpy_proj/dspawpy_tests/inputs/supplement/neb_unfinished"
monitor_force_energy(
    directory=unfinished_neb_folder,
    outdir="imgs",  # Output image path
)

Generates energy and force change trend charts:

8.9 Phonon Data Processing

Using the example of a phonon band structure and density of states calculation for MgO, using phonon.h5:

If phonopy is not installed, running the following script will result in the message no module named 'phonopy', but this does not affect the program's normal operation.

Phonon band data processing

Refer to 9phonon_bandplot.py:

# coding:utf-8
import os

from pymatgen.phonon.plotter import PhononBSPlotter

from dspawpy.io.read import get_phonon_band_data

band_data = get_phonon_band_data(
    "dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5",
)  # Read phonon band structure
bsp = PhononBSPlotter(band_data)
axes_or_plt = bsp.get_plot(ylim=None, units="thz")  # Y-axis range # Units
import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif isinstance(axes_or_plt, tuple):
    fig = axes_or_plt[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = "dspawpy_proj/dspawpy_tests/outputs/us/9phonon_bandplot.png"  # File name for the output phonon band plot
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

Executing the code yields a phonon band structure curve similar to the following:

警告

If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it's likely due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code to your Python script starting from the second line:

import matplotlib
matplotlib.use('agg')

Phonon Density of States Data Processing

Refer to 9phonon_dosplot.py:

# coding:utf-8
import os

from pymatgen.phonon.plotter import PhononDosPlotter

from dspawpy.io.read import get_phonon_dos_data

dos = get_phonon_dos_data("dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5")
dp = PhononDosPlotter(
    stack=False,  # True indicates drawing an area plot
    sigma=None,  # Gaussian blur parameter
)
dp.add_dos(
    label="Phonon", dos=dos
)  # Legend  # The phonon density of states to be plotted
axes_or_plt = dp.get_plot(
    xlim=[0, 20],  # x-axis range
    ylim=None,  # y-axis range
    units="THz",  # Unit
)
import matplotlib.pyplot as plt  # noqa: E402

if isinstance(axes_or_plt, plt.Axes):
    fig = axes_or_plt.get_figure()  # version newer than v2023.8.10
elif isinstance(axes_or_plt, tuple):
    fig = axes_or_plt[0].get_figure()
else:
    fig = axes_or_plt.gcf()  # older version pymatgen

filename = " dspawpy_proj/dspawpy_tests/outputs/us/9phonon_dosplot.png"  # Energy barrier plot output filename
os.makedirs(os.path.dirname(os.path.abspath(filename)), exist_ok=True)
fig.savefig(filename, dpi=300)

Executing the code yields a phonon density of states curve similar to the following:

警告

If you execute the script above by SSH connection to a remote server and encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Phonon Thermodynamic Data Processing

Refer to 9phonon_thermal.py:

# coding:utf-8
from dspawpy.plot import plot_phonon_thermal

plot_phonon_thermal(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.h5",  # phonon.h5 data file path
    figname="dspawpy_proj/dspawpy_tests/outputs/us/9phonon.png",  # Output phonon thermodynamics figure filename
    show=False,  # Whether to display the image
)

Executing the code yields phonon thermodynamic curves similar to the following:

API: get_phonon_band_data(), get_phonon_dos_data(), plot_phonon_thermal()

The get_phonon_band_data function is responsible for reading phonon band data:
dspawpy.io.read.get_phonon_band_data(phonon_band_dir: str, verbose: bool = False)
Reads phonon band data from an h5 or json file and constructs a PhononBandStructureSymmLine object

参数:

phonon_band_dir -- Path to the band structure file, phonon.h5 / phonon.json, or a folder containing these files

返回类型:

PhononBandStructureSymmLine

示例
>>> from dspawpy.io.read import get_phonon_band_data >>> band_data = get_phonon_band_data("dspawpy_proj/dspawpy_tests/inputs/2.16/phonon.h5") # Read phonon band data >>> band_data = get_phonon_band_data("dspawpy_proj/dspawpy_tests/inputs/2.16/phonon.json") # Read phonon band data

The get_phonon_dos_data function is responsible for reading the phonon density of states:

dspawpy.io.read.get_phonon_dos_data(phonon_dos_dir: str, verbose: bool = False)

Reads phonon density of states data from an h5 or json file, constructs a PhononDos object

参数:: phonon_dos_dir -- Path to the phonon DOS file, phonon_dos.h5 / phonon_dos.json, or a folder containing these files
返回类型:: PhononDos

示例

>>> from dspawpy.io.read import get_phonon_dos_data
>>> phdos = get_phonon_dos_data(phonon_dos_dir='dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.json')
>>> phdos = get_phonon_dos_data(phonon_dos_dir='dspawpy_proj/dspawpy_tests/inputs/2.16.1/phonon.h5')
>>> phdos.frequencies
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,
        1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.1,
        2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3. ,  3.1,  3.2,
        3.3,  3.4,  3.5,  3.6,  3.7,  3.8,  3.9,  4. ,  4.1,  4.2,  4.3,
        4.4,  4.5,  4.6,  4.7,  4.8,  4.9,  5. ,  5.1,  5.2,  5.3,  5.4,
        5.5,  5.6,  5.7,  5.8,  5.9,  6. ,  6.1,  6.2,  6.3,  6.4,  6.5,
        6.6,  6.7,  6.8,  6.9,  7. ,  7.1,  7.2,  7.3,  7.4,  7.5,  7.6,
        7.7,  7.8,  7.9,  8. ,  8.1,  8.2,  8.3,  8.4,  8.5,  8.6,  8.7,
        8.8,  8.9,  9. ,  9.1,  9.2,  9.3,  9.4,  9.5,  9.6,  9.7,  9.8,
        9.9, 10. , 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9,
       11. , 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12. ,
       12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13. , 13.1,
       13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14. , 14.1, 14.2,
       14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15. , 15.1, 15.2, 15.3,
       15.4, 15.5, 15.6, 15.7, 15.8, 15.9, 16. , 16.1, 16.2, 16.3, 16.4,
       16.5, 16.6, 16.7, 16.8, 16.9, 17. , 17.1, 17.2, 17.3, 17.4, 17.5,
       17.6, 17.7, 17.8, 17.9, 18. , 18.1, 18.2, 18.3, 18.4, 18.5, 18.6,
       18.7, 18.8, 18.9, 19. , 19.1, 19.2, 19.3, 19.4, 19.5, 19.6, 19.7,
       19.8, 19.9, 20. ])

The plot_phonon_thermal function is responsible for plotting phonon thermodynamic properties:
dspawpy.plot.plot_phonon_thermal(datafile: str = 'phonon.h5', figname: str = 'phonon.png', show: bool = True, raw: bool = False, verbose: bool = False)
Task completed for phonon thermodynamic calculations, plot curves of relevant physical quantities versus temperature

phonon.h5/phonon.json -> phonon.png
参数:
datafile -- Path to an h5 or json file or a folder containing any of these files, default 'phonon.h5'

figname -- Filename to save the image

show -- Whether to pop up an interactive interface

raw -- Whether to save the plotting data to rawphonon.csv file
返回:

Image path, default 'phonon.png'

返回类型:

figname
示例
>>> from dspawpy.plot import plot_phonon_thermal >>> plot_phonon_thermal('dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.h5', figname='dspawpy_proj/dspawpy_tests/outputs/doctest/phonon_thermal_h5.png', show=False) >>> plot_phonon_thermal('dspawpy_proj/dspawpy_tests/inputs/2.26/phonon.json', figname='dspawpy_proj/dspawpy_tests/outputs/doctest/phonon_thermal_json.png', show=False, raw=True)

警告

If you execute the above script by connecting to a remote server via SSH and encounter QT-related error messages, it's likely that the program you are using (e.g., MobaXterm) is incompatible with the QT libraries. You can either switch programs (e.g., VSCode or the system's built-in terminal) or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.10 aimd molecular dynamics simulation data processing

For a quick start, take the molecular dynamics simulation data of the \(H_{2}O\) molecular system, for example, the aimd.h5 file:

Convert the trajectory file format to .xyz or .dump.

Read data from the HDF5 file output by AIMD and generate trajectory files.

The generated .xyz or .dump files can be dragged and dropped into OVITO for visualization. You can open the OVITO visualization interface through Device Studio --> Simulator --> OVITO, and then drag and drop the .xyz or .dump files into OVITO.

See 10write_aimd_traj.py:

# coding:utf-8
from dspawpy.io.structure import convert

convert(
    infile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # Structure to be converted, if in the current path, only the filename can be written
    si=None,  # Filter the configuration number, if not specified, read all by default
    ele=None,  # Filter element symbol, default reads atomic information for all elements
    ai=None,  # Filter atomic indices (starting from 1), default to read all atomic information
    outfile="dspawpy_proj/dspawpy_tests/outputs/us/10aimdTraj.xyz",  # Can also generate .dump files (lower precision), currently only supports orthogonal unit cells
)

Executing the code will generate trajectory files in .xyz and .dump formats, which can be opened with OVITO. For more details on structure file conversion, refer to structure structure conversion

备注

OVITO and dspawpy do not support saving systems with non-orthogonal unit cells as dump files.

Changes in energy, temperature, etc. curves during the dynamics process.

Refer to 10check_aimd_conv.py:

# coding:utf-8
from dspawpy.plot import plot_aimd

plot_aimd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # Data file path
    show=False,  # Whether to pop up the image window
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10aimd.png",  # Output image file name
    flags_str="1 2 3 4 5",  # Select data types
)
# The meaning of flags_str is as follows
# 1. Kinetic energy
# 2. Total Energy
# 3. Pressure
# 4. Temperature
# 5. Volume

Executing the code will generate the following combined graph:

警告

If you execute the above script by SSH connection to a remote server and encounter QT-related error messages, it's possible that the program you are using (such as MobaXterm) is incompatible with the QT libraries. You can either switch programs (for example, VSCode or the system's built-in terminal command line), or add the following code starting from the second line of the Python script:

import matplotlib
matplotlib.use('agg')

Mean Squared Displacement (MSD) Analysis

See 10aimd_msd.py:

# coding:utf-8
from dspawpy.analysis.aimdtools import get_lagtime_msd, plot_msd

# If AIMD is not completed in one go, you can assign multiple h5 file paths to the datafile parameter in a list form
lagtime, msd = get_lagtime_msd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # Data file path
    select="all",  # Default selects all atoms
    msd_type="xyz",  # Default to calculate MSD in the xyz directions
    timestep=None,  # Default reads the timestep from the datafile
)
# Plot the graph using the obtained data and save it
plot_msd(
    lagtime,  # X-axis coordinate
    msd,  # vertical axis
    xlim=None,  # Set the display range of the x-axis
    ylim=None,  # Set the display range of the y-axis
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10MSD.png",  # Output image filename
    show=False,  # Whether to pop up the image window
    ax=None,  # Optional subplot specification
)

Executing the code will generate an image similar to the following:

警告

If you execute the above script by SSH connection to a remote server and encounter QT-related error messages, it might be due to incompatibility between the program you're using (like MobaXterm, etc.) and the QT libraries. Either switch to a different program (such as VSCode or the system's built-in terminal), or add the following code starting from the second line of the Python script:

import matplotlib
matplotlib.use('agg')

Root Mean Square Deviation (RMSD) Analysis

See 10aimd_rmsd.py:

# coding:utf-8
from dspawpy.analysis.aimdtools import get_lagtime_rmsd, plot_rmsd

# If AIMD is not completed in a single run, you can assign multiple paths of h5 files in list form to the datafile parameter
lagtime, rmsd = get_lagtime_rmsd(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",
    timestep=None,  # Data file path  # Default reads the time step from the datafile
)
plot_rmsd(
    lagtime,  # Horizontal coordinate
    rmsd,  # vertical coordinate
    xlim=None,  # Set the display range of the x-axis
    ylim=None,  # Set the display range of the y-axis
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10RMSD.png",  # Output image filename
    show=False,  # Whether to pop up the image window
    ax=None,  # Optional subplot specification
)

Executing the code will generate an image similar to the following:

警告

If you execute the above script by connecting to a remote server via SSH, and QT-related error messages appear, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT libraries. Either change the program (such as VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

Analysis of Atomic Pair Distribution Functions or Radial Distribution Functions (RDFs)

See 10aimd_rdf.py :

# coding:utf-8
from dspawpy.analysis.aimdtools import get_rs_rdfs, plot_rdf

# If AIMD is not completed in one go, you can assign multiple h5 file paths to the datafile parameter in the form of a list.
rs, rdfs = get_rs_rdfs(
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.18/aimd.h5",  # Data file path
    ele1="H",  # Central element
    ele2="O",  # Target element
    rmin=0.0,  # Minimum radius
    rmax=10.0,  # Maximum radius
    ngrid=1000,  # Number of grid points
    sigma=0.1,  # sigma value
)
plot_rdf(
    rs,  # x-axis values
    rdfs,  # Vertical coordinate
    "H",  # Central element
    "O",  # Object element
    figname="dspawpy_proj/dspawpy_tests/outputs/us/10RDF.png",  # Image save path
    show=False,  # Whether to pop up the image window
    ax=None,  # Subplot can be specified
)

Executing the code will generate an image similar to the following:

The statistical calculations involved in this section are complex; please refer to the function API for more details.

警告

If you connect to a remote server via SSH and execute the above script, and you encounter QT-related error messages, it's possible that the program you're using (such as MobaXterm) is incompatible with the QT libraries. Either change the program (for example, VSCode or the system's built-in terminal command line), or add the following code, starting on the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.11 Ferroelectric Polarization Data Processing

For a quick start, take the series of scf.h5 files obtained from ferroelectric calculations on the \(HfO_{2}\) system as an example:

See 11Ferri.py:

# coding:utf-8
from dspawpy.plot import plot_polarization_figure

plot_polarization_figure(
    directory="dspawpy_proj/dspawpy_tests/inputs/2.20",  # Path for iron polarization calculation
    repetition=2,  # Number of times to repeat the data points when plotting
    figname="dspawpy_proj/dspawpy_tests/outputs/us/11pol.png",  # Output polarization figure filename
    show=False,  # Whether to display the polarization figure
)  # --> pol.png

Executing the code will generate the following combined figure:

The ferroelectric polarization values for the head and tail configurations can be found below:

from dspawpy.plot import plot_polarization_figure

python
plot_polarization_figure(directory='.', annotation=True, annotation_style=1)  # Displays the ferroelectric polarization values for the initial and final configurations.

The code will generate the following combined plot:

Alternatively, a second annotation style can be used:

from dspawpy.plot import plot_polarization_figure

plot_polarization_figure(directory='.', annotation=True, annotation_style=2)  # Displays the ferroelectric polarization values for the initial and final configurations.

The code will generate the following combined plot:

警告

If you encounter QT-related error messages when executing the above script via SSH connection to a remote server, it may be due to incompatibility between the program used (e.g., MobaXterm) and the QT library. Either change the program (e.g., VSCode or the system's built-in terminal command line), or add the following code starting from the second line of your Python script:

import matplotlib
matplotlib.use('agg')

8.12 Zero-Point Vibrational Energy Data Processing

Taking the frequency.txt file obtained from a quick start \(CO\) system frequency calculation as an example, the zero-point vibrational energy is calculated based on the following formula:

\[ZPE=\sum_{i=1}^{3 N} \frac{h \nu_i}{2}\]

where \(\nu_i\) are the normal mode frequencies, \(h\) is Planck's constant (\(6.626\times10^{-34} J\cdot s\)), and \(N\) is the number of atoms.

See 12getZPE.py:

# coding:utf-8
from dspawpy.io.utils import getZPE

# Import the frequency.txt file obtained from frequency calculation
getZPE(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
)

The code execution results will be saved to the ZPE.dat file, and the file content is as follows:

Data read from D:\quickstart\2.13\frequency.txt
Frequency (meV)
284.840038

--> Zero-point energy,  ZPE (eV): 0.142420019

8.13 TS Hot Calibration Energy

Contribution of the entropy change of the adsorbate to the energy

Calculation is based on the following formula:

\[\Delta S_{a d s}\left(0 \rightarrow T, P^0\right)=S_{v i b}=\sum_{i=1}^{3 N}\left[\frac{N_{\mathrm{A}} h \nu_i}{T\left(e^{h \nu_i / k_{\mathrm{B}} T}-1\right)}-R \ln \left(1-e^{-h \nu_i / k_{\mathrm{B}} T}\right)\right]\]

Here, \(\Delta S_{a d s}\) represents the entropy change of the adsorbate, calculated based on the harmonic approximation. \(S_{v i b}\) represents the vibrational entropy. \(\nu_i\) is the normal mode frequency, \(N_A\) is Avogadro's constant (\(6.022\times10^{23} mol^{-1}\)), \(h\) is Planck's constant (\(6.626\times10^{-34} J\cdot s\)), \(k_B\) is the Boltzmann constant (\(1.38\times10^{-23} J\cdot K^{-1}\)), \(R\) is the ideal gas constant (\(8.314 J\cdot mol^{-1}\cdot K^{-1}\)), and \(T\) is the system temperature in units of \(K\).

See 13getTSads.py for reference:

# coding:utf-8
from dspawpy.io.utils import getTSads

# Import the frequency.txt file calculated from frequency, temperature can be modified arbitrarily
getTSads(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
    T=298.15,
)

The execution result will be saved to the TS.dat file, with the following content:

Data read from D:\quickstart\2.13\frequency.txt
Frequency (THz)
68.873994

--> Entropy contribution,  T*S (eV): 4.7566990201851275e-06

Ideal gas entropy contribution to energy

Calculations are based on the following formula:

\[\begin{split}\begin{aligned} S(T, P) & =S\left(T, P^{\circ}\right)-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \\ & =S_{\text {trans }}+S_{\text {rot }}+S_{\text {elec }}+S_{\text {vib }}-k_{\mathrm{B}} \ln \frac{P}{P^{\circ}} \end{aligned}\end{split}\]

Where:

\[S_{\text {trans }}=k_{\mathrm{B}}\left\{\ln \left[\left(\frac{2 \pi M k_{\mathrm{B}} T}{h^2}\right)^{3 / 2} \frac{k_{\mathrm{B}} T}{P^{\circ}}\right]+\frac{5}{2}\right\}\]

\[\begin{split}S_{\mathrm{rot}}= \begin{cases}0 & \text {, monatomic } \\ k_{\mathrm{B}}\left[\ln \left(\frac{8 \pi^2 I k_{\mathrm{B}} T}{\sigma h^2}\right)+1\right] & \text {, linear } \\ k_{\mathrm{B}}\left\{\ln \left[\frac{\sqrt{\pi I_A I_{\mathrm{B}} I_{\mathrm{C}}}}{\sigma}\left(\frac{8 \pi^2 k_B T}{h^2}\right)^{3 / 2}\right]+\frac{3}{2}\right\} & \text {, nonlinear }\end{cases}\end{split}\]

\[S_{\text {elec }}=k_{\mathrm{B}} \ln [2 \times(\text { total spin })+1]\]

\[S_{\text {vib }}=k_{\mathrm{B}} \sum_i^{\text {vib DOF }}\left[\frac{\epsilon_i}{k_{\mathrm{B}} T\left(e^{\epsilon_i / k_{\mathrm{B}} T}-1\right)}-\ln \left(1-e^{-\epsilon_i / k_{\mathrm{B}} T}\right)\right]\]

where \(I_A\) to \(I_C\) are the three principal moments of inertia for a non-linear molecule, \(I\) is the degenerate moment of inertia for a linear molecule, and \(\sigma\) is the symmetry number of the molecule. Furthermore, monatomic refers to a monatomic molecule, linear refers to a linear molecule, and nonlinear refers to a non-linear molecule. total spin is the total spin quantum number. vib DOF represents vibrational degrees of freedom.

Refer to the 13getTSgas.py script for processing:

# coding:utf-8
from dspawpy.io.utils import getTSgas

# Read elements and coordinates from the calculation result file (json or h5)
TSgas = getTSgas(
    fretxt="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt",
    datafile="dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.h5",
    potentialenergy=-0.0,
    geometry="linear",
    symmetrynumber=1,
    spin=1,
    temperature=298.15,
    pressure=101325.0,
)
print("Entropy contribution, T*S (eV)：", TSgas)

# If only the frequency.txt file is available, the calculation can be completed by manually specifying the elements and coordinates
# TSgas = getTSgas(fretxt='dspawpy_proj/dspawpy_tests/inputs/2.13/frequency.txt', datafile=None, potentialenergy=-0.0, elements=['H', 'H'], geometry='linear', positions=[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], symmetrynumber=1, spin=1, temperature=298.15, pressure=101325.0)

8.14 Appendix

Quickly download all scripts by clicking UserScripts.zip
dspawpy Changelog

8 Auxiliary Tool User Guide

8.1 Installation and Updates

Update dspawpy

8.2 structure structure conversion

8.3 Volumetric Data Processing

volumetricData Visualization

Differential volumetric data visualization

Volumetric data plane average

8.4 band data processing

Conventional Band Treatment

The band is projected onto each element separately, with the size of the data points representing the element's contribution to the orbital.

Band projections onto different elements' different orbitals

Projecting band structure onto different atomic orbitals

Band unfolding processing

band-compare band structure comparison figure processing

8.5 DOS Data Processing

Total Density of States

Project Density of States onto different orbitals

Projecting the density of states onto different elements

Projecting the density of states onto different orbitals of different atoms

Projecting the density of states onto the split d-orbitals (t2g, eg) of different atoms

d-centered analysis

8.6 bandDos: Displaying Band Structure and Density of States Together

Display band structure and density of states in a single figure.

Display band structure and projected density of states on a single plot.

8.7 optical data processing

8.8 neb data processing

Generating intermediate configurations for input files

Plotting the energy barrier diagram

neb.iniFin = true/false

neb.iniFin = true

Processing Data for Transition State Calculations

Observing the NEB Chain

Calculate the inter-configuration distance

Continued calculation with neb

Energy and maximum atomic force variation trend during NEB calculation

8 LOOP 1:

8.9 Phonon Data Processing

Phonon band data processing

Phonon Density of States Data Processing

Phonon Thermodynamic Data Processing

8.10 aimd molecular dynamics simulation data processing

Convert the trajectory file format to .xyz or .dump.

Changes in energy, temperature, etc. curves during the dynamics process.

Mean Squared Displacement (MSD) Analysis

Root Mean Square Deviation (RMSD) Analysis

Analysis of Atomic Pair Distribution Functions or Radial Distribution Functions (RDFs)

8.11 Ferroelectric Polarization Data Processing

8.12 Zero-Point Vibrational Energy Data Processing

8.13 TS Hot Calibration Energy

Contribution of the entropy change of the adsorbate to the energy

Ideal gas entropy contribution to energy

8.14 Appendix