Tip

All input files can be downloaded: Files.

pifep

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

Options

temp

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

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