Interfaces to external QM codes
Until now we have used the built-in MNDO module to carry out semiempirical QM calculations. However, ChemShell has interfaces to many different QM codes. This page describes how to perform calculations using other codes, with GAMESS-UK, NWChem and ORCA as examples.
Please see the QM interfaces manual pages for a full list of available interfaces.
Calling external codes from ChemShell
For most of ChemShell's QM code interfaces, the external code is executed as a binary using a system call in much the same way as a user would run it from the command line.
Normally this requires the external code to be in the user's $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 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 binaries respectively should be in the path.
Note 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 water_gamess-uk.chm, water_nwchem.chm and water_orca.chm in samples/qm contain simple water calculations using GAMESS-UK, NWChem and ORCA respectively.
As in the MNDO examples, the energy is retrieved from the matrix e and stored in a Tcl variable. The first calculation uses the default Hamiltonian (Hartree-Fock) and basis set (STO-3G) taken from the internal libraries of the relevant code. In the GAMESS-UK case this is run with:
energy coords=h2o.c theory=gamess energy=e
set e_sto3g [ get_matrix_element matrix=e indices= {0 0 } ]
For other codes, simply change theory= as appropriate. For example the NWChem case is:
energy coords=h2o.c theory=nwchem energy=e
set e_sto3g [ get_matrix_element matrix=e indices= {0 0 } ]
A second calculation then carried out to demonstrate how a different basis set can be chosen, e.g. 6-31G:
energy coords=h2o.c theory=gamess : {basis=6-31G } energy=e
set e_631g [ get_matrix_element matrix=e indices= {0 0 } ]
Finally a DFT calculation is performed with the BLYP functional using the same 6-31G basis set:
energy coords=h2o.c theory=gamess : { basis=6-31G hamiltonian=dft functional=blyp } energy=e
set e_blyp [ get_matrix_element matrix=e indices= {0 0 } ]
Note for technical reasons in the case of GAMESS-UK the functional should be specified with functional=, while for other codes it can be specified with hamiltonian=. Check before running to avoid errors.
At the end of the script a summary of the energies from each calculation is printed:
puts stdout "Energy (HF,STO-3G) was $e_sto3g" puts stdout "Energy (HF,6-31G) was $e_631g" puts stdout "Energy (BLYP,6-31G) was $e_blyp"
If you have access to more than one of the external QM codes you can compare results. For each of the three examples you should see good agreement between the codes.
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.
A subtlety in the water script is that the matrix object e is declared before use:
matrix e
This ensures that the matrix object is accessed from memory rather than disk each time it is used. For the tutorial example it does not make much difference to the time taken (try removing the line and re-running), but if you write a script that accesses an object many times (e.g. in a loop), you should always declare it in order to avoid inefficient disk I/O. See the manual for more information.
Test calculations for many other QM codes can be found in the examples/ directory of your ChemShell distribution. Apart from changing theory= you may also see changes to the theory options, especially for keywords such as basis which refer to the code's own basis set names. The interfaces aim to be consistent as far as possible in their handling of major options but you should be careful to check through the list whenever you change theory=.
More complex basis specifications
A common requirement is to specify a more complex basis set. There are two approaches possible with GAMESS-UK:
- Using the basisspec= argument, which reads in basis sets from ChemShell's internal basis set library. This ensures that the basis sets used are completely consistent when comparing different QM codes. It is also allows different basis sets to be applied to different atom types, as in the example basisspec.chm
- Providing basis sets and/or ECPs in separate files, using the basisfile= and ecpfile= options. This may be the simplest approach for users who are familiar with GAMESS-UK as you can use its standard format and all the options available within GAMESS-UK. See the example basisfile.chm, where the basis set is read in from basis.txt.
Similar options are available for other QM codes, depending on the details of the interface - please consult the relevant interface page.