calc_rate calculates the rate constant (k(T)) of an elementary reaction via combined umbrella sampling and recrossing trajectory calculations. This can be done for a variety of different reaction types, which are represented via collective variables during the biased samplings. Different number of RPMD beads can be used. All calculations can be done with all potential energy surfaces available in the program package. Since this program is so far the most complex one within the Caracal program package, several examples treating numerous different reactions and situations will be covered with a large amount of detail.

A rate constant calculation with calc_rate consists of five steps:

1. Generation of initial structures for the umbrella sampling. The system is sampled successively along the predefined reaction path and structures with the ξ of the current segment are stored. All this is done classically, i.e., with one RPMD bead.
2. Equilibration and sampling on all segments of the reaction path. Bias potentials are located on the current ξ values of the reaction coordinate, for each value, an ensemble of combined equilibration/sampling trajectories are calculated. This step is the most time-consuming one.
3. Calculation of the potential of mean force (PMF) profile along the reaction path. The results of all umbrella samplings are plugged together either via umbrella integration or via the weighted histogram analysis method (WHAM).
4. Determination of the recrossing factor. A constrained is placed on the &Xi value with the largest PMF and a trajectory is calculated from which child trajectories are started. The fraction of those reaching the product side determines the resulting recrossing coefficient.
5. Calculation of the rate constant by plugging PMF difference and recrossing factor into the suitable formula.

List of examples:

• 1: H+H2 on analytical PES (simplest possible reaction on simplest PES): calc_rate/h+h2
• 2: CF3O+CH4 with the TREQ method (Typical medium sized reaction with TREQ): calc_rate/cf3o+ch4_treq
• 3: Ethane dehydrogenation (Elimination with TREQ, rather experimental): calc_rate/c2h6_elimination
• 4: Mechanochemistry (Rate of C-C bond breaking with different forces): calc_rate/c-c_mechano
• 5: H+H2 with external call (The H+H2 system again, now with external program call): calc_rate/h+h2_external
• 6: P+O2 on a custom PES (The P+O2 reaction, on a PES integrated by you into Caracal): calc_rate/p+o2_custom

Folder: calc_rate/h+h2

The first example is the simplest one to think of: The bimolecular H+H2 exchange reaction. Within EVB-QMDFF, a number of analytical potential energy surfaces were added as benchmark cases for general tests or comparisons. One of those is the H+H2 surface by Truhlar and Horowitz (link).

Only two input files are needed to run this calculation: The rate.key file containing all keywords to control the calculation and ts.xyz containing the initial structure of the system, which usually should be the TS.

The rate.key file looks as follows:

pes ana_h3 
ts_struc ts.xyz
bead_number 16
deltat 0.1

print_polymer 1.0
tdump 1.0

nvt {
  temp 300
  thermostat andersen
  andersen_step 80
}

mecha {
  type bimolec
  reactant1 1 2
  reactant2 3
  dist_inf 16
  bond_form 2-3
  bond_break 1-2
  n_paths 2
}

umbrella {
  bias 0.05 
  bonds -0.05  1.05
  dist  0.01
  gen_steps 10000
  equi_steps 10000
  sample_steps 20000
  sample_trajs 10
}

pmf {
  xi_range -0.05 1.05
  bins 5000
  method integration
}

recross {
  equi_steps 50000
  child_total 2000
  child_interval 1000
  child_perpoint 100
  child_steps 500
  mpi
}

This file contains some general keywords, which are similar to those of the dynamic program. pes ana_h3 defines the potential energy surface (PES) to use. In this case, it is the preimplemented analytical H3 surface. ts_struc reads in the name of the xyz file containing the transition state (TS) of the reaction. The number of RPMD beads is defined with bead_number, in this case, 8 beads are used. The time step deltat is given as 0.1 fs, which is quite short (since only light hydrogen atoms are present) and might be set to slightly larger values in order to make the calculation cheaper. The keywords print polymer 1.0 orders the program to print the xyz trajectory of one umbrella sampling trajectory at a certain value of the reaction coordinate ξ to file. Trajectories are not printed by calc_rate by default since those are directly evaluated internally and would result in giant files, due to the huge amount of sampled structures. The ξ=1 value is by default the TS dividing surface, therefore, structures resulting from the sampling trajectory with the bias potential located at the transition state (with respect to the reaction coordinate) are printed to the file. tdump 1.0 results in one time frame written every 1 fs simulation time, e.g., every tenth structure in this case.

Next comes the NVT section, where the thermodynamical ensemble is defined.

nvt {
  temp 300
  thermostat andersen
  andersen_step 80
}

