Relative energy gradient (REG) analysis along a PES Scan (1D) or IRC

ADF Bonding Analysis QTAIM Interacting Quantum Atoms Relative Energy Gradient

Overview

This tutorial presents two examples of Relative Energy Gradient (REG) analysis: first, a one-dimensional PES Scan of the water dimer; second, an intrinsic reaction coordinate (IRC) calculation for an SN2 reaction between chloride and chloromethane.

In this tutorial, we present a Relative Energy Gradient (REG) analysis [1]-[2] based on a QTAIM/IQA energy decomposition. This REG analysis is capable of selecting energy terms such that they reflect the overall energy variations of the system. This is a crucial feature for gaining meaningful chemical insights from an energy partitioning method, such as IQA in this case. The efficacy of the method was originally demonstrated by the Popelier’s group in Manchester using the B3LYP-level description of the water dimer [1]. It was shown that a selected subset of energy terms can provide a semi-qualitative description of the system’s behavior. Importantly, the method is based on physically meaningful principles and does not require any prior knowledge of the system or arbitrary parameters. Overall, the REG analysis offers a reliable and intuitive metric for analyzing energy partitioning in chemical systems.

The Python/PLAMS script provided with this tutorial supports ADF results from both a one-dimensional PES Scan and an intrinsic reaction coordinate (IRC) calculation. It automatically detects the calculation type from the supplied .results directory and displays either PES Scan detected! or IRC detected!; no manual selection is required. Both calculation types require QTAIM/IQA results saved along the path.

Note

Using this script for published research requires citing the original REG paper [1]. Please note that another Python script (not compatible with ADF however) has been developed by Popelier’s group under a MIT license. It is freely available here

Contents of this tutorial:

• Example 1: PES Scan of the water dimer: set up the scan and QTAIM/IQA calculation, run the REG script, and interpret the results from the original paper.

• Example 2: IRC for chloride substitution in chloromethane: set up the IRC and QTAIM/IQA calculation, run the same script, and analyze the exchange energies of the breaking and forming C-Cl bonds.

Both examples use a Python/PLAMS script to extract IQA results along the reaction path and perform the REG analysis.

This tutorial assumes a basic familiarity with the GUI. If you are not familiar with using the GUI yet, you might take a look at the GUI introduction tutorials before continuing.

Example 1: PES Scan of the water dimer

The System

This tutorial studies a very simple system: two water monomers in their most stable dimer conformation as illustrated on this picture.

In the original paper, a PES Scan was first performed in 10 steps along the O…H hydrogen bond (from 1.25 Å to 3.25 Å) at the B3LYP/aug-cc-pVTZ level. You can find more details in Reference [1].

In this tutorial, we want to keep the same protocol. We will then choose the B3LYP functional and an equivalent STO basis set (TZ2P).

Importing the initial structure and setting the job parameters

You can download the initial structure from this XYZ file (or construct it directly on screen if you are familiar with the GUI building tools). In this file, we kept the same atom numbering with respect to Popelier’s original paper.

1. Download the XYZ file by clicking here

2. In AMSinput: File → Import Coordinates

3. Select the XYZ file you just downloaded using the file dialogue window

Now we can set the parameters for the PES Scan:

1. Select Task → PES Scan

2. On the right, click on

3. In the new panel, tick the Save results for all PES points check box

4. Select the two atoms of the hydrogen bond by keeping the shift key pressed

5. Click now on the sign next to H(3)…O(4) (path coordinate)

6. Check the initial value is set to 1.25 Å and set the final value to 3.25 Å

Your PES Scan settings should now correspond to this screen capture:

Now, follow the following steps for the QTAIM/IQA settings:

1. Switch to the QTAIM panel: Properties → QTAIM

2. In this new panel, tick the Perform QTAIM analysis check box

3. Select Full for the analysis level

4. Tick the Interacting Quantum Atoms (IQA) check box

5. Set Print to verbose

Finally, set the last general settings:

1. Switch back to the Main panel

2. Select XC Functional → Hybrid → B3LYP

