Path Integral Free-energy Perturbation

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:

Example 1: Proton transfer of \(\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 \(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:

\(k_1 : \mathrm{CH}_4 + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_4\)

\(k_2 : \mathrm{CH}_3\mathrm{D} + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_3\mathrm{D}\)

Here is an input script:

 1 basis
 2     6-31g(d)
 3 end
 4
 5 mol
 6      C      2.40521526000000     11.31482711000000     -5.72742519000000
 7      H      3.38961363000000     11.74694615000000     -5.95680724000000
 8      H      2.53807336000000     10.39290699000000     -5.21807680000000
 9      H      1.83040992000000     11.97651941000000     -5.09613657000000
10      H      1.97154650000000     11.15308162000000     -6.72243539000000
11      C      1.15156814000000      8.97757611000000     -3.51078100000000
12      H      0.70253665000000      8.42147626000000     -4.39345334000000
13      H      2.16310940000000      8.68858290000000     -3.12958401000000
14      H      0.50722084000000      9.84987092000000     -3.03707199000000
15 end
16
17 basinfo
18         angular_type       spherical
19         linear_dependence  1.E-8
20 end
21
22 scf
23     charge      0      # The net charge.
24     spin2p1     2      # 2S+1
25     temp        0      # The electronic temperature.
26     type        U      # Two options: R, U
27     max_it      100    # The maximum number of iterations.
28     energy_cov  1.E-6  # The energy convergence threshold.
29     density_cov 1.E-6  # The density convergence threshold.
30     schwarz     1.E-10 # The Schwarz screening threshold.
31 end
32
33 pifep
34     temp                    300.          # in Kelvin
35     action                  pa
36     sampling                bisection
37     bisection_level         3
38     classical_coord_file    rsconfig.xyz # for reactant state
39     #classical_coord_file    tsconfig.xyz # for transition state
40     checkpoint_fn           p8chk
41     num_pifep_samples       5
42     isotope_indices         3
43     isotope_num             2
44     quantised_atoms_list    1 3 6         # quantise C-H-C region
45     write_log               kie_test.log
46     random_seed             50000
47 end
48
49 grimmedisp
50         type bj
51 end
52
53 task
54         pi b3lyp
55 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 \(\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 \(2^N\), where \(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 \(\mathrm{k}2\) reaction, the reactive \(\mathrm{H}\) is exchanged into \(\mathrm{D}\); here we use isotope_indices 3 in conjunction with isotope_num  2 to indicate that atom index 3 of mol, \(\mathrm{H}\), is exchanged with its 2nd isotope, \(\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 \(\mathrm{C}-\mathrm{H}-\mathrm{C}\) and \(\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 \(k1\) as well as \(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

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:

 1#!/bin/bash
 2
 3# performing pifep calculations
 4for b in {3..5..1} # bisection level
 5do
 6    p=$((2**${b}))  # no. beads
 7cat > input.inp << EOF
 8basis
 9        6-31g(d)
10end
11
12mol
13       seed.xyz
14end
15
16basinfo
17        angular_type       spherical
18        linear_dependence  1.E-8
19end
20
21scf
22    charge      0      # The net charge.
23    spin2p1     2      # 2S+1
24    temp        0      # The electronic temperature.
25    type        U      # Two options: R, U
26    max_it      100    # The maximum number of iterations.
27    energy_cov  1.E-6  # The energy convergence threshold.
28    density_cov 1.E-6  # The density convergence threshold.
29    schwarz     1.E-10 # The Schwarz screening threshold.
30end
31
32pifep
33        temp                    300
34        action                  pa
35        sampling                bisection
36        bisection_level         ${b}
37        classical_coord_file    rs3config.xyz # RS, use ts3config.xyz for TS
38        checkpoint_fn           ${p}chk_s10
39        #restart                ${p}chk_s10
40        num_pifep_samples       5
41        isotope_indices         3
42        isotope_num             2
43        quantised_atoms_list    1 3 6
44        mpi_run
45end
46
47grimmedisp
48        type bj #zero
49end
50
51task
52        pi b3lyp
53end
54EOF
55    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
56    # 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
57done

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

1============== Total mean obtained from double average ===============
2delta G (cm -> qm) =      0.01089060 +/-      0.00018964 Hartree
3delta G (cm -> qm) =      0.00752452 +/-      0.00014664 Hartree

Thermodynamic quantum corrections, \(\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 \(2\sigma\) for \(k_1\) and \(k_2\) reactions, and here they are denoted as \(\Delta G_{\mathrm{L}}(\mathrm{cm} \rightarrow \mathrm{qm})\) and \(\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:

\(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 \(\Delta G^{\mathrm{RS}}\) and \(\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 \(\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,

\(k_1 : \mathrm{CH}_4 + \mathrm{CH}_3 \rightarrow \mathrm{CH}_3 + \mathrm{CH}_4\)

\(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:

 1basis
 2    6-31g(d)
 3end
 4
 5basinfo
 6    angular_type       spherical    # Two options: spherical; Cartesian
 7    linear_dependence  1.E-8
 8end
 9
10scf
11    charge      0      # The net charge.
12    spin2p1     2      # 2S+1
13    temp        0      # The electronic temperature.
14    type        U      # Two options: R, U
15    energy_cov  1.E-6  # The energy convergence threshold.
16    max_it      100    # The maximum number of iterations.
17    density_cov 1.E-6  # The density convergence threshold.
18    schwarz     1.E-10 # The Schwarz screening threshold.
19end
20
21scfguess
22    type       atmden
23end
24
25grimmedisp
26    type bj
27end
28
29charmm
30    parameters   par_all36_cgenff.prm
31    topology     nanotube_solute.psf
32    scaling14    1.0
33    rcutoff      11.0
34    rswitch      12.0
35end
36
37qmmm
38    qm_region  2425-2433 # Define the QM region with atom indices.
39end
40
41mol
42    temp.xyz
43end
44
45pifep
46    temp                    300
47    action                  pa
48    sampling                bisection
49    bisection_level         3
50    classical_coord_file    rs.xyz         # RS, use ts.xyz for TS
51    checkpoint_fn           chk
52    num_pifep_samples       5
53    isotope_indices         2429           # the reactive H atom
54    isotope_num             2              # H->D isotope exchange
55    quantised_atoms_list    2425 2429 2430 # C-H-C are quantised with PI beads
56    mpi_run
57end
58
59task
60    pi b3lyp/charmm
61end

Similar to example 1, \(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 \(\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 here.

PI-FEP simulation for \(k_2\) reaction is achieved via isotope exchange to the reactive \(\mathrm{H}\), atom 2429, into its 2nd isotope, \(\mathrm{D}\).

Here we assume \(\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

 1custom_pes_type
 2    pes_type    morse_pes
 3end
 4
 5custom_pes_parameters
 6    de          136.3
 7    alpha       2.2112
 8    r0          0.9166
 9end
10
11mol
12    H   0.000000    0.0000001   0.0000001
13    F   0.916600    0.0000001   0.0000001
14end
15
16pifep
17    temp                    300         # temperature in Kelvin
18    action                  pa          # pa, tia, sca, ca
19    sampling                bisection
20    bisection_level         3
21    seed_geo_from_file      hf-coords.xyz
22    checkpoint_fn           8chk_b
23    num_pifep_samples       10000
24    quantised_atoms_list
25    isotope_indices         1
26    isotope_num             2
27    write_log               hf.log
28end
29
30task
31        pi custom_pes
32end

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

\(V(r) = De \left(1 - e^{-\alpha (r-r_0)} \right)^2\),

its parameters, \(De\), \(\alpha\), and \(r_0\) are specified under custom_pes_parameters keyword block. The parameter set in use here describes a \(\mathrm{HF}\) system.

The isotope effect for exchanging the first atom, \(\mathrm{H}\), with its second isotope, \(\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 here