ChemShell Tutorial

Search the manual:

Calculation of a Potential Energy Surface

Potential energy surfaces can be scanned in two ways:

1. Using general coordinate systems, then the DL-FIND optimiser is the most appropriate to use.
2. In z-matrix coordinates using the surface command.

1. Potential energy surface in general coordinates

This example is contained in the file doc/tutorial/samples/pe_surface/surface_dl-find.chm in the ChemShell tutorial.

A Tcl "for" loop has to be used to iterate over the desired values of the reaction coordinate. The starting value is specified in begin, the last value in end, and with each iteration the value rises by the value of step:

set begin   0.0
set end   360.0
set step   20.

# loop over the potential energy surface
for {set tor $begin } {$tor < $end} { set tor [ expr $tor + $step ] } {

  # do the optimisation for this step

}

DL-FIND can only constrain a particular quantity to its starting value. Thus, the geometry has to be altered in order to fulfil the desired value of the reaction coordinate (a torsion in this case) $tor. The alteration is done by using restraints as "theory" (i.e. energy functional):

dl-find coords=c active_atoms= [ list 2 1 6 9 ] \
    theory=restraint : [ list \
      restraints= [ list [ list torsion 2 1 6 9 $tor 3.0 ] ] ] \
    result=moved.c

Only the atoms actually involved in the reaction coordinate should be active.

Now, the constrained geometry optimisation is carried out:

dl-find coords=moved.c theory= dl_poly :  [ list $args_mm ] \
    constraints= { { torsion 2 1 6 9 } } result=c coordinates= dlc

These two steps are done in each iteration of the Tcl loop described above. The resulting potential energy surface can be obtained by writing it to standard output or to a file, see pe_surface/surface_dl-find.chm.

2. Potential energy surface in z-matrix coordinates

The surface command computes a series of energy values on a grid defined by variables of a ChemShell zmatrix object. The grid may be defined from a subset of all variables. Optionally, variables not defining the grid may be optimised at each grid-point.

surface is invoked in the following general fashion:

% surface zmatrix=zmatrix_object theory=computational_method_for_energy \
result=name_of_matrix_holding_results \
control=tcl_list_defining_grid_points

Arguments to the surface command are:

1. zmatrix=

   specifies the zmatrix object containing the (mixed Cartesian/zmatrix) definition of the molecule for which a energies are to be calculated.

2. theory=

   specifies the computational method for calculation of the energy

3. result=

   specifies the name of a matrix object to hold the resulting energies on the grid specified

4. control=

   specifies the grid-points to be calculated through a Tcl list of control records. Each record controls one internal coordinate (z-matrix) parameter, and consists of four fields:

   name #pts start step 
   

   1. name indicates the name of the parameter, as specified in the z-matrix.
   2. #pts is the number of steps to take in the parameter value.
   3. start is the value of the parameter at the first point.
   4. step is the increment to apply to the parameter.

Examples

% surface theory=mopac : keywords=am1 zmatrix=test.zmatrix \
 result=test1.grid control= {
{ roh1  5 1.9  0.05 }
{ ahoh  3 90.0 5.0  }
}
% print_matrix matrix=test1.grid
listing of matrix surface result 
                1           2           3   
1     :      -0.0868     -0.0900     -0.0915
2     :      -0.0850     -0.0880     -0.0894
3     :      -0.0823     -0.0852     -0.0865
4     :      -0.0790     -0.0817     -0.0829
5     :      -0.0750     -0.0776     -0.0786

In this example 5 instances of the variable roh1 and 3 instances of the variable ahoh are explored. They start off at the value 1.9 (bohr) and 90 (degrees), respectively and are incremented by 0.05 and 5.0, respectively. The script behind the command first varies the last specified variable.

The grid of points generated is a general matrix of rank n, where n is the number of variables defining the grid. Matrices of rank 2 may be printed as such, but for matrices of higher rank this is not possible. It is therefore important that the user be aware of the order in which the points are generated.

% surface theory=mopac : keywords=am1 zmatrix=test.zmatrix \
 result=test1.grid control= {
{ roh1  2 1.9  0.05 }
{ ahoh  3 90.0 5.0  }
{ roh2  2 2.0  0.1 }
}
% exec cat test1.grid
block=matrix records=0
block=matrix_title records=1
surface result
block=dense_real_matrix records=12 dimensions=2 3 2 
-8.20157811890000e-02
-8.00696789000000e-02
-8.48827991790000e-02
-8.27988878570000e-02
-8.61552060200000e-02
-8.39491157020000e-02
-7.46723553400000e-02
-7.26340925310000e-02
-7.72504094760000e-02
-7.50814388500000e-02
-7.82752262190000e-02
-7.59914315280000e-02
%