3. Select Relativity → None

4. Select Basis set → TZ2P and Frozen Core → None

5. Set the Numerical Quality to Good for this tutorial

Note

Please note that we always recommend a high quality integration grid for IQA calculations (Good or even VeryGood).

Performing the PES Scan and the REG analysis

Save the settings and run the ADF calculation.

At the end of the calculation, we recommend to check the PES Scan in AMSMovie: SCM → Movie The total bond energy profile should be similar to the one we obtain on this picture:

You can now quit ADF (you do not need to save the changes): SCM → Quit All.

To perform the REG analysis, download our Python/PLAMS script ADFREG.py. This script has been tested with AMS2023 and should be compatible with later versions. You can now execute it using the following command:

$AMSBIN/plams ADFREG.py -v resultsdir=path/to/your/results/directory.results

For instance, if your AMS input was named dimer.ams, a dimer.results folder has been created. Consequently, we use the following command:

$AMSBIN/plams ADFREG.py -v resultsdir=dimer.results

After reading the results, the script displays the following analysis settings:

Use defaults values (recommended for beginners)?
* Original EDA: Intra: no decomposition; Inter: Vcl and VX
* Segmentation: extrema only
* R2 threshold: 0.9
* m(REG) limits: below -0.2 or above +0.2

Type yes (default) or no:

For this tutorial, we recommend typing yes and pressing Enter. Pressing Enter without typing anything also accepts these defaults:

• Original EDA: Intra: no decomposition; Inter: Vcl and VX: each atom is described by its total intra-atomic energy (Eintra), without splitting it into individual contributions. Inter-atomic interactions are split into Coulomb (VCl) and exchange (VX) terms, as in the original REG analysis.

• Segmentation: extrema only: the energy profile is divided into segments at its minima and maxima. For the water dimer in this tutorial, the minimum separates the compression and dissociation segments.

• R2 threshold: 0.9: only terms with a coefficient of determination R² greater than this threshold are included in the REG results tables and considered for plotting.

• m(REG) limits: below -0.2 or above +0.2: terms with slopes between these limits are omitted from the REG results tables and plotting menu.

If you answer no, the script asks you to configure these settings individually. You can choose a decomposition of the intra-atomic energy into kinetic (Vkin), Coulomb (VCl), and exchange (VX) contributions, include dispersion contributions (Vdisp) when a dispersion correction is active, segment the path at both extrema and inflection points, and change the R² threshold and m(REG) cutoff. Keep the recommended defaults to reproduce the analysis below.

The script reports dispersion when the input contains an active DISPERSION keyword. To include dispersion contributions (Vdisp), answer no and choose to include them when prompted.

Eventually, we obtain two output files: dimer_ADFREG.txt and dimer_ADFREG.csv. The first one contains the main REG results, and the second one exports the path coordinates, energies, and IQA descriptors for further processing.

The script then displays an interactive Plot options menu. The text and CSV files have already been written at this point:

Plot options:
1) Total energy vs path coordinate
2) Bond energy vs path coordinate
3) IQA term vs path coordinate
4) IQA term vs total energy
5) IQA term vs bond energy
6) Quit
Enter your choice (1-6):

• 1) Total energy vs path coordinate: plot the total energy against the path coordinate.

• 2) Bond energy vs path coordinate: plot the bond energy against the path coordinate.

• 3) IQA term vs path coordinate: select an IQA descriptor and plot its variation along the path.

• 4) IQA term vs total energy: plot an IQA descriptor against the total energy (a REG graph), including a linear regression for each selected segment. The legend gives the regression equation and R²; the slope is m_REG.

• 5) IQA term vs bond energy: produce the corresponding REG graph using the total bond energy instead of the total energy.

• 6) Quit: exit the script. Pressing Enter without entering a choice also quits.

For options 3, 4, and 5, a second menu asks you to select one segment or 0) All segments. The script then lists IQA terms whose R² exceeds the startup threshold and whose signed m(REG) is below the negative cutoff or above the positive cutoff on every selected segment. The defaults are R² > 0.9 and m(REG) < -0.2 or m(REG) > +0.2. Enter the number of the term you want to plot. A term present in the output files may therefore be absent from this menu. If no term meets both criteria, the script returns to Plot options.

