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
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
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
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
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
# 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
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.
frag= option specifies that we will use the water molecule already
read in, while the Hartree-Fock level of theory is specified by
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()
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
run method then starts the calculation.
Providing ChemShell has been installed correctly, you can run the calculation as follows (using NWChem as an example):
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.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¶
dft_water.nwchem.py and the GAMESS-UK and ORCA equivalents
illustrate how to carry out a density functional theory calculation with
ChemShell. The input is very similar to the Hartree-Fock calculation, except
method="dft" is selected and a functional should be specified with
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.
This 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 this case, 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")
After the second calculation is run, we 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.