M2: Properties with MRChem

M2: Properties with MRChem#

How to compute magnetic properties with MRChem#

Below you will find links to all all the relevant pages of the MRChem documentation. Please refer to these if you need more details on how to set up the calculations.

If you need additional help, don’t hesitate to ask us or take a peak at these notebooks:

Introduction#

In this exercise you will use the MRChem code to compute NMR shielding constants for a water molecule, using the linear response implementation.

We provide some GTO NMR data that can serve as a comparison (computed with the ORCA code using aug-pcJ-1, aug-pcJ-2, aug-pcJ-3, and aug-pcJ-4 London orbitals).

The MWn notation is a shorthand for a relative precision of \(1\times 10^{-n}\), and is commonly used to refer to calculations of different precisions (analogous to the DZ, TZ, QZ… notation for GTO basis sets).

# 3D viewer for the water molecule
import py3Dmol

view = py3Dmol.view(width=200, height=200)
view.addModel(open('geometries/water.xyz').read())
view.setStyle({'stick': {}})
view.zoomTo()
view.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

GTO data:

[('dz', '0o', 319.317, 52.281),
 ('dz', '1h', 30.918, 18.326),
 ('dz', '2h', 30.918, 18.326),
 ('tz', '0o', 319.156, 52.351),
 ('tz', '1h', 30.766, 18.54),
 ('tz', '2h', 30.766, 18.54),
 ('qz', '0o', 319.183, 52.345),
 ('qz', '1h', 30.715, 18.62),
 ('qz', '2h', 30.715, 18.62),
 ('5z', '0o', 319.158, 52.345),
 ('5z', '1h', 30.71, 18.627),
 ('5z', '2h', 30.71, 18.627)]

The MRChem input file#

Below is an outline of an MRChem input for an NMR properties calculation at the MW3 precision level. Try to use the documentation pages to find the relevant keywords for setting up the calculation. If you need additional help, don’t hesitate to ask :)

world_prec = 1.0e-3
world_unit = angstrom

Molecule {
  $coords
  ...
  $end
}

WaveFunction {
}

SCF {
}

Properties {
}

Response {
}

Computational details#

The GTO results were computed with the PBE functional, and so you should also use PBE if the results are to comparable.
An optimized geometry of water is supplied in ./geometries/water.xyz.
Use a KAIN SCF and response accelerator history of 5.
Make sure the SCF and response max_iter options are large enough (20 should be plenty).
As starting guess for the SCF, sad_dz works well.

Note Input blocks are case sensitive, while keywords inside blocks are not

Instructions#

Generate the necessary input files. You can either do this in your favorite command line text editor, or you can use our Python scripts provided below from within the notebook. We give some tips on the notebook route below.
Submit the calculations as described above (terminal route) or below (notebook route).
You are encouraged to take a look at the output file (jobname.out) to look at the information printed here, and also to check out the JSON output file (jobname.json) to see the same output in a very scripting friendly format.
Extract the relevant information from either the jobname.out file or the jobname.json file (help on how to do this within the notebook given below).
Using the MW5 data as a reference, compute the basis set errors for the other methods. Try to visualize these either as a table or by plotting the errors (tips given below).

Performing the calculations#

Run calculations for at the MW3, MW4, and MW5 precision levels

import os
import pandas as pd
import matplotlib.pyplot as plt

from utils.functions import MRChemOutput, makeNMRInput, submit

Step 0: Set some global variables to be used for this exercise#

BASENAME   = ''             # Subsequent input and output files will be named <BASENAME_mwn>
NPROCS     = ''             # Number of OMP threads to be used in the calculations
PRECISIONS = {'mw3': 1e-3, 
              'mw4': 1e-4, 
              'mw5': 1e-5}  # precisions and labels

Step 1: Generate input file#

makeNMRInput will generate an input file in the current directory. Check out the convenience_scripts.ipynb Notebook for more help on how to use the function.