For REG calculations, boundary points are excluded from both adjacent regression fits only when every affected segment still has at least four points after all its boundary exclusions. Otherwise, the boundary points are retained in the fits. All points remain visible in plots and included in the CSV output; the text results report how many points were used for each segment’s REG analysis.

The plots show the calculated points and smooth interpolating curves, with one color per segment. Close the plot window to return to the menu and choose another plot. For this tutorial, choose option 1 to display the segmented total energy profile shown below; option 3 lets you inspect individual IQA descriptors, and option 4 displays their REG graphs.

Script Troubleshoot

If the script does not run properly, please check the following points:

• Check the AMS version (2023 or later).

• Check the name and the path of your .results directory in the command line.

• Check if you did not forget to save the converged points of your PES Scan: Model → Geometry Constraints and PES Scan.

• In Details → Run Script, check if the runscript auto update check box is ticked.

Interpreting the results

1) Extracting chemical information from the PES Scan

The energetic profile is usually divided in several parts (segments), according to the extrema values of the total energy. In our example, our Morse potential exhibits a single minimum at 1.917 Å (hydrogen bond length). Therefore, we obtain two segments (see Fig. 5 in the original paper), the first one corresponding to a path coordinate varying from 1.250 Å to 1.917 Å while the total energy decreases and the second one to a path coordinate varying from 1.917 Å to 3.250 Å while the dimer dissociates towards a total energy corresponding to the sum of the two isolated fragment energies. The following graph, obtained with 1) Total energy vs path coordinate in the Plot options menu, illustrates this segmentation along the path coordinate.

The ADFREG script examines the variations of each IQA intra-atomic energy (Eintra) as well as IQA inter-atomic interaction energies (either VCl or VX for Coulomb and exchange parts respectively) along each segment. Then, a linear regression is performed for each of these IQA descriptors with respect to the variation of the total energy (or the total bond energy; only the intercept value will change). \(m_{REG,i}\) corresponds to the slope of each linear regression equation (one per IQA descriptor):

\[E_i(s) = m_{REG,i} \cdot E_{total}(s) + c_i\]

where \(E_{i}(s)\) corresponds to the values of each IQA descriptor i along the path coordinate (s) and \(E_{total}(s)\) corresponds to the values of the total energy along the path (or, again, the total bond energy since only the \(c_i\) will differ).

For each segment of the path, all the results of this REG analysis are saved in the first output textfile (job_name_ADFREG.txt).

A positive \(m_{REG}\) value (with a high value of R² or Pearson coefficient r) corresponds to a descriptor variation that correlates to the total energy variation and contributes positively to this variation. More precisely, if \(m_{REG}\) is greater than 1, the corresponding descriptor varies more rapidly than the total energy while a \(m_{REG}\) value inferior to 1 corresponds to a descriptor that varies more slowly than the total energy.

In this table (see job_name_ADFREG.txt for details), you find the positive \(m_{REG}\) values as well as statistical measures for the first segment (sorted in decreasing order).

IQA term	Atom(s)	R 2	r	m REG
Eintra	O4	0.975	0.988	1.926
VX	O1-H3	0.957	0.978	1.850
VCl	H3-H6	0.930	0.964	1.249
VCl	H3-H5	0.930	0.964	1.240
Eintra	H3	0.939	0.969	0.863
Eintra	H6	0.984	0.992	0.325
VX	O4-H6	0.966	0.983	0.286
VX	O4-H5	0.944	0.972	0.270

