Write Simpson simulation files¶
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, 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 and Quantum Espresso 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.
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 (spin_rate = 40e3) 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 due 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]
# Get spin system
spinsys = write_spinsys(H_subset,
use_ms = True,
grad = {'H' : -1},
ref = {'H' : 31.7})
# Prepare Simpson simulation
simp_in = SimpSim(
out_name = 'ethanol_sim',
out_format = 'spe',
spinsys = spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 400e6,
start_op = 'Inx',
detect_op = 'Inp',
crystal_file = 'rep256',
gamma_angles = 6,
sw = 10e3,
verbose = '01',
lb = 10,
zerofill = 4096,
method = "direct",
tsw = 1e4,
pulse_sequence = no_pulse,
)
# Write input file
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.
Read VASP NMR calculations¶
Unlike CASTEP and Quantum Espresso 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. The read_vasp function handles these conversions and adds the required information. Here is one example how to use it, and how the results compare with a similar CASTEP calculation.
# Castep calculation
castep = read('../../examples/write/AlPO-14.magres')
p_idx = [atom.index for atom in castep if atom.symbol == 'P']
p_subset = castep[p_idx]
al_idx = [atom.index for atom in castep if atom.symbol == 'Al']
al_subset = castep[al_idx]
p_spinsys = write_spinsys(p_subset,
use_ms = True,
grad = {'P' : -1},
ref = {'P' : 283.839})
# Al is a quadrupolar so it is important to consider the quadrupolar interaction
# Simpson allows simulation of quadrupolar broadening up to 2nd order
al_spinsys = write_spinsys(al_subset,
use_ms = True,
grad = {'Al' : -1},
ref = {'Al' : 538.78},
q_order=2)
p_inp = SimpSim(
out_name = 'castep_sim_p',
out_format = 'spe',
spinsys = p_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inp',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 30e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 3e4,
pulse_sequence = no_pulse,
)
al_inp = SimpSim(
out_name = 'castep_sim_al',
out_format = 'spe',
spinsys = al_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inc',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 20e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 2e4,
pulse_sequence = no_pulse,
)
p_inp.save('../../examples/write/castep_sim_p.in')
al_inp.save('../../examples/write/castep_sim_al.in')
# VASP calculation
vasp = read_vasp('../../examples/write/AlPO-14.OUTCAR', format='vasp-out')
p_idx = [atom.index for atom in vasp if atom.symbol == 'P']
p_subset = vasp[p_idx]
al_idx = [atom.index for atom in vasp if atom.symbol == 'Al']
al_subset = vasp[al_idx]
p_spinsys = write_spinsys(p_subset,
use_ms = True,
grad = {'P' : -1},
ref = {'P' : 294.50})
al_spinsys = write_spinsys(al_subset,
use_ms = True,
grad = {'Al' : -1},
ref = {'Al' : 542.96},
q_order=2)
p_inp = SimpSim(
out_name = 'vasp_sim_p',
out_format = 'spe',
spinsys = p_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inp',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 30e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 3e4,
pulse_sequence = no_pulse,
)
al_inp = SimpSim(
out_name = 'vasp_sim_al',
out_format = 'spe',
spinsys = al_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inc',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 20e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 2e4,
pulse_sequence = no_pulse,
)
p_inp.save('../../examples/write/vasp_sim_p.in')
al_inp.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.
# Prepare Simpson for each Al atom of CASTEP
al_idx = [atom.index for atom in castep if atom.symbol == 'Al']
for i, atom in enumerate(al_idx):
al_subset = castep[[atom]]
al_spinsys = write_spinsys(al_subset,
use_ms = True,
grad = {'Al' : -1},
ref = {'Al' : 538.78},
q_order=2)
al_inp = SimpSim(
out_name = f'castep_sim_al_{i}',
out_format = 'spe',
spinsys = al_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inc',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 20e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 2e4,
pulse_sequence = no_pulse,
)
al_inp.save(f'../../examples/write/split_simulation_al/castep_sim_al_{i}.in')
# Prepare Simpson for each Al atom of VASP
al_idx = [atom.index for atom in vasp if atom.symbol == 'Al']
for i, atom in enumerate(al_idx):
al_subset = vasp[[atom]]
al_spinsys = write_spinsys(al_subset,
use_ms = True,
grad = {'Al' : -1},
ref = {'Al' : 542.96},
q_order=2)
al_inp = SimpSim(
out_name = f'vasp_sim_al_{i}',
out_format = 'spe',
spinsys = al_spinsys,
spin_rate = 40e3,
np = 2048,
proton_freq = 800e6,
start_op = 'Inx',
detect_op = 'Inc',
crystal_file = 'rep168',
gamma_angles = 6,
sw = 20e3,
verbose = '01',
lb = 200,
zerofill = 4096,
method = "direct",
tsw = 2e4,
pulse_sequence = no_pulse,
)
al_inp.save(f'../../examples/write/split_simulation_al/vasp_sim_al_{i}.in')
Now let's combined all the 27Al spectra and plot the results.
# Read CASTEP Al simulations
real = []
imag = []
for i in range(0, len(al_idx)):
castep_al_out = SimpReader(f'../../examples/write/split_simulation_al/castep_sim_al_{i}.spe', b0='800MHz', nucleus='27Al', format='spe')
real.append(castep_al_out.data['real'])
imag.append(castep_al_out.data['imag'])
# Read one example and update the data
castep_al_out = SimpReader('../../examples/write/split_simulation_al/castep_sim_al_0.spe', b0='800MHz', nucleus='27Al', format='spe')
castep_al_out.data['real'] = sum(real)
castep_al_out.data['imag'] = sum(imag)
# Do the same of VASP
real = []
imag = []
for i in range(0, len(al_idx)):
vasp_al_out = SimpReader(f'../../examples/write/split_simulation_al/vasp_sim_al_{i}.spe', b0='800MHz', nucleus='27Al', format='spe')
real.append(vasp_al_out.data['real'])
imag.append(vasp_al_out.data['imag'])
vasp_al_out = SimpReader('../../examples/write/split_simulation_al/vasp_sim_al_0.spe', b0='800MHz', nucleus='27Al', format='spe')
vasp_al_out.data['real'] = sum(real)
vasp_al_out.data['imag'] = sum(imag)
castep_p_out = SimpReader('../../examples/write/castep_sim_p.spe', b0='800MHz', nucleus='31P', format='spe')
vasp_p_out = SimpReader('../../examples/write/vasp_sim_p.spe', b0='800MHz', nucleus='31P', format='spe')
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].plot(castep_p_out.data['ppm'], castep_p_out.data['real'])
ax[0].plot(vasp_p_out.data['ppm'], vasp_p_out.data['real'])
ax[0].set_xlabel('$^{31}$P (ppm)')
ax[0].set_xlim(-10, -45)
ax[0].legend(['CASTEP', 'VASP'])
ax[1].plot(castep_al_out.data['ppm'], castep_al_out.data['real'])
ax[1].plot(vasp_al_out.data['ppm'], vasp_al_out.data['real'])
ax[1].set_xlabel('$^{27}$Al (ppm)')
ax[1].set_xlim(45, 0)
ax[1].legend(['CASTEP', 'VASP'])
plt.show()
