Tip

All input files can be downloaded: Files.

pifep

This option controls the path integral free-energy perturbation (PI-FEP) simulation.

Options 

Value	A real number
Default	`298.15`

Simulation temperature in Kelvin.

sampling

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.

bisection_level

Value	An integer
Default	`6`

The number of PI beads under bisection algorithm which is defined as \(2^N\), where \(N\) denotes bisection level. This keyword must be used in conjunction with bisection sampling algorithm.

classical_coord_file

Value	A filename
Default	None

Filename of an .xyz file which stores classical nuclear coordinates.

checkpoint_fn

Value	A filename
Default	`chk`

Checkpoint filename.

restart

Value	A filename
Default	None

The checkpoint filename from which the simulation resumes.

num_pifep_samples

Value	An integer
Default	`100`

The number of PI beads samplings to perform for each classical nuclear coordinate.

quantised_atoms_list

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.

isotope_indices

Value	List of integers or range of integers
Default	None

The indices of atoms which are specified for isotope exchange.

isotope_num

Value	List of integers
Default	None

The isotope number based on table below (the atomic masses are hardcoded for now)

Isotopes table
Isotope Name	isotope_num
\(^1\mathrm{H}\)	1
\(^2\mathrm{H}\)	2
\(^3\mathrm{H}\)	3
\(^{12}\mathrm{C}\)	1
\(^{13}\mathrm{C}\)	2
\(^{14}\mathrm{C}\)	3
\(^{14}\mathrm{N}\)	1
\(^{15}\mathrm{N}\)	2
\(^{16}\mathrm{O}\)	1
\(^{17}\mathrm{O}\)	2
\(^{18}\mathrm{O}\)	3

write_log

Value	A file name
Default	`log.dat`

Log filename which stores summarised simulation output.

random_seed (optional)

Value	An integers
Default	None

A seed value used to initialize the pseudo-random number generator (PRNG).

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 classical nuclear coordinates.

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:

\[\begin{split}\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}\end{split}\]

where \(\beta\) denotes reciprocal temperature times Boltzmann constant. \(\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(\(\mathrm{L}\)) and heavy(\(\mathrm{H}\)), at a given state, \(\bar{z}\), can be written as a free energy perturbation:

\[\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 \(\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, \(^{1}H\), to be exchanged with \(^{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:

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:

\[\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 \(\bar{z}^{\mathrm{TS}}\) and \(\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 geom.xyz, CHARMM parameter (lig.prm, par_all36_cgenff.prm, toppar_water_ions.str), and topology (step2_solvator.psf) files are provided. The classical coordinate file, rs-2.xyz, which stores 2 classical nuclear coordinates of the reaction at the reactant state, is provided for demonstration purposes.

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:

 # 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.

pifep

Options

Theoretical Background

Input examples

Example 1: Proton/Deuterium transfer of carbon acid using AM1 model

Example 2: Proton/Deuterium transfer of carbon acid in aqueous using AM1/CHARMM model

Checkpoint and restart

Code Execution in Parallel