We clearly see that the first barrier (when shortening the O…H bond distance, compressing the dimer) is partly due to the increase of the Eintra(O4) energy term with the largest \(m_{REG}\) coefficient (1.926). This is in agreement with a paper [2] showing that steric repulsions correspond to an increase of intra-atomic energy terms. To a lesser extent, the H3 intra-atomic term is also correlated (R² = 0.939) to the total energy variation with a \(m_{REG}\) value close to 1 (0.863). The electrostatic repulsion VCl between O1 (HB donor) and O4 (HB acceptor) also increases as the distance decreases, in agreement with chemical intuition. Its R² value (0.864) is below the 0.9 cutoff, so it is not included in the table. Finally, we observe that H3 and H5 repel each other strongly (as well as H3 and H6) when the inter-fragment distance decreases too much. These terms also contribute to the barrier with a R² coefficient equal to 0.930 and a m_REG value greater than 1. Finally, the exchange (VX) between O1 and H3 decreases a lot at short distance (this O-H bond begins to break), contributing strongly (\(m_{REG}\) = 1.850) to the overall destabilization of the dimer.

The following table lists the negative \(m_{REG}\) coefficients with the corresponding statistical coefficients, from the most negative to the least negative value.

IQA term	Atom(s)	R 2	r	m REG
VCl	H3-O4	0.929	-0.964	-2.794
VX	H3-O4	0.952	-0.976	-1.408
VX	O1-O4	0.937	-0.968	-0.878
VCl	O4-H6	0.971	-0.985	-0.664
VCl	O4-H5	0.917	-0.957	-0.598
VX	O1-H2	0.967	-0.983	-0.241
VCl	H2-H3	0.954	-0.977	-0.209

The first terms of this list correspond to the IQA terms that work against the variation of the total energy. The main term (\(m_{REG}\) = -2.794 with R² = 0.929) corresponds to the coulombic interaction between H3 and O4 (hydrogen bond). Indeed, decreasing the distance between these atoms increases dramatically the attractive interaction between them. This stabilizes the dimer and counters the global destabilization. To a lesser extent, this hydrogen bond becomes stronger in terms of exchange (covalency). This term varies faster and against (\(m_{REG}\) = -1.408) the total energy, tempting to stabilize the dimer. Other terms may also be of interest but some others are meaningless due to poor correlations with the total energy variations (R² < 0.9; VCl(O1/H3) for instance) or slow variations (\(m_{REG}\) close to zero).

A similar analysis can be carried out on the second segment. You can find all the corresponding \(m_{REG}\) values as well as R² and Pearson coefficients in the same output file. Besides, a complete analysis of these segments can be found in the original work of Thacker and Popelier.

2) Plotting the variations of IQA descriptors

The script uses Matplotlib to generate the graphs shown in this tutorial. Use option 3 in the Plot options menu to plot an IQA descriptor along the path, or options 4 and 5 to produce REG graphs.

The CSV output file (job_name_ADFREG.csv) also makes the data available for further analysis, including terms not offered by the interactive plotting menu. It contains the following comma-separated values:

• the values of the path coordinates (all segments) in the first column.

• the values of the total bond energies along the path.

• the values of the total energies along the path.

• the values of each intra-atomic energy along the path (one per column).

• the values of each VCl inter-atomic term along the path (one pair per column).

• the values of each VX inter-atomic term along the path (one pair per column).

For example, choose 3) IQA term vs path coordinate, select 0) All segments, and then select VCl H3-O4 to plot the inter-atomic Coulomb energy of the hydrogen bond. The script produces the following graph:

Next, choose 4) IQA term vs total energy, select 0) All segments, and select the same VCl H3-O4 term. This produces a REG GRAPH, showing the correlations between the IQA descriptor and the total energy for each segment of the path.

The following graph corresponds to Fig.4 (bottom graph) in the original paper: VCl(H3-O4) against the total energy of the dimer. The script displays the linear regression equations and R² values for each segment in the legend.

Example 2: IRC for chloride substitution in chloromethane

We now apply the same REG analysis to the symmetric SN2 reaction:

\[\mathrm{Cl^- + CH_3Cl \longrightarrow ClCH_3 + Cl^-}\]

