From DFT Calculations to SIMPSON Spectra¶
A Simpson input file is typically divided into four main sections:
- Spin System: This section identifies the nuclei to be simulated and includes their NMR parameters, such as isotropic shift, anisotropy, asymmetry, and Euler angles.
- Parameters: This defines the settings for the "virtual spectrometer," such as the magnetic field strength (
proton_freq), spinning rate (spin_rate), number of points (np), pulse lengths, pulse powers, etc. - Pulse Sequence: This part specifies the pulse program used in the simulation.
- Processing: Defines the processing parameters, such as line broadening (
lb) or Fourier transformation of the FID signal.
Simpyson simplifies the preparation of these input files by interfacing data from DFT calculations with Simpson, using ASE. The spin system for magres files from CASTEP NMR calculations can be easily generated with Soprano, while Simpyson can handle the rest of the input file.
For a more detailed description of each parameter, please refer to the following Simpson papers Bak et al, Tošner et al., and Juhl et al.
from __future__ import annotations
import shutil
import matplotlib.pyplot as plt
from ase.io import read
from soprano.calculate.nmr.simpson import write_spinsys
from simpyson.calculator import SimpCalc
from simpyson.converter import read_vasp
from simpyson.io import read_simp
from simpyson.utils import add_spectra
1. Reading DFT Outputs¶
In the following example, we start by reading the structure from an ethanol CASTEP calculation, selecting the H atoms, referencing the chemical shifts (grad = {'H' : -1} and ref = {'H' : 31.7}), and construct the input file for a spectrum acquired at 400 MHz (proton_freq = 400e6) with a spinning rate of 40 kHz (spin_rate = 40e3). To keep the simulation simple, we begin with the magnetization in the xy-plane (start_op = 'Inx'), i.e. the detection plane, thus removing the need to do pulses. For this, we can use the no_pulse sequence.
# Read magres file
ethanol = read('../../examples/write/ethanol.magres')
# Select only the H atoms
h_idx = [atom.index for atom in ethanol if atom.symbol == 'H']
H_subset = ethanol[h_idx]
# Build spin system from CASTEP NMR output via Soprano
spinsys = write_spinsys(H_subset,
use_ms=True,
grad={'H': -1},
ref={'H': 31.7})
simp_in = SimpCalc(
spinsys=spinsys,
out_name='ethanol_sim',
out_format='spe',
spin_rate=40e3,
np=2048,
proton_frequency=400e6,
start_operator='Inx',
detect_operator='Inp',
crystal_file='rep256',
gamma_angles=6,
sw=10e3,
verbose=0,
lb=10,
zerofill=4096,
method='direct',
tsw=1e4,
pulse_sequence='no_pulse',
)
simp_in.save('../../examples/write/ethanol_sim.in')
The results can be easily plotted. The difference between this spectrum and the one from the previous tutorial is due to the absence of J-coupling in this simulation.
2. Read VASP NMR calculations¶
Unlike CASTEP magres files, VASP's OUTCAR file contains the magnetic shielding and electric field gradient tensors in units that are not in Hz, as required by Soprano and Simpson. Additionally, ASE's read function doesn't automatically add the ms and efg arrays to the Atoms object, which are necessary for Soprano to prepare the spin system.
Simpyson can parse VASP OUTCAR files to extract NMR parameters like chemical shielding tensors and Electric Field Gradients (EFG). Use read_vasp for this, that creates a ASE Atoms object containing the magnetic shielding and efg tensor under atoms.arrays['ms'] and atoms.arrays['efg'], similar to Magres files.
from simpyson.converter import read_vasp
# Path to example OUTCAR
vasp_file = '../../examples/write/AlPO-14.OUTCAR'
atoms = read_vasp(vasp_file, format='vasp-out')
# Print first four ms and efg tensors
print('Example of magnetic shielding tensors')
print(f"{atoms.arrays['ms'][:4]}")
print('\n')
print('Example of electric field gradient tensors')
print(f"{atoms.arrays['efg'][:4]}")
2.1 Compare VASP and CASTEP calculations¶
# CASTEP
castep = read('../../examples/write/AlPO-14.magres')
p_idx = [atom.index for atom in castep if atom.symbol == 'P']
al_idx = [atom.index for atom in castep if atom.symbol == 'Al']
p_spinsys = write_spinsys(castep[p_idx],
use_ms=True, grad={'P': -1}, ref={'P': 283.839})
# Al is quadrupolar: include second-order quadrupolar broadening (q_order=2)
al_spinsys = write_spinsys(castep[al_idx],
use_ms=True, grad={'Al': -1}, ref={'Al': 538.78},
q_order=2)
SimpCalc(
spinsys=p_spinsys, out_name='castep_sim_p', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inp',
crystal_file='rep168', gamma_angles=6, sw=30e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=3e4, pulse_sequence='no_pulse',
).save('../../examples/write/castep_sim_p.in')
SimpCalc(
spinsys=al_spinsys, out_name='castep_sim_al', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inc',
crystal_file='rep168', gamma_angles=6, sw=20e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=2e4, pulse_sequence='no_pulse',
).save('../../examples/write/castep_sim_al.in')
# VASP
vasp = read_vasp('../../examples/write/AlPO-14.OUTCAR', format='vasp-out')
p_idx = [atom.index for atom in vasp if atom.symbol == 'P']
al_idx = [atom.index for atom in vasp if atom.symbol == 'Al']
p_spinsys = write_spinsys(vasp[p_idx],
use_ms=True, grad={'P': -1}, ref={'P': 294.50})
al_spinsys = write_spinsys(vasp[al_idx],
use_ms=True, grad={'Al': -1}, ref={'Al': 542.96},
q_order=2)
SimpCalc(
spinsys=p_spinsys, out_name='vasp_sim_p', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inp',
crystal_file='rep168', gamma_angles=6, sw=30e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=3e4, pulse_sequence='no_pulse',
).save('../../examples/write/vasp_sim_p.in')
SimpCalc(
spinsys=al_spinsys, out_name='vasp_sim_al', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inc',
crystal_file='rep168', gamma_angles=6, sw=20e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=2e4, pulse_sequence='no_pulse',
).save('../../examples/write/vasp_sim_al.in')
Unless you have a laptop/node with crazy amounts of memory, simulating the full Al spin system, containing 8 Al atoms each with second-order quadrupolar broadening q_order = 2, can be quite challenging. Since we're not considering coupling between the Al atoms, a more efficient approach is to prepare an input file for each individual Al atom and then merge all the spectra into a combined spectrum. This reduces the computational load and makes the simulation much more manageable.
al_idx = [atom.index for atom in castep if atom.symbol == 'Al']
# CASTEP — one input file per Al site
for i, atom in enumerate(al_idx):
al_spinsys = write_spinsys(castep[[atom]],
use_ms=True, grad={'Al': -1}, ref={'Al': 538.78},
q_order=2)
SimpCalc(
spinsys=al_spinsys, out_name=f'castep_sim_al_{i}', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inc',
crystal_file='rep168', gamma_angles=6, sw=20e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=2e4, pulse_sequence='no_pulse',
).save(f'../../examples/write/split_simulation_al/castep_sim_al_{i}.in')
al_idx = [atom.index for atom in vasp if atom.symbol == 'Al']
# VASP — one input file per Al site
for i, atom in enumerate(al_idx):
al_spinsys = write_spinsys(vasp[[atom]],
use_ms=True, grad={'Al': -1}, ref={'Al': 542.96},
q_order=2)
SimpCalc(
spinsys=al_spinsys, out_name=f'vasp_sim_al_{i}', out_format='spe',
spin_rate=40e3, np=2048, proton_frequency=800e6,
start_operator='Inx', detect_operator='Inc',
crystal_file='rep168', gamma_angles=6, sw=20e3,
verbose=0, lb=200, zerofill=4096, method='direct',
tsw=2e4, pulse_sequence='no_pulse',
).save(f'../../examples/write/split_simulation_al/vasp_sim_al_{i}.in')
Now let's combined all the 27Al spectra and plot the results.
al_idx = [atom.index for atom in vasp if atom.symbol == 'Al']
# Read and combine all per-site 27Al spectra
castep_al_spectra = [
read_simp(f'../../examples/write/split_simulation_al/castep_sim_al_{i}.spe',
format='spe', b0='800MHz', nucleus='27Al')
for i in range(len(al_idx))
]
castep_al_out = add_spectra(castep_al_spectra)
vasp_al_spectra = [
read_simp(f'../../examples/write/split_simulation_al/vasp_sim_al_{i}.spe',
format='spe', b0='800MHz', nucleus='27Al')
for i in range(len(al_idx))
]
vasp_al_out = add_spectra(vasp_al_spectra)
# Read single-site 31P spectra
castep_p_out = read_simp('../../examples/write/castep_sim_p.spe', format='spe', b0='800MHz', nucleus='31P')
vasp_p_out = read_simp('../../examples/write/vasp_sim_p.spe', format='spe', b0='800MHz', nucleus='31P')
# Plot
fig, ax = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)
ax[0].plot(castep_p_out.ppm['ppm'], castep_p_out.ppm['real'], label='CASTEP')
ax[0].plot(vasp_p_out.ppm['ppm'], vasp_p_out.ppm['real'], label='VASP')
ax[0].set_xlabel('$^{31}$P (ppm)')
ax[0].set_xlim(-10, -45)
ax[0].legend()
ax[1].plot(castep_al_out.ppm['ppm'], castep_al_out.ppm['real'], label='CASTEP')
ax[1].plot(vasp_al_out.ppm['ppm'], vasp_al_out.ppm['real'], label='VASP')
ax[1].set_xlabel('$^{27}$Al (ppm)')
ax[1].set_xlim(45, 0)
ax[1].legend()
plt.show()