Currently, calc_rate can only do nonperiodic calculations with the NVT ensemble (assuming an infinite volume of vacuum around the reactive system). The desired temperature temp 300 is given as well as the Andersen thermostat it apply it (thermostat andersen). The probabilistic Andersen thermostat is usually a realiable choice for rate constant calculations of small systems.

The MECHA section is of central importance, since the mechanism of the reaction to study is defined here.

mecha {
  type bimolec
  educt1 1 2
  educt2 3
  dist_inf 16
  bond_form 2-3
  bond_break 1-2
  n_paths 2
}

All following settings depend on the keyword type, which selects one of the reaction types predefined by calc_rate. Please inform me if you want to calculate a reaction which is currently not covered by the list of types! In this case, a simple hydrogen exchange shall be studied, which is covered by the bimolec type. A bimolecular reaction starts with two separated reactants. Those must be defined here by the indices of the atoms belonging to them. In this case, the initial H2 molecule contains atoms 1 and 2 (educt1 1 2), and the single hydrogen atom attacking is number three (educt2 3). The asymptotic starting point of the reaction is defined as both reactants being infinitely separated. In reality, a distance should be chosen where the interaction between the reactants can be neglected (but also not too large, since the complete reaction path must be covered with neighbored umbrella samplings; building a longer path raises the needed number of samplings). Usually, 10 to 20 Angstroms are a good choice, 16 in this case: dist_inf 16. The type bimolec requires one bond to form and one to break. Those are defined by the bond_form and bond_break keywords, after which each bond is described by the indices of the involved atoms, connected with a minus sign: bond_form 2-3. Finally, the symmetry of the reaction, i.e., the number of equivalent reaction paths, must be given with n_paths. In the case of H+H2, the hydrogen atom can attack both atoms of the hydrogen molecule, therefore leading to a twofold symmetry (n_paths 2). This number does not affect the cost of the overall rate constant calculation, since it is only a multiplicative constant in the final rate constant formula.

The UMBRELLA section defines all settings of the first and second part in the calc_rate procedure, i.e., the generation of starting structures and the umbrella equilibrations/samplings.

umbrella {
  bias 0.05 
  bonds -0.05  1.05
  dist  0.01
  gen_steps 10000
  equi_steps 10000
  sample_steps 20000
  sample_trajs 10
}

A sensitive parameter is the strength of the umbrella force constant (bias 0.05). Its value in atomic units is multiplied by the temperature of the current simulation. If the force constant is too weak, windows on regions of the PES with large slopes will be sampled insufficiently. The program will detect this and throw error messages: ERROR! The actual xi value [number] of the structure is too much away from the ideal value [number]! The calculation will not terminate at this point but try to start a new trajectory. Only if the maximum number of allowed errors during the current calc_rate calculation (1000 as default) is reached, the program is terminated. If the force constant is too strong, where will be gaps between the distributions resulting from single samplings and the PMF curve will look noisy (if the value is absurdly large, the whole system might break during the dynamics, resulting in errors as well). Values between 0.02 and 0.1 usually proved to be reasonable. In general, the reaction path variable ξ is defined such that the reactants are located at ξ = 0 and the TS is located at ξ = 1. Therefore, the reaction path should at least be sampled between those borders. It is, however, strongly recommended to somewhat extend the treated range somewhat beyond those values, for example by 0.05 in both directions. Then, the full range of the reaction path is defined by bonds -0.05 1.05. The related keyword 'dist 0.01' defines that umbrella samplings will be performed in increments of Δ ξ = 0.01, resulting in 111 umbrella windows. During the generation of initial structures (part 1), 10000 MD steps are sampled for each window (gen_steps 10000). In part 2, the umbrella sampling, 10 trajectories are sampled for each window, where for each trajectory, starting from the classical initial structure, 10000 MD steps are calculated for equilibration with 8 RPMD beads (equi_steps 10000), followed by 20000 MD steps of actual samplings, whose results are stored for the subsequent part (sample_steps 20000).

The PMF section orders how the potential of mean force along the defined reaction coordinate shall be calculated.

pmf {
  xi_range -0.05 1.05
  bins 5000
  method integration
}

Similar to the umbrella sampling part, a range of the reaction coordinate ξ must be given, which should not larger than that within the umbrella sampling part (being the same range should work usually): xi_range -0.05 1.05. The PMF is calculated on 5000 bins along the coordinate range (bins 5000) with the umbrella integration method (method integration). This method is usually preferred for small systems, only for larger ones (with shorter samplings affordable), the more expensive WHAM method might be used, which needs less amounts of data to result in a somewhat reasonable curve.

