The following is a user documentation on (a) classical rate calculations, (b) getting the initial files for instanton optimisations, (c) optimise instantons, (d) calculate instanton rates. The procedure follows the one as implemented in the programm package DL-find. More details about the algorithms can be found here [83,36]. Note that for the interface to work properly one has to set the number of frozen atoms in odata and provide there a geomtry of one image as well to obtain the proper masses of the system. Before starting generate a executable with ./build.csh dlfoptim ifort. The ifort compiler is recommended to be used.
These files are required to run a classical rate calcualtion with CLASSICALRATES in OPTIM through an interface to the quantum rates module in DL-FIND. The output (stdout or the file arrhenius) can directly be used in an Arrhenius plot: 1000/T in Kelvin, log10 of the classical rates in s-1 (cm3s-1 in case of bimolecular rates) calculated completely classical, with quantised vibrations (which includes the zero-point vibrational energy) and including tunnelling approximatively via the simplified Wigner correction:
|κ(T) = 1 + (βωb)2 = 1 + , κ(Tc) = 1 + 2.645||(4)|
In the last column the exact analytical quantum rates for a symmetric Eckart barrier fitted to the particular system (barrier hight and ωb) are shown. All degrees of freedom perpendicular to the reaction coordinate are approximated as quantum mechanical harmonic oscillators.
KIEs can be calculated directly by first running OPTIM with the keyword CLASSICALRATES instead of INSTANTONSTARTDUMP on the Hessians for the light isotopologue. The file arrhenius of this run can be copied to the file rate_H in the directory where the rate with heavier isotopes is to be calculated. There, the same class.in as in the light case (at least the same temperature parameters) should be used. The rates obtained with the light nucleids are read and the KIE is directly calculated and written to a file called kie.
The keyword INSTANTONOPT nimageinst distortinst deltainst temperature in OPTIM specifies the optimisation of instantons.
A finite value of distortinst specifies how far the images will be spread along the unstable mode of the classical TS, see. Newton–Raphson optimisation is recommended[83,36] and currently called from OPTIM. The NR optimiser is modified to avoid convergence to higher-order saddle points . This avoids the collapse of the instanton path to the classical TS.
Instanton searches are performed in mass-weighted coordinates with masses consistent with atomic units (electron mass, me). That is, the mass of a hydrogen atom (1H) is 1837.15 me. This scales all distances up by a factor of 42.695 ( = (atomic mass unit/me)1/2) compared to mass-scaled coordinates. Thus, the tolerance criterion (tolerance) has to be smaller by the same factor to achieve equivalent convergence. A tolerance of 10-7 is usually sufficient, a tolerance of 10-8 is also often still possible. Since NR converges quadratically, the more stringent tolerance generally does not increase the number of steps dramatically.
If NR (or P-RFO) is used, the updated Hessian will be used to calculate a preliminary estimate of the rate (if qts_reactant.txt is available). In that case, qts_hessian_upd.txt will be written, which contains only the updated Hessian. qts_coords.txt will in any case be written. It acts as input for subsequent recalculation of the Hessian and a rate calculation.
Restarting instanton searches: Proper restart information (check files) is not written for the time being. Using NR, a restart is possible, though, by renaming qts_hessian_intermediate.txt (which is written after each step) to qts_hessian.txt and starting the simulation again. It will start from the Hessian and the geometry after the last full set of energies has been obtained.
Restarting of rate calculations is also only possible by using the Hessian information written for each image in qts_hessian_imageX.txt. For these files to be read, set inithessian=6. However, this requires at the moment a manual adjustment in the interface file to OPTIM.
|Pn = k Po - k + 1||(6)|
In an approximation (FPA) one can keep the instanton geometry fixed and just change the masses . This is done by calculating an instanton rate with inithessian=5 (read the Hessian from file rather than recalculating it) and changing the masses.
For lower temperature (compared to Tc) the number of images necessary can be kept at bay by adapting the integration grid (dtau) to the potential energy along the instanton path . This only makes sense if the instanton path has reached the reactant minimum. It can be achieved by setting nebk=1 in interfacetoOPTIM.f90 (this is not interpreted as the NEB force constant, but as a parameter which can vary from 0 to 1. One corresponds to a fully adaptive grid).
If KIEs should be calculated with INSTANTONRATE it is
not necessary to recalculate the Hessian for the heavier isotopologue. The
Hessians for the reactant and product are read in. If the masses provided via
the calling code (OPTIM) are different from the ones in the
Hessian files, the Hessians with the new masses will be calculated.