M1: Energies with MRChem

M1: Energies with MRChem#

How to compute energies 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 perform simple energy optimizations of a water molecule, using different MultiWavelet precision levels. We provide some GTO energies that can serve as a comparison (computed with the ORCA code using aug-pc-1, aug-pc-2, aug-pc-3, aug-pc-4 basis sets).

The MWn notation is a shorthand for a MW precision of \(1\times10^{-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 = {
    'aug-dz': -76.337833062097,
    'aug-tz': -76.384930572433,
    'aug-qz': -76.388407845518,
    'aug-5z': -76.388593505849
}

The MRChem input file#

Below is an outline of an MRChem input for a simple energy optimization at the MW3 precision level. Try to use the documentation pages to find the relevant keywords for setting up the calculation.

world_prec = 1.0e-3
world_unit = angstrom

Molecule {
  $coords
  ...
  $end
}

WaveFunction {
}

SCF {
}

Computational details#

Note: Most of these settings are default, but it is good practice to be explicit.

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 accelerator history of 5.
Make sure the SCF max_iter option is large enough (20 should be plenty).
As starting guess, sad_dz works well.

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

Instructions#

Generate the necessary input files necessary. You can either do this in your favorite command line text editor (vim, emacs, nano, etc), or use our supplied Python scripts from within the notebook. Some tips are given below to help you get you started if you choose the notebook route.
Submit the calculations as described above (terminal based) or below (notebook based).
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, and store the results in a python variable.
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.

Performing the calculations#

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

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

from utils.functions import MRChemOutput, makeEnergyInput, submit
from grid_plotter import grid_plotter
from vampyr import vampyr3d as vp

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#

makeEnergyInput 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.

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

makeEnergyInput(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 the SCF convergence#

You can use our convenience class MRChemOutput to see how the SCFs converged. The code assumes that your output files are named according to <jobname>_<mwn>.json (n=3,4,5), and that these are located in directories named according to <jobname>_<mwn>_calc. If you used the script providede above, then this will automatically be the case.

Just modify the base variable to be the same as what you named your input files.

In addition, the loop below stores the final SCF energies in a list called mw_data that will be used in the subsequent cell for data analysis and visualization.

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

Step 4: Data analysis#

The code below will combine the MW and GTO data into a pandas.DataFrame, compute the basis set errors using your MW5 results as a reference, and visualize the errors in a matplotlib bar plot.

gto_data = [
    ('dz', -76.337833062097),
    ('tz', -76.384930572433),
    ('qz', -76.388407845518),
    ('5z', -76.388593505849)
]

df = pd.DataFrame(gto_data+mw_data, columns=['Basis', 'Energy'])
ref = df.loc[df.Basis == 'mw5']['Energy'].values[0]
df['Error'] = df['Energy'] - ref
df.drop(df.loc[df.Basis == 'mw5'].index, inplace=True)

# Plot
fig, ax = plt.subplots(dpi=100)
ax.set_yscale('log')
ax.set_ylabel('Basis set error (a.u.)')

ax.bar(df.Basis.str.upper(), df.Error, edgecolor='black', color='salmon', lw=3)

Step 5: Visualize multiresolution grid for optimized orbital#

As a last step you can use the grid_plotter tool to visualize the optimized multiresolution grid of one of the orbital files. Just update the basename variable, and the code should find your orbital file (assuming you used the supplied Python functions for submitting your calculations).

def get_mra(jsonfile):
    """Return a vp.MultiResolutionAnalysis object with 
       parameters read from the JSON output file."""
    with open(jsonfile) as f:
        j = json.loads(f.read())
    
    boxes = j['input']['mra']['boxes']
    corner = j['input']['mra']['corner']
    order = j['input']['mra']['basis_order']
    scale = j['input']['mra']['min_scale']
    
    return vp.MultiResolutionAnalysis(box=vp.BoundingBox(nboxes=boxes, corner=corner, scale=scale), order=order)

# Use the supplied grid_plotter to visualize the grid
calcdir = jobname + '_calc'
orb = os.path.join(calcdir, 'orbitals', f'phi_p_scf_idx_0_re')

mra = get_mra(os.path.join(calcdir, jobname+'.json'))
tree = vp.FunctionTree(mra)
tree.loadTree(orb)
fig, ax = grid_plotter(tree=tree)

# fig.savefig(f'MultiResGrid_{jobname}.png')