Finally, the calculation of the recrossing coefficient κ is defined in the RECROSSING section.

recross {
  equi_steps 50000
  child_total 2000
  child_interval 1000
  child_perpoint 100
  child_steps 500
  mpi
}

A recrossing calculation consists of two parts: First, a parent trajectory is calculated, constrained to the transition state dividing surface (located at the maximum of the PMF curve within ξ). From that trajectory, short child trajectories are started in pairs of opposite randomly sampled initial momenta and evolved over short times without thermostatting. The recrossing (essentially the probability of a system arrived at the transition state dividing surface to react to the products) is proportional to the fraction of child trajectories ending on the product side. equi_steps states the number of MD steps the parent trajectory shall be equilibrated at the beginning. The total number of child trajectories started along the whole calculation is given by child_total 2000. Between two crowds of spawned child trajectories, the parent trajectory shall be propagated another 1000 steps (child_interval 1000), and for each dump of children, 100 of them will be started (child_perpoint 100). The propagation of child trajectories can usually be quite short, here, only 50 fs/ 500 MD steps (child_steps 500). The recrossing calculation is parallelized with MPI (mpi).

Especially the umbrella sampling part can take quite some time (at least for larger systems than H+H2), therefore it is parallelized with MPI as well. At least two cores must be used, since one core acts as master that collects the results, where all other cores propagate one umbrella trajectory each. The calculation can be called like this:

mpirun -np [number of cores] [path to bin]/calc_rate.x rate.key 

With 6 cores, the calculation roughly takes 4 minutes. After finishing, the following information is printed out:

 Type of reaction: BIMOLECULAR EXCHANGE

- The rate formula is:

k(T)=n_paths*kappa*4*pi*dist_inf**2*sqrt(1/(2*pi*beta*my_R))*exp(-beta*(W(xi_TS)-W(0)))
 where my_R=m1*m2/(m1+m2)

- The parameters are defined as (in atomic units):

  * n_paths =   2
  * kappa =   0.76669170
  * dist_inf =     30.23561799 bohr
  * beta =   1052.58224824 1/hartree
  * m1 =   3674.30529123 amu
  * m2 =   1837.15264562 amu
  * my_R =   1224.76843041 amu
  * W(xi_TS) =      0.01847336 hartee at xi=   0.99905981
  * w(0) =      0.00017084 hartree at xi=   0.00006001
  * deltaW =     48.05326811kJ/mol

----------------------------------------------
 The value of the reaction rate constant is:
     9.81630659E+07  cm^3/(mol*s)
     1.63003578E-16  cm^3/(molec*s)
----------------------------------------------

 Timings: 
 ----------
Read in of settings and initialization:         0.002 s.
Initial structure generation:                   2.147 s.
Umbrella equilibration and samplings:         176.160 s.
Evaluation of sampling results:                 0.003 s.
Calculation of the free energy surface:         0.019 s.
Recrossing factor calculation:                 35.828 s.
Final rate constant calculation:                0.000 s.

The calculation needed a total time of        214.160 seconds.
These are:      0 days,  0 hours,  3 minutes and 34 seconds.

 -.-. .- .-. .- -.-. .- .-.. 
Informations printed to calc_rate.log!
Other results are stored into folder 300K_16bead                                       
CARACAL rate constant calculation successfully finished!

The calculated rate constant agrees well with the value reported in Suleimanov et al..

Further details about the calculation can be found by looking into the 300K_16bead folder. The following files and folders can be found here:

• calc_rate_traj.xyz: only for visualization: an example trajectory, which was ordered to be printed with the print_polymer command. Here, one umbrella sampling trajectory at the TS region (xi=1.0) is printed. We can see the movement of the RPMD beads when opening it.
• equilibrated_ens.dat: The energies of the structures generated in the inital structure generation phase. This can serve as a debug for a newly built PES. If some energies are huge or strong fluctuations occur, something might be wrong.
• equilibrated_struc.xyz: The structures generated in the first part, their energies printed to equilibrated_ens.dat.
• pmf_integration.dat: One of the most interesting files, the calculated potential of mean field/free energy curve of the reaction, which is used for the rate constant calculation. It looks something like this:

If the curve looks very noisy, the longer samplings are needed. If several kinks or local minima are present, the PES might be bad.

• recrossing_time.dat: The time-dependent transmission coefficient kappa. It enters as a prefactor into the k(T) calculation and is determiend by the recrossing calculation. If this curve is very noisy, more child trajectories are needed, if it does not converge into the plotted area, longer child trajectories are needed.

