Single-point energy evaluations
Calculating an energy
For practical reasons, each computational chemistry code interfaced into ChemShell is accessed by a number of Tcl procedures, for example one to initialise the parameters, one to compute an energy, one to compute the gradient based on an energy, and a termination function to be called when no more evaluations are needed. These are the procedures that generate the input files for the interfaced codes and read the results from their output.
This complexity is hidden by a number of "driver" functions, which call the code interface procedures as necessary to accomplish the required task, such as an energy evaluation, geometry optimisation or an MD run.
The file mndo.chm in samples/energy contains an example single-point energy evaluation using the simplest driver, energy.
First, a water structure is stored in a ChemShell fragment. This is the geometry that will be used in the energy calculation.
c_create coords=h2o.pun {
coordinates angstrom
o 0.0 0.0 0.0
h 0.0 0.0 1.0
h 0.0 0.8 -0.5
}
Once the geometry is created we can print out a list of bond lengths and angles with the bang command.
bang coords=h2o.pun
(Note that bang is for information only and is not required for the energy calculation).
The energy is then calculated using the internal MNDO semiempirical quantum mechanics code in ChemShell, with an AM1 hamiltonian. The choice of QM code is specified by the theory= argument, in this case theory=mndo. We also need to specify the structure we want to perform the calculation on, which is the fragment h2o.pun we have just created. The full command is as follows:
energy coords=h2o.pun theory=mndo : hamiltonian=am1
The colon symbol (":") is used for nesting arguments in ChemShell commands. In the case above the colon indicates that the hamiltonian keyword is an argument to mndo, not to energy.
Using the result
In the mndo.chm example the energy driver simply prints out the energy in the ChemShell output. If you need to store the energy to use in a subsequent calculation, you can pass the name of a ChemShell matrix object which you wish to hold the result. An example of this is in the mndo2.chm script.
A ChemShell "matrix" can be used to hold matrices but is also used for storing general numerical data. In this case, we use a matrix object called e to hold a single value, the energy (i.e. it is a 1x1 "matrix").
energy energy=e coords=h2o.pun theory=mndo : hamiltonian=am1
Note that unlike ordinary Tcl variables, the information in the matrix object e is not lost when ChemShell closes. It is stored in the file e in ChemShell format, and can be read back in by other scripts.
The energy can be retrieved from the matrix object using the command get_matrix_element. Note that matrix indices count from zero, so we want the "0 0" matrix element. In the example, the energy is stored in the Tcl variable e_mndo, which can then be used like any other variable (see the Tcl basics section).
set e_mndo [ get_matrix_element matrix=e indices= {0 0} ]
In the mndo2.chm example, we compare the retrieved result to a reference energy to check that the calculation was successful.
set e_ref -0.08248797 puts stdout "Calculated energy is $e_mndo" puts stdout "Reference energy is $e_ref" set e_diff [ expr $e_mndo - $e_ref ] puts stdout "Difference is [ format "%15.8f" $e_diff ]"
In addition to the ChemShell output, you can also usually check the QM package input file generated by ChemShell and the output from the QM program if you wish. In the case of MNDO, the input is mndo.in and the output is mndo.out. Normally you only need to check these files if there is a problem with the calculation.
Please note that the internal MNDO module supplied with ChemShell is a cut-down version intended for testing and training. To obtain the fully-featured package please contact the Thiel group.