Basis set based orbital generation with NWChem

In this tutorial, you will get a brief introduction to basis set based orbital generation with NWChem. First, the parameters for the multiwavelet representation are defined, and then a NWChem calculation is performed. The resulting atomic orbitals (AOs) and molecular orbitals (MOs) are read with the FrayedEnds NWChem converter and translated into multiwavelets which can be used for further calculations. Finally, the orbitals are plotted with py3Dmol.

Parameters

distance = 2.5
iteration_energies = []
iterations = 6
molecule_name = "beh2"
box_size = 50.0
wavelet_order = 7
madness_thresh = 0.0001
basisset = "cc-pvdz"

Run NWChem calculation

Create NWChem input and run NWChem calculation. If the FrayedEnds devcontainer or the singularity image is used, NWChem is already installed. Otherwise, NWChem has to be installed and the path has to be adjusted.

import subprocess as sp

nwchem_input = (
    """
title "molecule"
memory stack 1500 mb heap 100 mb global 1400 mb
charge 0  
geometry noautosym nocenter
  Be 0.0 0.0 0.0
  H 0.0 0.0 """
    + distance.__str__()
    + """
  H 0.0 0.0 """
    + (-distance).__str__()
    + """
end
basis  
  * library """
    + basisset
    + """
end
scf  
 maxiter 200
end   
task scf  
"""
)

with open("nwchem", "w") as f:
    f.write(nwchem_input)

programm = sp.call(
    "/opt/anaconda3/envs/frayedends/bin/nwchem nwchem",
    stdout=open("nwchem.out", "w"),
    stderr=open("nwchem_err.log", "w"),
    shell=True,
)

Convert NWChem AOs and MOs to MRA-Orbitals

Read the atomic orbitals (AOs) and molecular orbitals (MOs) from a NWChem calculation and translate them into multiwavelets.

import frayedends as fe

world = fe.MadWorld3D(L=box_size, k=wavelet_order, thresh=madness_thresh)

converter = fe.NWChem_Converter(world)
converter.read_nwchem_file("nwchem")
aos = converter.get_normalized_aos()
mos = converter.get_mos()

MADNESS runtime initialized with 9 threads in the pool and affinity OFF

Plotting of orbitals

First, a molecule object with the same geometry as in the NWChem calculation is created and then orbital(s) (here MO 5) are exported to a .cube file.

molecule = fe.MolecularGeometry(units="bohr")
molecule.add_atom(0.0, 0.0, 0.0, "Be")
molecule.add_atom(0.0, 0.0, distance, "H")
molecule.add_atom(0.0, 0.0, -distance, "H")

world.cube_plot("orbital", mos[5], molecule)

fe.cleanup(globals())

The .cube file can then be loaded and plotted interactively with py3Dmol, for example.

import py3Dmol

orbdata = open("orbital.cube", "r").read()
v = py3Dmol.view()
v.addVolumetricData(orbdata, "cube", {"isoval": -0.001, "color": "red", "opacity": 0.75})
v.addVolumetricData(orbdata, "cube", {"isoval": 0.001, "color": "blue", "opacity": 0.75})
v.addModel(orbdata, "cube")
v.setStyle({"sphere": {}})
v.zoomTo()
v.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.