• statistics: This folder basically contains the results of the umbrella sampling trajectories, when the umbrella integration method is used for evaluation. Averages and variances of the constructed reaction coordinate are printed here and are read in, if the calculation shall be restarted or continued after a premature canceling (continuation is possible, all completed windows are skipped in the new calculation).
• xi_distributions: This folder is empty and only filled, if the PMF profile shall be determined with the WHAM method.

Folder: calc_rate/h+h2

The potential energy surface of this reaction was created in a black box manner using the black_box program. In this program, the Transition Region Corrected Reaction Path Empirical Valence Bond Quantum Mechanically Derived Force Field (TRC-RP-EVB-QMDFF) (TREQ) method (link) is used to describe the potential energy surface. Its parameters are specified in the second section of the rate.key file, after the TS structure and the number of beads:

pes treq

qmdff {
   ffnames educts.qmdff products.qmdff
   eshift -452.74250344759     -452.7631216076
}

treq {
  points 28
  irc_struc irc.xyz
  irc_ens irc_ens.dat
  rp_exp_coeff 35.0
  rp_evb_mid 0.6  0.1
  read_coord
}

As we can see, TREQ is essentially an EVB-QMDFF method, therefore, it needs two QMDFFs, one for the reactant side (educts.qmdff) and one for the products side (products.qmdff) of the relevant part of the potential energy surface. The energy shift of both QMDFFs is given as well: eshift -452.74250344759002 -452.76312160763899, within black_box the energies are also calculated with the higher level method, if a dual-level approach is desired.

