5. Locate minima of IRI-π by basin analysis and topology analysis¶
In this section I illustrate how to locate minima of IRI-\(\pi\) by using basin analysis and topology analysis modules of Multiwfn. If you want to perform topology analysis for IRI rather than IRI-\(\pi\), please check example given in Section 4.2.11 of Multiwfn manual.
We will first perform basin analysis on IRI-\(\pi\), because in principle this function can find all extrema without omission.
Boot up Multiwfn and input:
thiophene-pi.mwfn
17 // Basin analysis
-1 // Select the method for generating basins
2 // After selecting this, minima rather than maxima (attractors) will be searched
1 // Generate basins and locate attractors
24 // Corresponding to IRI-π in this context
2 // Medium quality grid (Note: If you hope to reproduce Fig. 11 in IRI original paper, you should select “high-quality grid” here)
After calculation, we find the following information of attractors on screen (they correspond to IRI-\(\pi\) minima):
The attractors after clustering:
Index Average X,Y,Z coordinate (Angstrom) Value
1 -0.55563607 -0.97897784 -0.63965642 0.13482492
2 -1.82566138 -0.02645886 2.00622963 1.62806775
3 0.55563607 -0.97897784 -0.63965642 0.13482492
4 -0.29104747 -1.24356645 -0.00464377 0.44488094
5 1.82566138 -0.02645886 2.00622963 1.62806775
6 0.29104747 -1.24356645 -0.00464377 0.44488094
7 -0.55563607 0.02645886 -1.22175136 0.24197372
8 0.55563607 0.02645886 -1.22175136 0.24197372
9 -0.29104747 -0.71438923 -1.27466908 0.39114927
10 0.29104747 -0.71438923 -1.27466908 0.39114927
11 -0.50271835 -0.55563607 1.15954609 1.76399196
12 0.50271835 -0.55563607 1.15954609 1.76399196
13 -0.60855379 -0.02645886 0.41869800 1.78222197
14 0.60855379 -0.02645886 0.41869800 1.78222197
15 -0.55563607 0.97897784 -0.63965642 0.13482515
16 0.55563607 0.97897784 -0.63965642 0.13482515
17 -0.29104747 0.71438923 -1.27466908 0.39114929
18 0.29104747 0.71438923 -1.27466908 0.39114929
19 -0.29104747 1.24356645 -0.00464377 0.44488083
20 0.29104747 1.24356645 -0.00464377 0.44488083
21 -0.00000000 0.02645886 1.10662837 4.05937287
22 -0.50271835 0.55563607 1.15954609 1.76398209
23 0.50271835 0.55563607 1.15954609 1.76398209
24 -0.50271835 -0.02645886 1.68872330 0.60663899
25 0.50271835 -0.02645886 1.68872330 0.60663899By choosing option 0 we can visualize distribution and indices of located IRI-\(\pi\) minima:
The attractor indicated by the red arrow corresponds to a minimum of IRI-\(\pi\) on \(\text{C}_\alpha\text{-}\text{C}_\beta\) bond of the thiophene. From the attractor information we found its X, Y, Z coordinate is (0.55563607, -0.97897784, -0.63965642) Å, and the IRI-\(\pi\) value at this point is 0.13482492 a.u.
Notice that the position of the IRI-\(\pi\) minima obtained in the basin analysis module is not quite accurate, since each attractor is located at a grid point while the grid spacing is finite (0.1 Bohr in present case), and this is why the IRI-\(\pi\) value (0.1348 a.u.) at the located attractor is not quite close to 0. We can use topology analysis module to refine the position of IRI-\(\pi\) minima. For example, to refine attractor 3 in the above map, we input the following commands:
-10 // Return to main menu
2 // Topology analysis
-11 // Reselect the real space function whose critical points (CPs) will be searched
24 // IRI (corresponding to IRI-π in this context)
2 // Type of IRI is IRI-π
You will note from screen prompt that algorithm for locating CPs is automatically changed to Barzilai-Borwein, this is because it is found that the default Newton method is usually unable to successfully locate minima of IRI and IRI-\(\pi\) due to their unusual variation behaviors even if initial guessing position is good enough.
Then input:
1 // Search CPs from given starting points
1 // Manually input a starting point
0.55563607, -0.97897784, -0.63965642 A // The coordinate of the minimum of IRI-π in Å located by basin analysis module as mentioned above
Multiwfn reports that a CP has been successfully located:
Index Coordinate Type
1 1.02050664 -1.82601991 -1.21830162 (3,-1)Do not worry about it being of type (3,-1) rather than (3,+3), IRI-\(\pi\) is not a very well-behaved function. If you plot an isosurface map for IRI-\(\pi\), you can easily confirm that it is indeed a minimum.
Now we visualize it. Input 0 to return to upper level of menu and then choose option 0 to open GUI window, you will see (note that you can click "CP labels" check box to display CP indices):
The small sphere is the refined IRI-\(\pi\) (3,+3), namely minimum. In the text window, you can find its coordinate:
To check IRI-\(\pi\) value at this point, click "RETURN" to close the GUI window and then select option 7, then input index of this IRI-\(\pi\) CP, namely input 1, you will see:
Position (Bohr): 1.020506644029 -1.826019905376 -1.218301618380
Position (Angstrom): 0.540028859595 -0.966288120580 -0.644697452453
CP type: (3,-1)
Density of all electrons: 0.4915317130E-01
Density of Alpha electrons: 0.2457658565E-01
Density of Beta electrons: 0.2457658565E-01
Spin density of electrons: 0.0000000000E+00
Lagrangian kinetic energy G(r): 0.2208846083E-02
G(r) in X,Y,Z: 0.1727111900E-06 0.7891031504E-03 0.1419570221E-02
Hamiltonian kinetic energy K(r): 0.3052019403E-01
Potential energy density V(r): -0.3272904011E-01
Energy density E(r) or H(r): -0.3052019403E-01
Laplacian of electron density: -0.1132453918E+00
Electron localization function (ELF): 0.9864577594E+00
Localized orbital locator (LOL): 0.8955444670E+00
Local information entropy: 0.5405753062E-01
Interaction region indicator (IRI): 0.4039348777E-04
Reduced density gradient (RDG): 0.1318597754E-04
... [ignored]Since the IRI-\(\pi\) value at this point is extremely small (0.4039348777E-04 a.u.), it can be simply regarded as zero, indicating that this is the accurate position of IRI-\(\pi\) minimum. Also note that \(\pi\) electron density at this point is 0.049 a.u., which is in line with that labelled in Fig. 11 of original paper of IRI method.
By the way, in fact this point also corresponds to critical point of \(\pi\) electron density. This can also be easily numerically verified in Multiwfn by inputting following commands:
-10 // Return to main menu
1 // Output all properties at a point
1.020506644029, -1.826019905376, -1.218301618380 // The X, Y, Z coordinates of the IRI-π minimum searched by the topology analysis module
1 // The unit of inputted coordinate is Bohr
You can find the following information on screen:
Note: The following information is for electron density
Components of gradient in x/y/z are:
-0.8026242570E-09 0.8564212720E-06 0.1193509793E-05
Norm of gradient is: 0.1468987293E-05Indeed, gradient norm of electron density (corresponding to \(\pi\) electron density in the present case) at this point is essentially zero.
It is also possible to simultaneously refine a batch of IRI/IRI-\(\pi\) minima via topology analysis module. To this aim, after locating attractors using basin analysis module, you should choose option "-4 Export attractors as pdb/pqr/txt file", and then select corresponding option to export the located attractors as "attractors.txt". Then in the option "1 Search CPs from given starting points" of topology analysis module, select "4 Using starting points from a .txt file" and then input the path of the "attractors.txt".