Single-point calculations

The simplest task that can be executed by ChemShell is an energy (and optionally gradient) evaluation on a single structure at a given level of theory. ChemShell contains interfaces to many external programs which can perform such calculations and presents a common interface to them, which makes switching from one to another straightforward.

Calling external codes from ChemShell

For most of ChemShell’s interfaces, the external code is executed using a system call in much the same way as you would run it from the command line.

Normally this requires the external code to be in your $PATH environment variable, although it is also usually possible to specify a full executable path in the input script. Occasionally other environment settings are also needed; for details see the relevant interface page in the manual.

In this tutorial, our examples use the NWChem, GAMESS-UK and ORCA electronic structure packages, but they can be easily modified to use other programs with ChemShell interfaces.

In the case of GAMESS-UK, you should put the location of your distribution’s rungamess script in your $PATH. Similarly for NWChem and ORCA the location of the nwchem and orca executables respectively should be in the path.

NB: for some external codes which can be compiled from source (such as GAMESS-UK and NWChem), it is possible to directly link the code into ChemShell as a library at compile time. This has advantages when running large-scale parallel jobs, as the two programs can share the same MPI environment, which means they run more efficiently and advanced techniques such as task-farming can be used. The installation guide has instructions for linking-in external codes to ChemShell.

A simple calculation on water

The example inputs hf_water.nwchem.py, hf_water.gamess-uk.py and hf_water.orca.py in samples/sp contain simple Hartree-Fock calculations on water using NWChem, GAMESS-UK and ORCA respectively.

Before running the calculation let’s look at the input files in detail.

The water geometry is contained in a standard XYZ file water.xyz:

3
Water
O       0.0000000      0.0000000      0.0000000
H      -0.7546067      0.5900326      0.0000000
H       0.7546067      0.5900326      0.0000000

This is often the simplest way to provide a geometry to ChemShell, although other ways will be looked at in later tutorials.

In the ChemShell input hf_water.nwchem.py, the geometry is read as follows:

# Read in water molecule geometry from disk
my_mol = Fragment(coords="water.xyz")

This command creates a Fragment object, which is the ChemShell data structure used to hold molecular information. The coords= option specifies the file to read the geometry from, which can take a number of formats including XYZ, and the resulting Fragment object is called my_mol.

The level of theory for the calculation can then be specified:

# Set up the QM calculation
my_theory = NWChem(frag=my_mol, method="hf", basis="3-21g")

Here, we create an NWChem object containing the options for the calculation. The frag= option specifies that we will use the water molecule already read in, while the Hartree-Fock level of theory is specified by method="hf". Additional options can then be specified for the HF calculation, such as specifying a basis set, here 3-21G.

Notice that the name of the object changes in the GAMESS-UK and ORCA inputs, but the settings remain unchanged. Wherever options are available across multiple software packages, ChemShell aims to present a consistent interface to allow easy switching from one package to another.

Finally, the calculation is run:

# Request a single-point calculation
my_task = SP(theory=my_theory)
my_task.run()

or equivalently:

my_task = SP(theory=my_theory).run()

Here we have specified that a single-point (SP) calculation should be run by creating an SP object, with the theory= option referencing the level of theory we have set already in my_theory.

Calling the run method then starts the calculation.

Providing ChemShell has been installed correctly, you can run the calculation using the command in a Unix shell session (e.g., bash) as follows (using NWChem as an example):

$ chemsh hf_water.nwchem.py

In the output, you will see a ChemShell banner giving information on the installation, followed by some information on the command being run. If the calculation is successful, the calculated energy will be reported under the heading Final SP energy.

It is also possible to see the NWChem input generated by ChemShell and output in the files _nwchem.inp and _nwchem.out. This is particularly useful if you encounter problems with the calculation and need to look at how it has progressed in more detail.

If you have access to more than one of the QM programs, try comparing the energies reported from each. You should find they agree to around 7 decimal places.

Density functional theory

The file dft_water.nwchem.py and the GAMESS-UK and ORCA equivalents illustrate how to carry out a density functional theory (DFT) calculation with ChemShell. The input is very similar to the Hartree-Fock calculation, except that method="dft" is selected and a functional should be specified with the functional= option, in this case BLYP:

# Set up the QM calculation
my_theory = NWChem(frag=my_mol, method="dft", functional="blyp", basis="3-21g")

The calculation is then run exactly as in the Hartree-Fock case:

# Request a single-point calculation and store the result
my_task = SP(theory=my_theory)
my_task.run()

ChemShell stores the resulting energy in a Result object, which can be accessed (in this case) as my_task.result.energy. This value can be stored in a variable for later use as follows:

energy1 = my_task.result.energy

The DFT example also shows how to carry out multiple calculations within a single input file. This flexibility is one of the key advantages of using ChemShell scripts over writing input files for the external codes directly. In the example, we simply create a second level of theory and task for the second calculation. Here we choose to run another BLYP calculation on the same molecule, but with a larger 6-31G basis set:

# Set up a second QM calculation with a different basis set
my_theory2 = NWChem(frag=my_mol, method="dft", functional="blyp", basis="6-31g")

# Request a new single-point calculation and store the result
my_task2 = SP(theory=my_theory2)
my_task2.run()
energy2 = my_task2.result.energy

Note

An alternative approach, if you do not need to keep the Theory options set up for the first calculation, is to modify the options of my_theory and run it again instead of creating a new my_theory2. The code for this, replacing the block above, would be:

# Change the basis set (only) and run calculation again
my_theory.basis = "6-31g"
my_task.run()
energy2 = my_task.result.energy

Note that this approach will overwrite the result in my_task.result.energy, but the original value will remain stored in the variable energy1 for future reference. Try experimenting with both ways of performing the second calculation - you should find the results are the same.

After the second calculation is run, we can print out a summary of the results:

# Report the results
print("\n Summary of energies")
print(" 3-21g energy: ", energy1)
print(" 6-31g energy: ", energy2)

Try running this example and note the difference in energies between the two calculations. Is the 6-31G energy higher or lower than the 3-21G energy?

If you have access to more than one of the QM programs, try comparing the reported energies from each again. DFT energies will not agree as closely as Hartree-Fock due to differences in the integration grids used in each code, but you should still see good agreement to around 4-5 decimal places.

If you compare codes in your own calculations and notice discrepancies, there can be several reasons. Handling of basis sets may differ (e.g. have you selected cartesian or spherical harmonic basis functions?), DFT functionals may not be defined consistently between codes (especially and notoriously B3LYP), or convergence criteria may be different. It is usually possible to adjust for such differences using ChemShell interface options if necessary.