In the treq { block itself, parameters describing the transition between the two QMDFFs are specified. points 28 gives the number of gradient+Hessian reference points along the reaction path. The related data is stored within the grad_hess.dat file. The structures and energies of the reaction path that serves as foundation of the TREQ PES are were stored in the files irc.xyz and irc_ens.dat; the file names can be defined with the keywords irc_struc' and irc_ens`.

The keywords rp_exp_coeff 35.0 and rp_evb_mid 0.6 0.1 are somewhat more technical. The first one is the coefficient of a multiplied Gaussian function that damps the coupling orthogonal to the path. The larger that value, the faster the coupling goes to zero. The second one defines width of the direct interpolation region along the path (in units of the reaction coordinate s, which goes from 0 to 1. The direct interpolation interval is always centered on the transition state, in this case, it covers 60% of the reaction path, where the outer 10 % (5% in both directions) are smooth transition regions towards usual RP-EVB coupling. Fortunately, these parameters can usually be kept at values similar to those listed here, their ideal values seem to be quite independent of the actual reaction the TREQ term shall describe.

Finally, the keyword read_coord again orders the program to read in a suitable set of internal coordinates for description of the TREQ PES from the file coord_def.inp. In this case, the black_box program generated the coordinates automatically by comparing the internal coordinates of both QMDFF's (essentially, all bonds or angles occurring in at least one of the QMDFFs is added to the set, a procedure that can of course be done manually as well).

The remaining settings are mostly similar to those of the first example. Most relevant differences occur in the mecha { section. Here, both reactants and the formed and broken bonds must be defined depending on the indices of the atoms. Further, there are four equivalent hydrogen atoms that could be abstracted by the CF3O molecule, therefore, the symmetry of the reaction is 4.

mecha {
  type bimolec
  educt1    1   2   3   4   6
  educt2    5   7   8   9  10
  dist_inf 20.0
  bond_form  7-6
  bond_break  6-1
  n_paths 4
}

The calculation can again be invoked with MPI, where the umbrella samplings and the recrossing trajectories are parallelized, the number of cores could be raised up to the number of sampled umbrella windows.

mpirun -np [number of cores] [path to bin]/calc_rate.x rate.key 

Depending on the PC and the number of cores, the calculation time might fluctuate considerably. If six cores on an usual office PC are used, it should take roughly 20 minutes. After the live-output during the sampling trajectories, the following information should be printed out:

Type of reaction: BIMOLECULAR EXCHANGE

- The rate formula is:

k(T)=n_paths*kappa*4*pi*dist_inf**2*sqrt(1/(2*pi*beta*my_R))*exp(-beta*(W(xi_TS)-W(0)))
 where my_R=m1*m2/(m1+m2)

- The parameters are defined as (in atomic units):

  * n_paths =   4
  * kappa =   0.77838284
  * dist_inf =     18.89726125 bohr
  * beta =    789.43668618 1/hartree
  * m1 =  29223.27239962 amu
  * m2 = 154927.50125870 amu
  * my_R =  24585.77002710 amu
  * W(xi_TS) =      0.01553261 hartee at xi=   0.99326665
  * w(0) =      0.00101985 hartree at xi=   0.00007602
  * deltaW =     38.10323776kJ/mol

----------------------------------------------
 The final value of the rate coefficient is: 
     4.93673769E+10  cm^3/(mol*s)
     8.19764440E-14  cm^3/(molec*s)
----------------------------------------------

 Timings: 
 ----------
Read in of settings and initialization:         0.279 s.
Initial structure generation:                  30.727 s.
Umbrella equilibration and samplings:        1012.140 s.
Evaluation of sampling results:                 0.003 s.
Calculation of the free energy surface:         0.015 s.
Recrossing factor calculation:                236.884 s.
Final rate constant calculation:                0.000 s.

The calculation needed a total time of       1280.048 seconds.
These are:      0 days,  0 hours, 21 minutes and 20 seconds.

.. .----. -- -.. --- -. . 
Informations printed to calc_rate.log!
Other results are stored into folder 400K_8bead                                        
Rate constant calculation successfully finished!

If we compare the calculated rate constant with the experimental value at 400 K (8.01E-14, NIST-Database, paper), an almost perfect agreement was achieved with this calculation! Since the velocities are initialized by chance and the settings are still rather cheap, the actual value will fluctuate somewhat in each run. More rigorous samplings and rate constant calculations at different temperatures on this TREQ surface, however, revealed that the experimental rates and Arrhenius parameters can indeed be reproduced quite well with this particular PES.

During the calculation, probably a number of error notifications will be dropped, e.g.:

ERROR! The energy of the actual structure (  -3589.4973010392500       Eh) is much
higher than the maximum allowed energy! (  -3609.6678390537641      Eh!
Add the keyword RPMD_EN_TOL to change the tolerance energy per bead (kJ/mol)
Actual value:   4000.0000000000000      kJ/mol
The error occured after        2436 verlet steps

Depending on the size of the set of internal coordinates in usage, some dozens of those errors might occur during the umbrella sampling part. Since, however, the trajectories are completely restarted from the presampled structures if an error occurs, it is highly unlikely that they influence the final rate constants. Their origin probably lies within the projection algorithm of the actual structure to the reaction path, which at some time introduces discontinuous steps in the potential energy surface if the nearest point on the reaction path changes significantly between two MD steps. The resulting large gradient components might then lead to fast movements and explosions. These problems will (hopefully) be solved with patches in the near future.

Folder: calc_rate/c2h6_elimination

The PES of the next example is also described with the TREQ function. The elimination of hydrogen from Ethane, forming ethylene, is of course one of the most important processes for industry. Usually it will be conduced with the assistance of a catalyst. Here, we look at the simple gas phase reaction. It is a unimolecular reaction. Currently, no exact prefactor of unimolecular reaction rates (which depends on the actual mechanism!) is implemented, therefore, the simple TST formula will be used. This is an approximation and might lead to errors of factor 2 to 4 in the calculated rates!

The keyfile for this calculation is elimination.key. Its structure is mostly similar to the one of the last calculation (since TREQ is used here as well).

• bead_number 4: We calculate the reaction at a rather high temperature (in order to have experimental values), therefore, tunneling effects will be almost negligible and a smaller number of RPMD beads is sufficient.
• type elimination: The elimination mechanism is activated, which expects two bonds to be broken (C-H) and one to be formed (H-H).
• reactant1 1 2 3 4 5 6 7 8: All atoms of the system must be listed here, since all are part of the reactant
• bond_form 7-8: The H-H bond will be formed
• bond_break 7-1 8-4: The indices of the two broken bonds
• n_paths 9: The number of equivalent reaction paths. In this case, there are nine different combinations of hydrogen atoms that could be eliminated (three per carbon atom).
• reactants_struc reactants.xyz: The structure of the reactants minimum. It is needed to setup the collective variables for the reaction path (interpolation between reactants and TS).
• bias 0.20: The force constant is set to a larger value than usual since the barrier is quite high (ca. 200 kJ/mol).
• bonds -0.05 1.10: The umbrella sampling coordinate range was slightly increased, since the overall reaction path is quite short and possible displacements of the PMF maximum will become larger in terms of the reaction coordinate.

The keyword read_coord is missing, since the internal coordinates for the TREQ term are now not read in from file but determined automatically by Caracal when setting up the system.

Invoke the calculation as usual by typing:

mpirun -np [number of cores] [path to bin]/calc_rate.x elimination.key 

The calculation will take some minutes, the final output should look something like this:

...
Type of reaction: Elimination (two bonds broken and one built).
This is an UNIMOLECULAR REACTION!

So far, no general rate constant formula has been implemented for this case.
For the different mechanisms, different prefactors will be developed in the 
 in the next version of CARACAL.
Until then, the simple TST rate formula (with kappa) will be used:

k(T)=kappa*n_paths*k_B*T/h*exp(-(W(xi_TS)-W(xi_min))/(k_B*T)) 
 where T=temperature, k_B=Boltzmann constant, h=Plancks constant
  W(xi_TS), W(xi_min)=free energies at the TS and the reactants minimum.
- The parameters are defined as (in atomic units):

  * kappa =   0.43495760
  * n_paths =   9
  * T =    873.00000000 K
  * W(xi_TS) =    275.09748003 kJ/mol at xi=   0.98764153
  * W(xi_min) =      0.00000000 kJ/mol at xi=   0.51490698

----------------------------------------------
 The final value of the rate coefficient is: 
     2.47071734E-03  s^(-1)
----------------------------------------------
...

If we compare this value with the experiments on the (NIST database), the calculated value is roughly an order of magnitude larger than the experimental one, which might be due to the QM reference, the rather short samplings or the non-exact prefactor in the rate constant equation.

For further analysis, PMF profile and time-dependent recrossing factor can be plotted. The recrossing-time profile looks quite common (convergent behavior), the PMF profile, however, looks quite special if we compare it to the profile from example 1. Besides the broad maximum at the TS (due to the shorter reaction path), the minimum is not located at Xi=0 (reactants minimum) but somewhere in between! By looking at the equilibrated_struc.xyz file, the structures at xi=0.5 still look rather "unbiased", therefore, the left maximum might correspond to an arrangement of the hydrogen atoms (torsional landscape) not optimal for elimination.

Another feature of Caracal is the possibility to add external forces to MD trajectories, both in dynamic (see Disulfide Mechanochemistry as well as in calc_rate. We can see how the rate constant of a calculation changes if two force vectors (of opposite directions, in order to avoid translations of the system) are applied to two arbitrary atoms in the system of interest. In this exemplaric case, which might serve as template for more realistic applications, we study the dissociation of a C-C bond in propane under influence of different external forces.

The molecule is described by a single QMDFF, since QMDFF uses Morse functions as bond potentials and is thus able to describe dissociations (semireactive force field).

Main input is the rate.key file, which now has an additional force section, activated by the add_force command in the header.

force {
  vec1 1 1E-9 -6.3 2.5 0
  vec2 8 1E-9 6.3 -2.5 0
}

Force vectors are applied to the atoms with indices 1 and 8 (first numbers after the vec1/2 keywords). The magnitude is 1 nN (1E-9 N), the directions (in cartesian coordinates) are -6.3 2.5 0 and 6.3 -2.5 0 (antiparallel).

The reaction mechanism is a decomposition, taking place along one bond (decom_1bond):

mecha {
  type decom_1bond
  reactant1  1 2 3 4 5 6 7 8 9 10 11
  bond_form
  bond_break  1-5
  dist_inf 10
  n_paths 1
  reactants_struc minimum.xyz
}

As in the elimination example, all atoms belong to one reactant (propane). The C-C bond between atoms 1 and 5 is broken; the external force is applied on the atoms 1 and 8 (on the outer atoms)! See the next picture.

Since recrossing calculations are pointless with bond breakings, where no real barrier exist (asymptotic reaction), the recorossing settings were reduced to a minimum. The kappa coefficient will almost be zero in most cases and be replaced by a default value of 0.5.

In order to study the influence of the external force on the calculated rates and PMF profiles, the script start.sh was prepared, which automatically starts calc_rate calculations of forces between 0 and 3 nN. Open the script and enter your path to calc_rate and MPI there! Then, start the script:

bash start.sh

The total calculation might take half an hour, then, folders with results for seven different applied forces are generated. First, we can compare the calculated rate constants (your results will be slightly different, due to the randomized samplings):

• 0 nN: 1.02 E-70 1/s
• 0.5 nN: 4.54 E-47 1/s
• 1 nN: 3.02E-32 1/s
• 1.5 nN: 7.00E-21 1/s
• 2.0 nN: 8.77E-12 1/s
• 2.5 nN: 3.96E-4 1/s
• 3.0 nN: 6.06E+2 1/s

The rate constants get significantly larger with stronger force attached, which is of course no surprise. As it is for a rope, the stronger the force attached to it, the faster it will be broken. We can compare the calculated PMF profiles as well, they were copied into the main folder by start.sh and are named pmf_0.0nN.dat to pmf_3.0nN.dat. The combined plot will look something like this:

Folder: calc_rate/h+h2_external

This example uses the same methodology (and the same PES) as the related example in the dynamic program: dynamic/h3_external Each time a gradient is needed for a structure (during equilibration, one bead of one of the umbrella trajectories or during the recrossing calculation), the external program is called and the resulting gradient is read in again. Note, however, that this method of external calling is quite inefficient. Especially if very fast programs or routines shall be connected, the needed time for file IO will be much longer than the actual external calculation time! If, e.g., a publicly available analytical PES shall be used with Caracal that is not already implemented, please contact me! I might add this routine to the program package itself.

Here, the H+H2 system already presented in example 1 is calculated again. It uses the well known H+H2 function by Truhlar and Horowitz in an updated version (DOI).
Since this surface is already implemented within Caracal (ANA_H3), this is rather pointless and shall mainly serve as a tutorial how to connect an external program/subroutine to Caracal.

The rate.key file looks essentially as in example 1:

pes external

ts_struc ts.xyz
rpmd_beads 8
deltat 0.1

print_polymer 1.0
tdump 1.0

nvt {
  temp 300
  thermostat andersen
  andersen_step 80   
}

mecha {
  type bimolec
  reactant1 1 2
  reactant2 3
  dist_inf 16
  bond_form 2-3
  bond_break 1-2
  n_paths 2  
}

umbrella {
  bias 0.05
  bonds -0.05  1.05
  dist  0.01
  gen_steps 10000
  equi_steps 10000
  sample_steps 20000
  sample_trajs 10  
}

pmf {
  xi_range -0.05 1.05
  bins 5000
  method integration
}

recross {
  equi_steps 50000
  child_total 10000
  child_interval 1000
  child_perpoint 100
  child_steps 500
  mpi
}

external {
   symlink /home/trebonius91/caracal/calc_rate/h3_external/external.x
}

Only two changes were made: First, the pes external keyword activates the external call.
Second, the external section with the symlink keyword was added. In this case, the external program is called with the command /home/trebonius91/caracal/calc_rate/h3_external/external.x (change this to your own path!) Note: in contrast to the dynamic.x program, where a relative path is valid as well, an absolute path must given here, since the program will also be called in several subfolders for the equilibration and umbrella trajectories!

Details concerning the Fortran90 program can be found in the dynamic section: dynamic/h3_external

In the current example, the external program (call_external.f90 + h3_external.f) must be compiled first, e.g., with gfortran:

gfortran -o external.x call_external.f90 h3_external.f

Then, the calc_rate.x program can be executed as usual:

mpirun -np [number of cores] [path to bin]/calc_rate.x  rate.key

This calculation will take orders of magnitude longer than that in example 1! Therefore, this external call method should only be used for large systems or expensive neural network potentials, where the overhead for interaction with the external program is small compared with the actual external energy/gradient calculation.

Folder: calc_rate/po2_custom

In this example, the addition of a phosphorus atom to an oxygen molecule to form a PO2 molecule is studied. Similar to example 5, no integrated PES of Caracal is used for that, but instead a PES from the literature (Chen et al.). This time, however, the PES is included directly into the Caracal code such that the whole process of calling it is much more efficient, especially since this PES routine only calculates the energy of a structure such that the gradient must be calculated numerically!

This tutorial uses the same custom PES as example 7 of the dynamic wiki page. If you already did this and included the code into your local Caracal, you can directly jump to step 4 of this example.

The following steps need to be taken:

1: Add the PES source code file to the Caracal code

Copy the file pes_po2.f to the src folder of Caracal and add it to the Makefile.

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

to

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   pes_po2.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

2: Implement an energy/gradient call to the routine

Open the pes_po2.f file. The main subroutine is placed directly at the top:

      SUBROUTINE PES_PO2(R1,R2,R3,ENER)
c       R1=roo;R2=rpo;R3=rpo
      IMPLICIT  REAL * 8 (A-H,O-Z)
      DOUBLE PRECISION :: R1,R2,R3,V3
      DOUBLE PRECISION, DIMENSION(3) :: R
      DOUBLE PRECISION :: E12,E13,E23


      R(1)=R1
      R(2)=R2
      R(3)=R3
      CALL CHIPR_O2(R1,E12)
      CALL CHIPR_PO(R2,E13)
      CALL CHIPR_PO(R3,E23)
      CALL TRIAT(R1,R2,R3,V3)
      ENER=E12+E13+E23+V3
      RETURN
      END

It can be seen that the subroutine requires three interatomic distances as input (R1, R2 and R3) and gives the current energy as output (ENER). Thus, there is no gradient calculated. Further, it does not accept a xyz structure as provided by Caracal as input! Therefore, the subroutine need to be changed somewhat into the following (one possibility):

      SUBROUTINE PES_PO2(xyz,ENER)
c       R1=roo;R2=rpo;R3=rpo
      IMPLICIT  REAL * 8 (A-H,O-Z)
      DOUBLE PRECISION :: R1,R2,R3,V3
      DOUBLE PRECISION, DIMENSION(3) :: R
      DOUBLE PRECISION, DIMENSION(3,3) :: xyz
      DOUBLE PRECISION :: E12,E13,E23

      R1 = sqrt(dot_product(xyz(:,2)-xyz(:,3),xyz(:,2)-xyz(:,3)))
      R2 = sqrt(dot_product(xyz(:,1)-xyz(:,2),xyz(:,1)-xyz(:,2)))
      R3 = sqrt(dot_product(xyz(:,1)-xyz(:,3),xyz(:,1)-xyz(:,3)))

      R(1)=R1
      R(2)=R2
      R(3)=R3
      CALL CHIPR_O2(R1,E12)
      CALL CHIPR_PO(R2,E13)
      CALL CHIPR_PO(R3,E23)
      CALL TRIAT(R1,R2,R3,V3)
      ENER=E12+E13+E23+V3
      RETURN
      END

First, the arguments of the subroutine were changed, R1, R2, R3 were replaced by xyz. This array is then defined in the declaration section: DOUBLE PRECISION, DIMENSION(3,3) :: xyz. The first dimension always are the x-, y- and z coordinates of the current atom, the second dimension is the index of the atom. Now, the needed interatomic distances must be calculated from the given structure. This is done directly after the variable declaration. According to the given comment describing the distances (with P being the first atom), the O-O distance is calculated by: R1 = sqrt(dot_product(xyz(:,2)-xyz(:,3),xyz(:,2)-xyz(:,3))) and so on.

Then, open the file custom_grad.f90 and add the gradient routine. There, add the subroutine to the code snippet for numerical gradients:

!   call your_routine(xyz2,e_evb)

   step=num_grad_step
   do i=1,natoms
      do j=1,3
         do k=1,2
            if (k.eq.1) then
               xyz2(j,i)=xyz2(j,i)+step
            else if (k.eq.2) then
               xyz2(j,i)=xyz2(j,i)-2*step
            end if
!
!     Activate this call here and adapt it to your routine!
!
!               call your_routine(xyz2,e_tmp)


            if (k.eq.1) then
               e_upper=e_tmp
            else if (k.eq.2) then
               e_lower=e_tmp
            end if
         end do
         g_evb(j,i)=(e_upper-e_lower)/(2*step)
         g_evb(j,i)=g_evb(j,i)*hartree*bohr
         xyz2(j,i)=xyz2(j,i)+step
      end do
   end do
end if

Two calls can be seen here. The first call calculates the energy of the undistorted structure (e_evb). It must be changed accordingly:

call pes_po2(xyz2,e_evb)

The second call inside the do-loop calculates the actual energy for each geometry distortion and stores it into e_tmp. Those values are then postprocessed to obtain the gradient components. The call need to be changed to:

call pes_po2(xyz2,e_tmp)

The subroutine works with atomic units, therefore, no unit conversions need to be done here.

3: Recompile Caracal

After all changes are implemented, recompiling the code and making the executables should be no problem:

make 

4: Run the calculation

Now, the actual rate constant calculation can be started. Open the rate.key file to see the settings for the calculation (only upper half shown here):

 pes custom 

 ts_struc ts.xyz
 rpmd_beads 1
 deltat 0.2
 num_grad

 print_polymer 1.0
 tdump 1.0

 custom {
   pes_number 1
 }

 nvt {
   temp 500
   thermostat andersen
   andersen_step 80
 }

 mecha {
   type addition
   reactant1 1
   reactant2 2 3
   dist_inf 8
   bond_form 1-2 1-3
   bond_break 2-3
   n_paths 1
 }

We use a custom PES that was built into Caracal manually (pes custom). Since this PES has no analytical gradient, it must be calculated numerically (num_grad). It is assumed that this PES is the first one you added to Caracal, therefore, the keyword pes_number 1 is given into the custom section. If you already did example 6 of the dynamics program, you need to alter this to pes_number 2 since then the O4 PES is the first one in your local list.

During the reaction, the P atom will insert into the O-O bond of the oxygen atom. Formally, this can be seen as an addition, thus type addition, with reactant as well as the forming and breaking bonds is given accordingly.

Now start the calculation as usual by typing

mpirun -np [number of cores] [path to bin]/calc_rate.x rate.key 

The calculation should take roughly 10 minutes on 6 cores of a local machine and gives the usual output (see other examples above).