Although the reactants and products are chemically equivalent, the two chlorine atoms exchange their roles: the incoming chloride forms a C-Cl bond while the other C-Cl bond breaks. An IRC calculation follows the reaction path downhill from a transition state in both directions. Unlike the water dimer scan, its coordinate is the IRC path length rather than a selected inter-atomic distance.

Preparing the structure and the ADF calculation

Start from the transition-state geometry supplied below, not from the separated reactants. The atom numbering is Cl1, C2, Cl3, H4, H5, H6; keep this order to identify the IQA terms later. The two C-Cl distances in this starting structure are approximately 2.36 Angstrom, and the three hydrogens lie approximately in a plane through the carbon atom.

1. Download the SN2 transition-state XYZ file.

2. In a new AMSinput window, select File → Import Coordinates and import this file.

3. In the Main panel, select the ADF engine and Task → IRC.

4. Set the total Charge to -1 and use a closed-shell calculation.

5. Select XC Functional → Hybrid → B3LYP.

6. Select Basis set → TZP, Frozen Core → None, and Relativity → None.

7. Set Numerical Quality to Good.

Keep the default IRC settings: AMS calculates the initial Hessian and follows both directions from the transition state. The converged results must be retained for the subsequent IQA analysis. In the IRC settings, ensure that KeepConvergedResults is Yes; the corresponding AMS input block is:

IRC
    Direction Both
    KeepConvergedResults Yes
End

See the IRC documentation for the other settings. For a different system, first locate a transition state and verify that its imaginary mode corresponds to the intended reaction before starting the IRC.

Enable the same QTAIM/IQA analysis as for the water dimer:

1. Open Properties → QTAIM.

2. Tick Perform QTAIM analysis and select Full for the analysis level.

3. Tick Interacting Quantum Atoms (IQA) and set Print to verbose.

4. Save the job as Cl_CH3Cl_IRC.ams and run it.

When the calculation has finished, open SCM → Movie to inspect both branches. Check that one C-Cl distance decreases while the other increases, and that the energy decreases away from the transition state in each direction. The resulting Cl_CH3Cl_IRC.results directory should contain ams.rkf, TS.rkf, and the converged IRC_Backward_* and IRC_Forward_* result files. Keep the entire directory: the REG script reads the IQA results at each converged geometry, not only at the transition state.

Running the REG analysis for the IRC

Use the same ADFREG.py script as for the PES Scan. From the directory containing the script and the calculation results, run:

$AMSBIN/plams ADFREG.py -v resultsdir=Cl_CH3Cl_IRC.results

You can also supply the full path to the results directory. The script displays IRC detected! automatically. Answer yes to accept the recommended analysis defaults, including segmentation at extrema only. The output files are Cl_CH3Cl_IRC_ADFREG.txt and Cl_CH3Cl_IRC_ADFREG.csv, written next to the results directory.

In Plot options, choose 1) Total energy vs path coordinate. In the reference calculation, the script identifies two segments containing 27 points each, sharing the transition state at coordinate zero: segment 1 runs from approximately -0.882 to 0, and segment 2 from 0 to 0.884. These are the negative and positive branches of the IRC path length, not C-Cl bond distances. The total energy increases toward the transition state on segment 1 and decreases toward the product complex on segment 2.

The following plot shows the approximately symmetric energy barrier, with the transition state at the maximum. The blue curve corresponds to segment 1 and the orange curve to segment 2; the markers represent the calculated geometries.

Interpreting the C-Cl exchange terms

First, open Cl_CH3Cl_IRC_ADFREG.txt and examine the positive and negative REG rankings for each segment, together with their R² values. In segment 1, VX C2-Cl3 has the largest positive REG coefficient (10.549), while VX Cl1-C2 has the largest negative coefficient in absolute value (-5.046). In segment 2, their roles are exchanged. Both terms have high R² values, making them relevant descriptors for interpreting the energy variations. The corresponding results extracted from the text file are:

Segment	IQA term	R 2	m REG
1	VX C2-Cl3	0.994	10.549
1	VX Cl1-C2	0.961	-5.046
2	VX Cl1-C2	0.994	10.499
2	VX C2-Cl3	0.962	-5.014

