The calculations of the NMR shielding density were done in the DIRAC software. For details please visit DIRAC tutorial.
This directory contains:
After unpacking DIRAC_data.tar.gz, the directory DIRAC_data contains the data grouped by system under study: h2se_h2o, h2te_h2o and h2po_h2o. In each of these directories, there are subdirectories specific to the method used.
For instance the structure of h2se_h2o directory is as follows:
Analogous directories are created for two other systems, h2te_h2o and h2po_h2o.
In this work we were interested in the NMR shielding density of the heavy nuclei X (X = Se, Te, Po) and the hydrogen nuclei Hb, which participates in the hydrogen bond. For both nuclei we calculated the xx, yy and zz components of the NMR shielding tensor density separately, therefore inside each visual_* directory listed above, there are the directories corresponding to one of the components of the NMR shielding tensor density calculated for X or Hb nucleus:
inside visual_h2seh2o_super and visual_h2se_me there are:
inside visual_h2o_me there are:
inside visual_h2se_fde4 and visual_h2se_fde0 there are:
plot_vis_shield_0_xx_se, plot_vis_shield_0_yy_se, plot_vis_shield_0_zz_se
plot_vis_shield_v_xx_hb, plot_vis_shield_v_yy_hb, plot_vis_shield_v_zz_hb
plot_vis_shield_v_xx_se, plot_vis_shield_v_yy_se, plot_vis_shield_v_zz_se
plot_vis_shield_v_w11_xx_hb, plot_vis_shield_v_w11_yy_hb, plot_vis_shield_v_w11_zz_hb
plot_vis_shield_v_w11_xx_se, plot_vis_shield_v_w11_yy_se, plot_vis_shield_v_w11_zz_se
plot_vis_shield_v_w11_w12_xx_hb, plot_vis_shield_v_w11_w12_yy_hb, plot_vis_shield_v_w11_w12_zz_hb
where 0 denotes that no FDE-LAO terms were added to the property gradient in the linear response calculations, v - only the FDE-LAO embedding potential is present, v_w11 - the FDE-LAO embedding potential and uncoupled kernel and v+w11_w12 - all FDE-LAO terms are present, including the FDE-LAO coupling kernel contributions (please see the description in the manuscript).
inside visual_h2o_fde4 and visual_h2o_fde0 there are:
plot_vis_shield_nics_0_xx_se, plot_vis_shield_nics_0_yy_se, plot_vis_shield_nics_0_zz_se
plot_vis_shield_nics_v_xx_hb, plot_vis_shield_nics_v_yy_hb, plot_vis_shield_nics_v_zz_hb
plot_vis_shield_nics_v_xx_se, plot_vis_shield_nics_v_yy_se, plot_vis_shield_nics_v_zz_se
plot_vis_shield_nics_v_w11_xx_hb, plot_vis_shield_nics_v_w11_yy_hb, plot_vis_shield_nics_v_w11_zz_hb
plot_vis_shield_nics_v_w11_xx_se, plot_vis_shield_nics_v_w11_yy_se, plot_vis_shield_nics_v_w11_zz_se
plot_vis_shield_nics_v_w11_w12_xx_hb, plot_vis_shield_nics_v_w11_w12_yy_hb, plot_vis_shield_nics_v_w11_w12_zz_hb
where 0, v, v_w11 and v+w11_w12 abbreviations were described above.
In directories with the data for H2O molecule (visual_h2o_me, visual_h2o_fde4, visual_h2o_fde0), the abbreviation nics in the directories names refers to the fact that the current density induced by an external magnetic field in H2O molecule was contracted with the position vector pointing either to Se (plot_*_se directories) or to Hb (plot_*_hb directories) nuclei, as would be done in NICS calculations.
All these directories contain plot.3d.scalar files copied directly from the DIRAC calculations.
The data on these plot.3d.scalar files can be visualized as a scalar field on a 3D grid. It would show the density corresponding to a particular component of the NMR shielding tensor (xx, yy or zz).
To construct the isotropic NMR shielding density data for visualization we can use a small python script, here shown for the NMR shielding density of Se nucleus in H2Se molecule calculated with mechanical embedding:
import pandas as pd
data_xx = pd.read_csv('DIRAC_data/h2se_h2o/visual_h2se_me/plot_vis_shield_xx_se/plot.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_yy = pd.read_csv('DIRAC_data/h2se_h2o/visual_h2se_me/plot_vis_shield_yy_se/plot.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_zz = pd.read_csv('DIRAC_data/h2se_h2o/visual_h2se_me/plot_vis_shield_zz_se/plot.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_iso = 0.333333*(data_xx["s"]+data_yy["s"]+data_zz["s"])
Now, data_iso is a vector of the isotropic NMR shielding density in each grid point. This operation was possible as all densities are calculated on the same grid (created for the corresponding supermolecule) for a particular system under study.
Now, we may print this vector and compare it with the last column on a plot.vis_shield.me.se.iso.3d.scalar file in the DIRAC_data/h2se_h2o/modified directory
data_iso
This isotropic shielding density vector is stored on file in the directory modified in each system directory (h2se_h2o, h2te_h2o, h2po_h2o).
The names of files in modified directory are self-explanatory, for instance:
etc.
In the next step, we may be interested in the differential isotropic NMR shielding density.
For instance to show the difference between the NMR shielding density of Se in a supermolecule and the same density in H2Se and H2O subsystems calculated with FDE(4) method and FDE-LAO embedding potential included in the property gradient, we slightly modify the script above:
import pandas as pd
data_super = pd.read_csv('DIRAC_data/h2se_h2o/modified/plot.vis_shield.super.se.iso.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_h2se = pd.read_csv('DIRAC_data/h2se_h2o/modified/plot.vis_shield.fde4_v.se.iso.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_h2o = pd.read_csv('DIRAC_data/h2se_h2o/modified/plot.vis_shield.nics.fde4_v.se.iso.3d.scalar', sep=r"\s+", names=["x", "y", "z", "s"])
data_diff = data_super["s"]-data_h2se["s"]-data_h2o["s"]
Again, we may check the result and compare it with the last column on the file plot.vis_shield.super-fde4_v-nics.se.iso.3d.scalar in DIRAC_data/h2se_h2o/modified directory
data_diff
This data is saved in a text format (*.scalar) and in a cube format (*.cube) in modified directory.
The cube files can be visualized with many available programs used in chemistry. The plots presented in the manuscript were generated from the text files (*.scalar) using mayavi scripts available in the scripts directory (as an example run the plot.py file attached there; requirements: python > 3.0 and mayavi 4.5.0) and here.
All plots in png format are stored in plots directory. Each plot presents a 3D surface corresponding to the +- isovalue, which is indicated in the file name, for instance: plot_vis_shield_fde4_v_w11-fde4_v_se_iso_3d_scalar_0_01.png file presents the isosurfaces of +0.01 and -0.01 isovalues (0_01) obtained as a difference between the shielding density of Se (se) in the SeH2 molecule calculated with the fde(4) method and v+w11 FDE-LAO terms in the property gradient and the corresponding density obtained with the fde(4) method and v FDE-LAO term in the property gradient (fde4_v_w11-fde4_v).