Fill in the missing data below to generate an input file in the current working directory.

Take a look at the generated input file before continuing.

Modify the precision and rerun to collect all your data.

prec_label = ''  # mw3, mw4, or mw5
jobname = f'{BASENAME}_{prec_label}'
xyzfile = ''

makeNMRInput(world_prec=PRECISIONS[prec_label], fname=jobname, xyzfile=xyzfile)

Step 2: Call MRChem and run the job#

submit will generate a new directory to which it moves the previously generated input file (which must exist). It then sets up the OMP_NUM_THREADS environment variable, and calls MRChem.

submit(nprocs=NPROCS, inputfile=jobname)

Step 3: Plot SCF convergence#

for label, prec in PRECISIONS.items():
    f = os.path.join(f'{jobname}_calc', f'{jobname}.json')
    calc = MRChemOutput(f)
    
    calc.plotSCFConvergence()
    calc.plotResponseConvergence()

Step 4: Collect data#

The code below will loop over your JSON output files and collect some relevant NMR data, and store it as a pandas.DataFrame.

mw_data = []
for label, prec in PRECISIONS.items():
    f = os.path.join(f'{jobname}_calc', f'{jobname}.json')
    calc = MRChemOutput(f)
    
    isos = calc.getNMRShieldingIsotropicAverage()
    anisos = calc.getNMRShieldingAnisotropy()
    
    for atom in isos.keys():
        mw_data.append((label, atom[4:].upper(), isos[atom], anisos[atom]))
        
gto_data = [('dz', '0O', 319.317, 52.281),
 ('dz', '1H', 30.918, 18.326),
 ('dz', '2H', 30.918, 18.326),
 ('tz', '0O', 319.156, 52.351),
 ('tz', '1H', 30.766, 18.54),
 ('tz', '2H', 30.766, 18.54),
 ('qz', '0O', 319.183, 52.345),
 ('qz', '1H', 30.715, 18.62),
 ('qz', '2H', 30.715, 18.62),
 ('5z', '0O', 319.158, 52.345),
 ('5z', '1H', 30.71, 18.627),
 ('5z', '2H', 30.71, 18.627)]

df = pd.DataFrame(mw_data + gto_data, columns=['Basis', 'Atom', 'Iso', 'Aniso'])
ref = df.loc[df.Basis == 'mw5']
df = df.merge(ref, how='inner', on='Atom', suffixes=('', '_ref'))
df.drop(df.loc[df.Basis == 'mw5'].index, inplace=True)

df['Error_iso'] = abs(df.Iso - df.Iso_ref) / df.Iso_ref * 100
df['Error_aniso'] = abs(df.Aniso - df.Aniso_ref) / df.Aniso_ref * 100

df = df.sort_values(by=['Error_iso']).reset_index(drop=True)

Step 5: Plot errors in isotropic averages#

fig, axes = plt.subplots(dpi=100, ncols=3, figsize=(10, 3))

for ax, atom in zip(axes.flat, df.Atom.unique()):
    sub_df = df.loc[df.Atom == atom]
    
    ax.bar(sub_df.Basis, sub_df.Error_iso, edgecolor='black', color='skyblue', lw=2)
    
    ax.set_ylabel('Relative Unsigned Error (%)')
    ax.set_title(f'Atom: {atom}')
    
plt.tight_layout(h_pad=1)

Step 6: Plot errors in anisotropies#

fig, axes = plt.subplots(dpi=100, ncols=3, figsize=(10, 3))

for ax, atom in zip(axes.flat, df.Atom.unique()):
    sub_df = df.loc[df.Atom == atom]
    
    ax.bar(sub_df.Basis, sub_df.Error_aniso, edgecolor='black', color='skyblue', lw=2)
    
    ax.set_yscale('log')
    ax.set_ylabel('Relative Unsigned Error (%)')
    ax.set_title(f'Atom: {atom}')
    

plt.tight_layout(h_pad=1)