Path Integral Free-energy Perturbation ================================================== .. contents:: :local: Path integral free-energy perturbation (PI-FEP) is an efficient method for computing KIEs of chemical reactions. It has been implemented into **Qbics** and can be used in conjunction with electronic structure theory, density functional theory, **CHARMM** force field method, QM/MM, and some self-defined potential models. Here we intent to deliver a brief yet comprehensive tutorial with a few examples. For any theoretical detail related to PI-FEP simulation, here is the list of articles one should refer to: - Major, D. T.; Gao, J. `Implementation of the bisection sampling method in path integral simulations. `_ *J. Mol. Graphics Modell.* **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. Example 1: Proton transfer of :math:`\mathrm{CH}_4 + \mathrm{CH}_3` described using b3lyp-d3/6-31+G(d) ------------------------------------------------------------------------------------------------------------------------------- Input file +++++++++++ PI-FEP simulations are performed on the basis of classical reaction coordinates are known, which, in general, can be obtained using classical molecular dynamics(MD) simulation. In this reaction, we define :math:`r_1(\mathrm{C-H}) - r_2(\mathrm{H-C})` as the reaction coordinate. Classical MD simulations with **wham** were pre-performed at ``300K`` to obtain the classical potential of mean force. Classical nuclear coordinates describing reactant and transition states of the reaction were extracted from the bins and stored into ``rs3config.xyz`` and ``ts3config.xyz``, respectively, which serve as our input for further PI-FEP simulations. Note, generally hundreds to thousands of those classical nuclear coordinates are required in order to achieve a good numerical convergence, here we only choose 3 for each state for demonstration. In this example, we aim to compute the KIE for the following reactions: :math:`k_1 : \mathrm{CH}_4 + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_4` :math:`k_2 : \mathrm{CH}_3\mathrm{D} + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_3\mathrm{D}` Here is an input script: .. code-block:: bash :linenos: basis 6-31g(d) end mol C 2.40521526000000 11.31482711000000 -5.72742519000000 H 3.38961363000000 11.74694615000000 -5.95680724000000 H 2.53807336000000 10.39290699000000 -5.21807680000000 H 1.83040992000000 11.97651941000000 -5.09613657000000 H 1.97154650000000 11.15308162000000 -6.72243539000000 C 1.15156814000000 8.97757611000000 -3.51078100000000 H 0.70253665000000 8.42147626000000 -4.39345334000000 H 2.16310940000000 8.68858290000000 -3.12958401000000 H 0.50722084000000 9.84987092000000 -3.03707199000000 end basinfo angular_type spherical linear_dependence 1.E-8 end scf charge 0 # The net charge. spin2p1 2 # 2S+1 temp 0 # The electronic temperature. type U # Two options: R, U max_it 100 # The maximum number of iterations. energy_cov 1.E-6 # The energy convergence threshold. density_cov 1.E-6 # The density convergence threshold. schwarz 1.E-10 # The Schwarz screening threshold. end pifep temp 300. # in Kelvin action pa sampling bisection bisection_level 3 classical_coord_file rsconfig.xyz # for reactant state #classical_coord_file tsconfig.xyz # for transition state checkpoint_fn p8chk num_pifep_samples 5 isotope_indices 3 isotope_num 2 quantised_atoms_list 1 3 6 # quantise C-H-C region write_log kie_test.log random_seed 50000 end grimmedisp type bj end task pi b3lyp end The PES of two reactions are described using ``b3lyp-d3/6-31+G(d)``, the SCF and basis-set settings can be found under ``scf`` and ``basinfo`` keywords, respectively. For such small basis set, Grimme's ``D3`` dispersion correction is used and specified under ``grimmedisp`` keyword block. An arbitrary geometry of :math:`\mathrm{CH}_4 + \mathrm{CH}_3` is provided under ``mol``, which allows to construct molecular data, `i.e.` atomic masses, molecular orbitals. Under ``task`` keyword, ``pi`` is used to trigger a path integral simulation, all PI-FEP simulation related parameters are specified within ``pifep`` keyword block. Keyword ``temp`` is used to set simulation temperature in Kelvin, a positive real number is taken. High-order PI action formalisms are available given only electronic energy and its gradients are required due to its in conjunction with FEP theory. Keywords for specifying high-order action schemes and their subsequential additional parameters can be found in manual; Here we use low-order primitive action scheme, ``pa``, because no evaluations of gradients are required. As the name suggests, ``bisection`` is used for sampling PI beads in this example, consequently, the number of beads that defines a PI ring-polymer equals :math:`2^N`, where :math:`N` is a positive integer and it is termed bisection level, ``bisection_level``. In this reaction, the classical nuclear coordinates are provided in ``classical_coord_file``, the geometry must be provided in an identical way in a single ``.xyz`` file as it was defined in ``mol``. Checkpoint feature may be enabled with ``checkpoint_fn``, a file name may be given or ``chk`` is used by default. The thermodynamic quantum corrections to PMF are obtained using double averaging method, i.e. the inner average over the number of bisection samplings performed based on each classical nuclear coordinate, and the outer average over total number of classical nuclear coordinates. The sample size of inner average, that is the number of bisection samplings in this example, is specified using ``num_pifep_samples``. In order to perform PI-FEP simulation for the :math:`\mathrm{k}2` reaction, the reactive :math:`\mathrm{H}` is exchanged into :math:`\mathrm{D}`; here we use ``isotope_indices 3`` in conjunction with ``isotope_num 2`` to indicate that atom index 3 of ``mol``, :math:`\mathrm{H}`, is exchanged with its 2nd isotope, :math:`\mathrm{D}`. (Note, ``isotope_indices`` takes a number of input formats, refer to manual for detail.) For computational efficiency concerns, one often define a quantum mechanically important region for the system of interest, that is to `quantise` a group of atoms with PI beads, this is handled by ``quantised_atoms_list`` keyword. Here for demonstration we assume :math:`\mathrm{C}-\mathrm{H}-\mathrm{C}` and :math:`\mathrm{C}-\mathrm{D}-\mathrm{C}` are the regions which are quantum mechanically important, and specify atoms 1 3 6, respectively, to be quantised with PI beads. Leaving the argument blank results in all atoms to be quantised. The inner average, which is the mean correction to the potential energy of a given classical nuclear configuration for :math:`k1` as well as :math:`k2` reactions, are printed in neat in a text based file via ``write_log`` keyword, file named ``log.dat`` is used by default. The reproducibility of the work is guaranteed by initialising the pseudo-random number generator a constant integer, which is handled by `random_seed` keyword followed by a integer value. Without the keyword the program will take a time-based random integer instead. Execute the program ++++++++++++++++++++ PI-FEP implementation in **Qbics** supports **MPI** feature, two versions are available: 1. Evaluations of potential and its gradients w.r.t. each PI bead index, i.e. a nuclear configuration, is assigned to a MPI process. 2. A classical nuclear coordinate is assigned to one individual MPI process, by which all bisection samplings and subsequential mean potential evaluations for the ring-polymers are computed. .. note:: The following keyword .. code-block:: bash pifep mpi_run end is required to enable the second MPI implementation. The second MPI implementation is ~30% faster than the first. However, ``write_log`` is not supported currently in second implementation. When performing PI-FEP simulations, one needs to increase the number of PI beads until the expectation values converge within defined criteria. Here is an standard shell script for performing such calculation: .. code-block:: bash :linenos: #!/bin/bash # performing pifep calculations for b in {3..5..1} # bisection level do p=$((2**${b})) # no. beads cat > input.inp << EOF basis 6-31g(d) end mol seed.xyz end basinfo angular_type spherical linear_dependence 1.E-8 end scf charge 0 # The net charge. spin2p1 2 # 2S+1 temp 0 # The electronic temperature. type U # Two options: R, U max_it 100 # The maximum number of iterations. energy_cov 1.E-6 # The energy convergence threshold. density_cov 1.E-6 # The density convergence threshold. schwarz 1.E-10 # The Schwarz screening threshold. end pifep temp 300 action pa sampling bisection bisection_level ${b} classical_coord_file rs3config.xyz # RS, use ts3config.xyz for TS checkpoint_fn ${p}chk_s10 #restart ${p}chk_s10 num_pifep_samples 5 isotope_indices 3 isotope_num 2 quantised_atoms_list 1 3 6 mpi_run end grimmedisp type bj #zero end task pi b3lyp end EOF mpirun -np 4 qbics input.inp -n 4 > pifep_${p}beads.dat # Execution of PI-FEP using 4 parallel processes (1 master), b3lyp/6-31+G(d) calculation using 4 OMP threads # qbics input.inp -n 16 > pifep_${p}beads.dat # Execution of PI-FEP in single processing, b3lyp/6-31+G(d) calculations using 16 OMP threads done **MPI** and **OpenMP** can be used in conjunction to perform a PI-FEP simulation, see line 56 and 57. Another few sample scripts are available :download:`here ` Analyse PI-FEP results +++++++++++++++++++++++ Executing above script produces 3 output files, `i.e.` ``pifep_8beads.dat``, ``pifep_16beads.dat``, and ``pifep_32beads.dat``, for PI-FEP simulations using 8, 16 and 32 beads, respectively. Crucial results for calculations of KIEs are printed in the end. Here's an example .. code-block:: text :linenos: ============== Total mean obtained from double average =============== delta G (cm -> qm) = 0.01089060 +/- 0.00018964 Hartree delta G (cm -> qm) = 0.00752452 +/- 0.00014664 Hartree Thermodynamic quantum corrections, :math:`\Delta G(\mathrm{cm} \rightarrow \mathrm{qm})` to classical PMF are obtained using double averaging. Line 2 and 3 are thermodynamic quantum corrections to classical PMF as well as their :math:`2\sigma` for :math:`k_1` and :math:`k_2` reactions, and here they are denoted as :math:`\Delta G_{\mathrm{L}}(\mathrm{cm} \rightarrow \mathrm{qm})` and :math:`\Delta G_{\mathrm{H}}(\mathrm{cm} \rightarrow \mathrm{qm})`, respectively. For the sake of conveniences, here we drop the braces. The KIEs are obtained using the following expression: :math:`k_1 / k_2 = e^{-\beta \left[ (\Delta_{\mathrm{L}}^{\mathrm{TS}} - \Delta_{\mathrm{L}}^{\mathrm{RS}}) - ( \Delta_{\mathrm{H}}^{\mathrm{TS}} - \Delta_{\mathrm{H}}^{\mathrm{RS}} )\right]}` where superscripts :math:`\Delta G^{\mathrm{RS}}` and :math:`\Delta G^{\mathrm{TS}}` denotes the quantum corrections for reactant and transition states, respectively. .. note:: There are two source of errors: 1. Sample size for outer average: the number of classical nuclear coordinates. 2. Sample size for inner average: the number of bisection samplings performed based on each classical nuclear coordinates. Quantum corrections along the reaction coordinates can be obtained via performing PI-FEP simulations along the classical coordinates. Example 2: Proton transfer of :math:`\mathrm{CH}_4 + \mathrm{CH}_3` in (6,6)-carbon-nanotube described using QM/MM -------------------------------------------------------------------------------------------------------------------- In this example, we demonstrate how to perform PI-FEP simulation with QM/MM method, this example can also be referred to for PI-FEP simulation with MM force field method. The KIE of the following two reactions, :math:`k_1 : \mathrm{CH}_4 + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_4` :math:`k_2 : \mathrm{CH}_3\mathrm{D} + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_3\mathrm{D}` taking place in (6,6)-carbon-nanotube(CNT) are studied. The potential energy surface of the system is described using QM/MM treatment, where the reaction is treated using b3lyp-d3/6-31+G(d), and the interaction with the environment is treated using Charmm force field method. Input file +++++++++++++++ Given the ``66cnt_ch4ch3.inp`` file below: .. code-block:: bash :linenos: basis 6-31g(d) end basinfo angular_type spherical # Two options: spherical; Cartesian linear_dependence 1.E-8 end scf charge 0 # The net charge. spin2p1 2 # 2S+1 temp 0 # The electronic temperature. type U # Two options: R, U energy_cov 1.E-6 # The energy convergence threshold. max_it 100 # The maximum number of iterations. density_cov 1.E-6 # The density convergence threshold. schwarz 1.E-10 # The Schwarz screening threshold. end scfguess type atmden end grimmedisp type bj end charmm parameters par_all36_cgenff.prm topology nanotube_solute.psf scaling14 1.0 rcutoff 11.0 rswitch 12.0 end qmmm qm_region 2425-2433 # Define the QM region with atom indices. end mol temp.xyz end pifep temp 300 action pa sampling bisection bisection_level 3 classical_coord_file rs.xyz # RS, use ts.xyz for TS checkpoint_fn chk num_pifep_samples 5 isotope_indices 2429 # the reactive H atom isotope_num 2 # H->D isotope exchange quantised_atoms_list 2425 2429 2430 # C-H-C are quantised with PI beads mpi_run end task pi b3lyp/charmm end Similar to example 1, :math:`r_1(\mathrm{C-H}) - r_2(\mathrm{H-C})` is defined as the reaction coordinate, MD simulation was performed at 300K, with **wham** algorithm classical nuclear coordinates at reactant and transition states were extracted into ``rs.xyz`` and ``ts.xyz`` files, separately. An arbitrary geometry of :math:`\mathrm{CH}_4 + \mathrm{CH}_3` in (6,6)-CNT was written in ``temp.xyz`` for ``mol``. The SCF, basis-set, and dispersion correction settings are found under ``scf``, ``basinfo``, and ``grimmedisp`` keyword blocks, respectively. Charmm force field settings are specified within ``charmm`` keyword block, any further detail should be referred to the corresponding tutorial. .. note:: Charmm force field parameter file, ``par_all36_cgenff.prm``, and topology file, ``nanotube_solute.psf``, for describing (6,6)-CNT, as well as above mentioned ``.xyz`` files are available :download:`here `. PI-FEP simulation for :math:`k_2` reaction is achieved via isotope exchange to the reactive :math:`\mathrm{H}`, atom 2429, into its 2nd isotope, :math:`\mathrm{D}`. Here we assume :math:`\mathrm{C-H-C}` is the region which is quantum mechanically important, and therefore they are quantised with PI beads. ``pi b3lpy/charmm`` is used to trigger PI-FEP simulation with QM/MM method under ``task`` keyword block. MPI feature is also available for QM/MM method, one can refer to example 1 in regard to executing a PI-FEP simulation. Example 3: PI-FEP with a user-defined model PES ------------------------------------------------ **Qbics** provides socket for users to add model PES. A simple Morse potential is used for demonstration: Input file +++++++++++ .. code-block:: bash :linenos: custom_pes_type pes_type morse_pes end custom_pes_parameters de 136.3 alpha 2.2112 r0 0.9166 end mol H 0.000000 0.0000001 0.0000001 F 0.916600 0.0000001 0.0000001 end pifep temp 300 # temperature in Kelvin action pa # pa, tia, sca, ca sampling bisection bisection_level 3 seed_geo_from_file hf-coords.xyz checkpoint_fn 8chk_b num_pifep_samples 10000 quantised_atoms_list isotope_indices 1 isotope_num 2 write_log hf.log end task pi custom_pes end Under ``task`` keyword block, ``custom_pes`` is used to trigger a simulation with a user-defined model PES. A morse PES model is triggered via providing ``pes_type morse_pes`` under ``custom_pes_type`` keyword block. A morse PES has the expression of :math:`V(r) = De \left(1 - e^{-\alpha (r-r_0)} \right)^2`, its parameters, :math:`De`, :math:`\alpha`, and :math:`r_0` are specified under ``custom_pes_parameters`` keyword block. The parameter set in use here describes a :math:`\mathrm{HF}` system. The isotope effect for exchanging the first atom, :math:`\mathrm{H}`, with its second isotope, :math:`\mathrm{D}`, is considered in this example. All atoms are quantised with PI beads. .. note:: MPI and OMP are not supported by ``custom_pes`` tasks unless otherwise explicitly implemented. The program can be executed in the same way as above examples. .. hint:: Morse PES is analytically smooth, second order derivatives are available in **Qbics** program. It is computationally much more affordable than QM methods. It is recommended to explore numerical convergence efficiencies using high-order action schemes. A collection of sample scripts are available :download:`here `