.. tip:: All input files can be downloaded: :download:`Files `. pifep ========= .. contents:: :local: This option controls the path integral free-energy perturbation (PI-FEP) simulation. Options ---------- .. option:: temp .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``298.15`` Simulation temperature in Kelvin. .. option:: sampling .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - ``bisection`` Bisection beads sampling algorithm. * - Default - ``bisection (must be provided)`` Define the beads sampling method. Chin approximation scheme, ``ca``, does not support bisection sampling algorithm. .. option:: bisection_level .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``6`` The number of PI beads under bisection algorithm which is defined as :math:`2^N`, where :math:`N` denotes bisection level. This keyword must be used in conjunction with ``bisection`` sampling algorithm. .. option:: classical_coord_file .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A filename * - Default - None Filename of an ``.xyz`` file which stores classical nuclear coordinates. .. option:: checkpoint_fn .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A filename * - Default - ``chk`` Checkpoint filename. .. option:: restart .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A filename * - Default - None The checkpoint filename from which the simulation resumes. .. option:: num_pifep_samples .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``100`` The number of PI beads samplings to perform for each classical nuclear coordinate. .. option:: quantised_atoms_list .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - List of integers or range of integers * - Default - None (meaning all atoms are quantised) The indices of atoms which are quantised to PI beads. .. option:: isotope_indices .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - List of integers or range of integers * - Default - None The indices of atoms which are specified for isotope exchange. .. option:: isotope_num .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - List of integers * - Default - None The isotope number based on table below (the atomic masses are hardcoded for now) .. table:: Isotopes table :widths: auto +---------------------------+-----------------+ | Isotope Name | isotope_num | +===========================+=================+ | :math:`^1\mathrm{H}` | 1 | +---------------------------+-----------------+ | :math:`^2\mathrm{H}` | 2 | +---------------------------+-----------------+ | :math:`^3\mathrm{H}` | 3 | +---------------------------+-----------------+ | :math:`^{12}\mathrm{C}` | 1 | +---------------------------+-----------------+ | :math:`^{13}\mathrm{C}` | 2 | +---------------------------+-----------------+ | :math:`^{14}\mathrm{C}` | 3 | +---------------------------+-----------------+ | :math:`^{14}\mathrm{N}` | 1 | +---------------------------+-----------------+ | :math:`^{15}\mathrm{N}` | 2 | +---------------------------+-----------------+ | :math:`^{16}\mathrm{O}` | 1 | +---------------------------+-----------------+ | :math:`^{17}\mathrm{O}` | 2 | +---------------------------+-----------------+ | :math:`^{18}\mathrm{O}` | 3 | +---------------------------+-----------------+ .. option:: write_log .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A file name * - Default - ``log.dat`` Log filename which stores summarised simulation output. .. option:: random_seed (optional) .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integers * - Default - None A seed value used to initialize the pseudo-random number generator (PRNG). Theoretical Background ---------------------------------- Path integral free-energy perturbation (PI-FEP) is an efficient method for computing kinetic isotope effects (KIEs) of chemical reactions. It can be used in conjunction with **QM** methods (e.g., DFT, AM1, xTB), **MM** methods (e.g., **CHARMM**), and, most importantly, combined QM/MM methods. The PI-FEP module within **Qbics** is implemented based on the works below: - Major, D. T.; Gao, J. `Implementation of the bisection sampling method in path integral simulations. `_ *J. Mol. Graphics Model.* **2005**, *24*, 121-127. - Major, D. T.; Garcia-Viloca, M.; Gao, J. `Path Integral Simulations of Proton Transfer Reactions in Aqueous Solution Using Combined QM/MM Potentials. `_ *J. Chem. Theory Comput.* **2006**, *2*, 236–245. - Major, D. T.; Gao, J. `An Integrated Path Integral and Free-Energy Perturbation-Umbrella Sampling Method for Computing Kinetic Isotope Effects of Chemical Reactions in Solution and in Enzymes. `_ *J. Chem. Theory Comput.* **2007**, *3*, 949–960. Input examples --------------------------- The classical reaction path is obtained using classical MD simulation *a priori*. The reaction coordinates in Cartesian coordinates are stored in a ``.xyz`` file and provided using the ``classical_coord_file`` keyword. The molecular definition of classical coordinates provided in the ``.xyz`` file must be identical to that in the ``mol`` section. In the examples below, proton transfer reactions of carbon acid will be studied using PI-FEP simulation. Reactions in the gas phase and in aqueous solution are given in Examples 1 and 2, respectively, for teaching purposes. Two examples used NDDO and NDDOMM methods, HF/DFT, MM, and xTB methods are also supported and used similarly. Example 1: Proton/Deuterium transfer of carbon acid using AM1 model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We begin by assuming the classical mechanical PMF (CM-PMF) for the reaction has already been determined. The CM-PMF is then divided into small windows with a number of classical nuclear coordinates. For demonstration purposes, here we use the transition state window, which contains only 2 :download:`classical nuclear coordinates `. .. code-block:: bash :linenos: :caption: :download:`pifep-ts.inp ` nddo charge -1 spin2p1 1 end mol C -0.91443957 -0.60832058 0.13234559 C 0.56752523 -0.88854599 -0.01251568 C -1.38098404 0.57350638 -0.39929395 C -1.86640306 -1.49230171 0.66658200 C 1.40989454 0.39927020 0.10774010 H 0.90764801 -1.61825893 0.79377293 H 0.75135566 -1.47355548 -1.02103320 C -2.81217464 0.76351767 -0.60469285 O -0.61312320 1.51868550 -0.90259881 C -3.23488869 -1.24001922 0.62699159 H -1.41778574 -2.42988305 1.05788219 H 0.47186573 1.21971490 -0.39894757 C 1.58524382 0.80585715 1.58203285 N 2.80702371 0.28695417 -0.48459260 C -3.72646798 -0.05196086 0.02804376 H -3.05945315 1.70341266 -1.14865278 H -3.95495955 -2.00324151 0.97518366 H 2.50342196 1.42577625 1.79057995 H 0.70024044 1.30710258 1.99617930 H 1.82497321 -0.12313539 2.15506963 O 3.72131801 0.73970515 0.16512141 O 3.08057940 -0.23586136 -1.53039467 H -4.79156884 0.24245093 0.02696987 end pifep sampling bisection bisection_level 3 # 2^3 beads are used classical_coord_file ts2.xyz checkpoint_fn chk #restart chk num_pifep_samples 50 quantised_atoms_list 1-3 5 8 9 12-14 # -O(Ph)-H-C(donor)- and 3 adjacent atoms isotope_indices 12 # 12th atom, H, is exchanged to D isotope_num 2 # 1H->2D write_log cc2-s50.log end task pi am1 end By default, if one leaves ``quantised_atoms_list`` empty, all atoms will be quantized into PI beads. Consequently, the non-reactive atoms will also contribute to the quantum canonical partition function, and therefore contribute to statistical uncertainty. We can limit the quantization to atoms that are only relevant to the reactive region, **i.e.** the donor C, the migrating proton, the acceptor O(Ph), and the three atoms adjacent to the primary isotopic substitution site are quantized to PI beads. In PI-FEP theory, the quantum mechanical potential of mean force (QM-PMF) has the following expression: .. math:: \begin{align} W_{\mathrm{QM} } &= W_{ \mathrm{CM} } + \Delta G_{\mathrm{QM} } \\ &= W_{ \mathrm{CM} } - \beta ln \langle \delta(\bar{z}) \langle e^{-\beta \Delta \bar{U} } \rangle_{\mathbf{FP} } \rangle_\mathrm{CM} \end{align} where :math:`\beta` denotes reciprocal temperature times Boltzmann constant. :math:`\bar{z}` represents the centroid of PI beads being constrained to the corresponding classical nuclear coordinate. The outer average is based on the number of classical nuclear coordinates contained in the ``.xyz`` file (specified under ``classical_coord_file``), and the inner average is based on the number of samplings of PI beads specified under ``num_pifep_samples``. The ratio of quantum partition functions of the two isotopes, i.e. light(:math:`\mathrm{L}`) and heavy(:math:`\mathrm{H}`), at a given state, :math:`\bar{z}`, can be written as a free energy perturbation: .. math:: \begin{align} e^{-\Delta G^{\mathrm{H} \rightarrow \mathrm{L}}_{\mathrm{QM}}} = \frac{Q^{\mathrm{L}}_{\mathrm{QM}}(\bar{z})}{Q^{\mathrm{H}}_{\mathrm{QM}}(\bar{z})} = \frac{ \langle \delta( \bar{z} ) \langle e^{-\beta \Delta \Delta \bar{ U }^{ \mathrm{H} \rightarrow \mathrm{L} } } e^{-\beta \bar{ U }^{ \mathrm{H} } } \rangle_{ \mathrm{FP} } \rangle_{\mathrm{CM}} }{ \langle \delta( \bar{z}) \langle e^{-\beta \bar{U}^{\mathrm{H}} } \rangle_{\mathrm{FP}} \rangle_{\mathrm{CM}} } \end{align} where :math:`\Delta \bar{U} = \frac{1}{P} \sum_{k}^{P} V_k^{\mathrm{QM}} - V^{\mathrm{CM}}`. This method allows the ratio of partition functions to be computed within a single PI-FEP simulation, which avoids statistical uncertainties that would arise from conducting two independent simulations. In the current reaction, the reactive proton is exchanged with deuterium via ``isotope_indices 12`` and ``isotope_num 2``, which together specify the 12th atom, :math:`^{1}H`, to be exchanged with :math:`^{2}H`. The PI beads distribution of the deuterium atom will be generated by mass-perturbation, which allows the quantum partition functions for H- and D-isotope systems to be computed within a single simulation. In a PI-FEP simulation output, the QM-PMF as well as the ratio of quantum partition functions of the light and heavy isotope systems are displayed at the end. The log file, ``cc2-s50.log``, stores the essential output: .. code-block:: bash :linenos: :caption: cc2-s50.log 1 0.0218207855 +/- 0.0011852039 0.0215333804 +/- 0.0010824629 2 0.0214098885 +/- 0.0014015810 0.0200354439 +/- 0.0016564722 mean delta G (cm -> qm^L) = 0.0215931591 +/- 0.0002022663 Hartree mean delta G (cm -> qm^H) = 0.0205141115 +/- 0.0006233911 Hartree exp(-beta deltaG(H->L)) = 4.2956460285 +/- 0.0317993487 ``ts2.xyz`` contains 2 classical nuclear coordinates. The first two lines display the QM-PMF for light and heavy isotopes at each classical nuclear coordinate, respectively. The total QM-PMF values for light and heavy isotope systems are reported as ``mean delta G(cm->qm^L)`` and ``mean delta G(cm->qm^H)``, respectively. The ratio of quantum partition functions, calculated using the FEP expression, is displayed as ``exp(-beta deltaG(H->L))``. 2 standard errors are used to represent all uncertainties. The KIE can then be easily obtained via: .. math:: \begin{align} \mathrm{KIE} = \frac{k^\mathrm{L}}{k^{\mathrm{H}}} = \left[ \frac{ Q^{ \mathrm{L} }_{ \mathrm{QM} } (\bar{z}^{ \mathrm{TS} } )}{ Q^{ \mathrm{H} }_{ \mathrm{qm} } (\bar{z}^{\mathrm{TS}} ) } \right] / \left[ \frac{ Q^{ \mathrm{L} }_{ \mathrm{QM} } (\bar{z}^{ \mathrm{RS} } )}{ Q^{ \mathrm{H} }_{ \mathrm{qm} } (\bar{z}^{\mathrm{RS}} ) } \right] \end{align} where :math:`\bar{z}^{\mathrm{TS}}` and :math:`\bar{z}^{\mathrm{RS}}` denote the reaction coordinates at the transition and reactant states, respectively. Example 2: Proton/Deuterium transfer of carbon acid in aqueous using AM1/CHARMM model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ We can also use PI-FEP simulations in conjunction with QM/MM methods to study proton transfer reactions in aqueous solution. In this example, we provide a PI-FEP simulation of proton/deuterium transfer of carbon acid in water solvent. The quantum mechanically important carbon acid is described using the AM1 model, and the water solvent is described using the TIP3P model. The carbon acid in a TIP3P waterbox can be built using CHARMM-GUI with appropriate periodic boundary condition (PBC) setup, and equilibrated with a number of steps of MD simulations using the NVT and NPT ensembles. We assume the CM-PMF is already obtained from classical MD simulation with umbrella-sampling/metadynamics. The molecular definition of carbon acid in a 39×39×39 Å water cube is provided in :download:`geom.xyz `, CHARMM parameter (:download:`lig.prm `, :download:`par_all36_cgenff.prm `, :download:`toppar_water_ions.str `), and topology (:download:`step2_solvator.psf `) files are provided. The classical coordinate file, :download:`rs-2.xyz `, which stores 2 classical nuclear coordinates of the reaction at the reactant state, is provided for demonstration purposes. .. code-block:: bash :linenos: :caption: pifep-rs.inp nddo charge -1 spin2p1 1 end mol geom.xyz end charmm parameters lig.prm par_all36_cgenff.prm toppar_water_ions.str topology step2_solvator.psf # rs-opt scaling14 1.0 rcutoff 11.0 rswitch 10.0 use_PBC Lbox 39 39 39 end qmmm qm_region 1-23 end pifep sampling bisection bisection_level 3 classical_coord_file rs-2.xyz checkpoint_fn chk_cc2_p8_s10 #restart chk_cc2_p8_s10 num_pifep_samples 10 quantised_atoms_list 1-3 5 8 9 12-14 # 1-23 if left out in blank isotope_indices 12 isotope_num 2 write_log ie-cc2-p8-s10.log end task pi am1/charmm end The PI beads quantization and isotope exchange settings are identical to those in Example 1. The simulation output format and interpretation follow the same description as provided in Example 1. Checkpoint and restart ------------------------- A text-based checkpoint file (named ``chk`` by default) is written as the simulation runs, if the simulation terminates due to any unexpected reasons, it can be restarted by ``restart`` option with the appropriate checkpoint filename. Code Execution in Parallel ------------------------------- The PI-FEP module supports both serial and MPI parallel computation, and can be used in conjunction with OpenMP. In the MPI parallel approach, each PI bead is assigned to a dedicated MPI worker process, while OpenMP enables further parallelization within each process by distributing computational tasks across multiple threads. Note that the number of MPI processes should not exceed the number of PI beads + 1. Assuming MPI is properly installed and configured, one can execute the program using the following: .. code-block:: bash :linenos: # 8 beads PI-FEP simulation executed using 9 MPI processes, 1 master thread is used to distribute tasks mpirun -np 9 qbics-linux-cpu-mpi pifep-rs.inp # 8 beads PI-FEP simulations executed using 5 MPI processes, each MPI process utilises 4 OpenMP threads mpirun -np 5 qbics-linux-cpu-mpi pifep-rs.inp -n 4 # PI-FEP simulation executed using 24 threads via OpenMP qbics-linux-cpu pifep-rs.inp -n 24 For PI-FEP simulations of any kinds of systems, parallel execution via MPI is strongly recommended for optimal performance.