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.