This program manages several different static calculations starting from a single structure of a system, similar to the typical functionality of a quantum chemistry package. Currently, the following functionalities are implemented:

• Energy/gradient calculation (single point)
• Frequency calculation (Hessian + Normal modes in molden format
• Geometry optimization to a local minimum
• Transition state (TS) optimization
• Reaction path optimization (IRC), starting from a TS

List of examples:

• 1: Egrad - Norbornene (energy and gradient calculation): explore/egrad_norbornene
• 2: Egrad with Orca (external call to Orca QM program): explore/egrad_orca
• 3: Optfreq - Norbornene (geometry optimization and frequency calculation): explore/optfreq_norbornene
• 4: TS + IRC - CH3+H2O (TS and IRC optimization): explore/ts_irc_ch4oh
• 5: Egrad - PES topoplogy (only determination of internal coordinates): explore/egrad_pes_topol

Folder: evbopt/egrad_norbornene

Energy and gradient calculations with explore.x can be performed for an arbitrary number of structures at once. They must be given in a trajectory file (xyz format) and need to have the same atom ordering as in the used potential energy functions. A simple example is the reaction path of the norbornene cycloaddtion/reversion (EVB optimization located here).

The dE-EVB-QMDFF has already been parametrized, it is defined by the following files:

• min1.qmdff: QMDFF of the reactants
• min2.qmdff: QMDFF of the products
• evb_pars.dat: Optimized dE-EVB coupling parameters

Besides those files, the keyfile scan.key is the most important one, defining what is to do:

pes de_evb

qmdff {
   ffnames min1.qmdff  min2.qmdff
   eshift    -272.22659664     -272.167983678
}

de_evb{
   coupling sp
}

xyzstart path.xyz
job egrad

Besides the PES specific keywords (coupling sp defining the exact coupling function for the dE coupling term, which parameters are listed in evb_pars.dat), only two keywords are needed for the calculation: xyzstart path.xyz, where the trajectory file containing the structures is defined and job egrad calling the calculation of energies and gradients.

After executing

explore.x scan.key

the energies and gradients are calculated in a fraction of a second:

...
EGRAD: Energies and gradients for a number of structures
 in the XYZSTART file will be calculated.
 . ...- -... --.- -- -.. ..-. ..-.
Calculating energy and gradient ...
Energies written to 'energies.dat', gradients written to 'gradients.dat'.
Energies of QMDFFs written to 'single_qmdff.dat'.
.. .----. -- -.. --- -. . 
Calculation successfully finished!
Output was written to explore.log
The calculation needed a time of     0.008 seconds.

As stated, energies and gradients can be found in energies.dat and gradients.dat (in Hartrees or Hartree/bohr). In this case, two QMDFFs are coupled to a EVB-QMDFF surface, therefore, the energies of the single QMDFFs are calculated as well and stored into single_qmdff.dat.

In this example, it is shown how the orca QM program can be used to calculate energies and gradients when they are needed. This procedure can be used as well for MD trajectories with dynamic.x or rate constant calculations with calc_rate.x. Keep in mind, however, that QM calculations of thousands of time steps can become very expensive! The example here is quite cheap and should be completed within a few minutes.

Main input is again the egrad.key file:

pes orca

xyzstart path.xyz
job egrad

orca {
   header_file header_inp.txt
   charge 0
   multi 1
   symlink orca.sh
}

The pes orca keyword activates the call to the Orca program. In the orca section, four keywords are needed to specify the external call:

• header_file: gives the name of the file in which the general keyword lines for the orca.inp file above the charge+multiplicity section are listed.
• charge: The total charge of the system
• multi: The spin multiplicity of the system
• symlink: The path to the orca executable (either absolute path or a command, if stored in the $PATH variables. Enter the Orca path in your system here!

In this case, the header_inp.txt file contains the following lines:

! PAL4
! PBE D3 cc-pVDZ engrad

Four processors are used, the PBE functional with Grimme D4 empirical dispersion correction and the cc-pVDZ basis is chosen. NOTE: It is important to add the engrad command, since the egrad.engrad file is needed to read in energy and gradient, which is only written by orca if this command is present!

After executing the explore.x program, the calculation will be finished after some minutes:

...
Used PES: direct call to orca QM package
Each structure will be submitted to an external orca
 calculation for each bead, separately!
----------------------------------------------------------

EGRAD: Energies and gradients for a number of structures
 in the XYZSTART file will be calculated.
 . ...- -... --.- -- -.. ..-. ..-.
Calculating energy and gradient ...
Energies written to 'energies.dat', gradients written to 'gradients.dat'.
.. .----. -- -.. --- -. . 
Calculation successfully finished!
Output was written to explore.log
The calculation needed a time of   121.988 seconds.

In the background, a new egrad.inp file is written for each structure, orca is called and after completing, energy and gradient are read from the egrad.engrad file. Energies and gradients can be read from the usual output files (see above).

Folder: evbopt/optfreq_norbornene

A simple example for a geometry optimization with a subsequent frequency calculation can be performed on the EVB-QMDFF surface of norbornene (EVB optimization located here). Both QMDFF files (min1.qmdff and min2.qmdff) as well as the EVB coupling parameters (evb_pars.dat) are needed for the setup of the potential energy surface, the initial structure (start.xyz) as well as the keyfile (optfreq.key) are needed for the geometry optimization calculation. Within optfreq.key, the following commands are listed:

pes de_evb

qmdff {
   ffnames  min1.qmdff  min2.qmdff
   eshift    -272.2265966470   -272.1679836787
   shift_manual
}

de_evb{
   coupling sp
}

xyzstart start.xyz
job optfreq

Besides the definition of the PES in the first part, only two keywords are needed for the optimization and frequency calculation. xyzstart start.xyz defines the xyz coordinate file of the initial structure, job optfreq calls the combined geometry optimization and frequency calculation.

The calculation is invoked by typing:

explore optfreq.key

All is done in a fraction of a second. General information is printed to screen:

OPTFREQ: A local minimum optimization followed by a 
    frequency calculation will be done.
Initializing geometry optimization...
Finished!
Succession of optimized structures written to geoopt.xyz
Final optimized structure written to opt_final.xyz
Calculate gradient and energy of optimized structure...
Calculate frequencies of the structure...
File explore.molden with normal mode spectrum written.
Finished!
.. .----. -- -.. --- -. . 
Calculation successfully finished!
Output was written to explore.log
The calculation needed a time of     0.071 seconds.

As can be seen, two output files are written. In explore.log, the energy, gradient, hessian and normal mode frequencies of the optimized structure are listed. Further, the progress of the optimization is printed:

iteration :    1  E= -272.14700411  Gnorm= 0.39056  displacement= 1.00000  alpha=  0.00  angle=  0.00
iteration :   10  E= -272.16230296  Gnorm= 0.04072  displacement= 0.13378  alpha=  0.80  angle=  0.49
iteration :   20  E= -272.16981593  Gnorm= 0.03834  displacement= 0.29551  alpha=  1.00  angle=  0.29
iteration :   30  E= -272.17049378  Gnorm= 0.00505  displacement= 0.01179  alpha=  1.00  angle=  0.34
iteration :   40  E= -272.17053522  Gnorm= 0.00055  displacement= 0.00348  alpha=  0.70  angle=  0.25
iteration :   50  E= -272.17053525  Gnorm= 0.00000  displacement= 0.00005  alpha=  1.00  angle=  0.35
The geometry optimization seems to be stagnating!
Since no further progress can be made, the optimization
will be ended here...
The geometry optimization has converged!
iteration :   56  E= -272.17053525  Gnorm= 0.00000  displacement= 0.00000  alpha=  0.60

In this case, the stagnating behavior of the optimization is no bad sign, the gradient norm (Gnorm) is already below the threshold, such that convergence is signaled. That indeed a minimum is reached can also be seen by looking at the normal mode frequencies at the bottom of explore.log, the largest frequency is roughly -4 cm-1. The progress of the optimization can be inspected by opening the file geoopt.xyz with a molecule viewer, the final optimized structure is also written to opt_final.xyz. Normal mode freqencies and the IR spectrum can be inspected with molden, by opening the file explore.molden. Since no dipole moments are calculated, no IR intensities can be calculated likewise; instead, equalized intensities are written.

It should be noted that several additional keywords can be altered to influence the progress of the optimization procedure. Type

explore.x -h

to list all of them. Usually, however, the standard settings should be sufficient.

Folder: evbopt/ts_irc_ch4oh

In this example, a TS optimization and a IRC generation will be done for the CH3+H2O system. This is described by the included analytical PES by Espiosa-Garcia (link. First, the TS optimization will be done. Its keyfile is ts.key:

pes ana_ch4oh
xyzstart ts_start.xyz

job opt_ts
opt {
  read_coord
}

The initial structure needs to be already a reasonable guess for the TS structure. This is ts_start.xyz in this case. The TS optimization is invoked by job opt_ts. Usually, the optimization performs better if a reasonable set of internal coordinates is submitted by the user (read_coord). In this case, the coordinates are stored in coord_def.inp:

1 2
2 3
2 5
2 4
4 6
6 7
1 2 5
3 2 4
4 6 7

Essentially, all bonds present in the presumed TS structure are added, further, some of the angles. The TS optimization can be invoked by typing:

explore.x ts.key 

The calculation finishes within a fraction of a second. As can be seen from the output, a small collection of files is written, containing different aspect of the obtained calculation results:

• ts_traj.xyz: Trajectory file with all structures obtained along the optimization algorithm (can be opened with a molecule viewer)
• ts_opt.xyz: The optimized structure of the transition state (can be opened with a molecule viewer)
• ts_opt.molden: Molden file containing structure and normal modes resulting from a final frequency calculation at the TS.
• explore.log: As usual, several details of the calculation (energy, gradient, convergence behavior of the optimization) are written to this file.

In the next part, a IRC reaction path will be calculated starting with the optimized TS structure. The main input file is irc.key:

pes ana_ch4oh
xyzstart ts_opt.xyz

job irc

irc {
   read_coord
   maxstep 200
   steplen 0.05  
}

The IRC calculation is ordered via job irc, starting directly from the optimized TS structure (xyzstart ts_opt.xyz), further details are defined within the irc { section. The same set of internal coordinates as for the TS optimization is read from coord_def.inp with the command read_coord. Although this is not strictly needed, the maximum number of IRC steps per direction and the length of a single step (in Angstroms) are given in the irc section as well.

Start the calculation by again invoking:

explore.x irc.key

After rapid completion of the calculation, three input files are written:

• irc.xyz: Trajectory file with the structures of the optimized IRC (from left to right)
• irc_ens.dat: Energies of the structures of the optimized IRC
• explore.log: Similar as for the TS optimization: details of the progress of the calculation

The energies within irc_ens (given in Hartrees, vs the reaction coordinate) can be plotted (here with the structures of reactants, TS and products as picture inlets):

Folder: evbopt/egrad_pes_topol

A special case for the application of explore.x is its usage with the pes topol setting. With this, no actual potential energy surface (PES) is set up but instead a pure topological analysis of the internal coordinates of the system is done, based on an extended Hückel theory-based Wiberg-Mayer bond order calculation or the comparison of the actual atomic distances to precomputed pairwise van-der-Waals radii.

All covalent bonds, angles, dihedral angles and out-of-plane coordinates are determined. If desired, also the corresponding Wilson matrix elements and Wilson matrix derivatives can be calculated.

In the example, the reaction path for a H_2 elimination from ethane is given in the file path.xyz. The keyfile topol.key contains the following lines:

pes topol

topol {
   bonds eht
   eht_scale 1.0
}

xyzstart path.xyz
job egrad

egrad {
   print_wilson 1
}

The pure coordinate evaluation mode is activated with pes topol. Note that this mode can only be used in combination with a job egrad calculation with explore.x! Further, the topol section is given, in which the extended Hückel theory (eht) based calculation of Wiberg-Mayer bond orders is called and the default cutoff for bond orders to define a bond is kept (eht_scale 1.0). The larger this value, less bonds will be detected since a larger bond order is needed (internal, a bond order of 0.5 is set to be sufficient, if the eht_scale is set to 1.0). This section is fully optional and redundant in this case, since the default settings for topol are the same. You can, however, change the method to the van-der-Waals radii comparison (bonds vdw) and set a multiplicator for those radii (vdw_scale). The larger this radii multiplicator, the more bonds will be detected.

The EHT method with default cutoff seems to work well in most cases. For very large systems, however, the EHT calculation might become expensive. There, the VDW method might be an alternative.

In the given example, the Wilson matrix elements are also calculated: print_wilson 1.

Perform the calculation with:

explore.x topol.key

The calculation should take no longer than a fraction of a second. Now, two output files were written: pes_topol.dat contains a list of all internal coordinates for each trajectory frame. The definitions of the coordinates (which types and which atoms are involved) as well as their actual values are given there, and for EHT calculations, the Wiberg-Mayer bond orders of all bonds are given as well. wilson_mat.dat contains the elements of the wilson matrix for the set of internal coordinates defined in pes_topol.dat. All those values are given in bohr and radians. If the derivatives of the Wilson matrix elements shall also be calculated, change the respective keyword to print_wilson 2. The calculation will then take approximately one second and an additional file wilson_deriv.dat containing the derivatives is written.