These inter-atomic exchange energies describe the exchange contribution to the C-Cl interactions: a more negative VX corresponds to a more stabilizing exchange interaction.

On segment 1, the large positive m_REG for VX(C2-Cl3) identifies the loss of exchange stabilization in the original bond as a major contribution to the energy rise toward the transition state. The negative coefficient for VX(Cl1-C2) shows that exchange stabilization in the forming bond opposes that rise.

On segment 2, the total energy decreases. VX(Cl1-C2) now has a large positive coefficient: its increasingly stabilizing exchange interaction contributes to the energy decrease as the new bond forms. Conversely, VX(C2-Cl3) has a negative coefficient because the loss of exchange stabilization in the departing bond opposes that decrease. Thus, a positive REG coefficient does not always mean destabilization: its interpretation depends on the direction of the total energy change.

The nearly matching coefficients after exchanging Cl1 and Cl3 reflect the symmetry of this identity reaction. Small numerical differences are expected. The two C-Cl exchange terms contribute in opposite directions on each segment, but only partially offset each other. The remaining intra-atomic and Coulomb terms also contribute to the total energy variation.

Note

The plotting menu also helps identify potentially interesting terms by listing those with R² > 0.9 and \(|m_{REG}| > 0.5\) on every selected segment, together with their statistics. These thresholds are arbitrary screening criteria used by the script to shorten the list of terms offered for plotting; they are not requirements of the REG method and do not by themselves establish chemical relevance. Both C-Cl exchange terms satisfy these criteria on both segments.

Having identified VX C2-Cl3 and VX Cl1-C2 from the REG results, we can now plot them to inspect how the breaking and forming bonds evolve along the IRC:

1. Choose 3) IQA term vs path coordinate in Plot options.

2. Select 0) All segments.

3. Select VX C2-Cl3 from the numbered term list.

4. Close the plot window and repeat for VX Cl1-C2.

On the negative-coordinate side, C2-Cl3 is the initially bonded pair and Cl1 is the incoming chloride. Toward the transition state, VX(C2-Cl3) becomes less negative as this bond weakens, while VX(Cl1-C2) becomes more negative as the incoming chloride forms a bond. Beyond the transition state, C2-Cl1 becomes the bonded pair and Cl3 departs. The two curves therefore illustrate the exchange of bonding partners.

As an example, the following plot shows VX(C2-Cl3) increasing from approximately -0.233 a.u. in the reactant complex to -0.092 a.u. at the transition state and -0.013 a.u. in the product complex. This progressively less negative exchange energy illustrates the loss of exchange stabilization as the original C2-Cl3 bond breaks. The corresponding VX(Cl1-C2) curve follows the opposite trend as the new bond forms.

To relate these changes to the total energy, choose 4) IQA term vs total energy, select 0) All segments, and plot each of the same two terms. In these REG graphs, the horizontal axis is the total energy, not the reaction coordinate. The dashed lines are the linear regressions; their equations and R² values appear below each graph.

For the breaking C2-Cl3 bond, the blue branch has a positive slope (10.549), while the orange branch has a negative slope (-5.014). VX(C2-Cl3) increases along both segments, but the total energy first increases and then decreases, explaining the change in the sign of the REG coefficient.

For the forming Cl1-C2 bond, the graph is approximately symmetric to the preceding one, with the segment colors exchanged. The blue branch has a negative slope (-5.046) and the orange branch a positive slope (10.499): exchange stabilization in this bond opposes the energy rise before the transition state and contributes to the energy decrease after it.

References

[1] (1,2,3,4)

Thacker, J.C.R., Popelier, P.L.A., The ANANKE relative energy gradient (REG) method to automate IQA analysis over configurational change, Theoretical Chemistry Accounts, 136, 86 (2017)

[2] (1,2)

Symons, B., Williamson, D.J., Brooks, C.M., Wilson, A., & Popelier, P.L.A., Does the Intra-atomic Deformation Energy of IQA Represent Steric Energy?, ChemistryOpen 8, 560 (2019).