*************************
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 <../../../manual/build/html/install.html>`_ 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 <https://en.wikipedia.org/wiki/Shell_(computing)>`__
session (e.g., `bash <https://en.wikipedia.org/wiki/Bash_(Unix_shell)>`__)
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.

