4. Examples¶
In this section some examples of using sobEDA.sh are given. All files for each example can be found from the file package of this tutorial, only some key points will be described below. Gaussian 16 C.02 and Multiwfn 3.8(dev) updated on 2023-Jun-24 are used in these examples.
4.1 sobEDAw analysis for H2O···NH3 dimer at B3LYP-D3(BJ)/6-311+G(2d,p) w.CB level¶
Structure of this system:
All files of this example have been provided in “HOH-NH3” folder, please check them to understand how to set up this sobEDAw analysis. As you can see from the sobEDA.sh in this example folder:
- iCP is set to 1 to use complex basis functions for fragment calculations.
- sobEDAw is set to 1 to perform sobEDAw analysis after sobEDA analysis.
- The line defining c, a, r parameters for the present level, namely parm_c=0.550;parm_a=0.037;parm_r=2.750, has been uncommented to make them take effect.
Make sure that Gaussian and Multiwfn have been properly installed, then copy sobEDA.sh, system.xyz, fragment.txt and template.gjf from the example folder to current folder, then run:
to add executable permission and execute it. The information will be printed on both screen and result.txt in current folder. The result is:
*************************
***** Final results *****
*************************
Total interaction energy: -7.18 kcal/mol
Physical components of interaction energy derived by sobEDA:
Electrostatic (E_els): -12.15 kcal/mol
Exchange (E_x): -6.40 kcal/mol
Pauli repulsion (E_rep): 20.28 kcal/mol
Exchange-repulsion (E_xrep = E_x + E_rep): 13.88 kcal/mol
Orbital (E_orb): -5.38 kcal/mol
DFT correlation (E_DFTc): -2.70 kcal/mol
Dispersion correction (E_dc): -0.83 kcal/mol
Coulomb correlation (E_c = E_DFTc + E_dc): -3.53 kcal/mol
Variant of sobEDA for weak interaction (sobEDAw):
Note: 96.26% DFT correlation is combined with dispersion correction to yield a SAPT-like dispersion term
Electrostatic (E_els): -12.15 kcal/mol
Exchange-repulsion (including scaled DFT correlation): 13.78 kcal/mol
Orbital (E_orb): -5.38 kcal/mol
Dispersion (E_disp): -3.43 kcal/molAs you can see, both sobEDA and sobEDAw results are shown. ΔEint is -7.18 kcal/mol, the sobEDAw terms ΔEels, ΔExrep, ΔEorb, and ΔEdisp are -12.15, 13.78, -5.38 and -3.43 kcal/mol, respectively.
The final.gjf/out/chk, promol.gjf/out/chk/fch, fragment1.gjf/out/fch and fragment2.gjf/out/fch in current folder are intermediate files generated by sobEDA.sh, commonly you can simply delete them, while you can also check them to try to identify possible reasons if you encountered any error during executing the script or you found the result is strange.
4.2 sobEDA analysis for OC-BH3 at B3LYP-D3(BJ)/6-311G* level¶
This section illustrates how to perform a very simple sobEDA analysis for chemical bond interaction, OC-BH3 is taken as an instance. Structure of this system:
All files of this example have been provided in “OC-BH3” folder, please check them to understand how to set up this analysis. Because in this example we do not use complex basis functions for fragment calculations, therefore iCP in sobEDA.sh in this case is 0. Because this example involves strong interaction, it is meaningless to conduct sobEDAw analysis, so sobEDAw in sobEDA.sh is set to 0.
Execute sobEDA.sh like last example, you will obtain result very soon:
Total interaction energy: -43.40 kcal/mol
Physical components of interaction energy derived by sobEDA:
Electrostatic (E_els): -69.50 kcal/mol
Exchange (E_x): -65.15 kcal/mol
Pauli repulsion (E_rep): 221.03 kcal/mol
Exchange-repulsion (E_xrep = E_x + E_rep): 155.88 kcal/mol
Orbital (E_orb): -117.36 kcal/mol
DFT correlation (E_DFTc): -10.89 kcal/mol
Dispersion correction (E_dc): -1.52 kcal/mol
Coulomb correlation (E_c = E_DFTc + E_dc): -12.41 kcal/molAs you can see, ΔEint is -43.40 kcal/mol, the sobEDA terms ΔEels, ΔEx, ΔErep, ΔEorb, ΔEDFTc and ΔEdc are -69.50, -65.15, 221.03, -117.36, -10.89 and -1.52 kcal/mol, respectively. For simplifying discussion, the ΔExrep combined from ΔEx and ΔErep as well as the ΔEc combined from ΔEDFTc and ΔEdc, are also shown (155.88 and -12.41 kcal/mol, respectively).
4.3 sobEDA analysis for Cu+-H2O at TPSSh-D3(BJ)/SDD&6-311G* level¶
This section illustrates sobEDA analysis for a charged system, and we will use mixed basis set and pseudopotential this time. Structure of the example system Cu+-H2O is shown below:
All files of this example have been provided in “Cu+H2O” folder, please check them to understand how to set up this analysis.
The fragment.txt of this example is shown below. It can be seen that fragment 2 (atom 4), namely the Cu(I) cation, is calculated with net charge of 1 and spin multiplicity of 1.
The template.gjf of this example is evidently more complicated than the previous examples:
#p TPSSh/genecp em=gd3bj IOp(3/174=1000000,3/175=2238200,3/177=452900,3/178=4655000) ExtraLinks=L608 nosymm
opted geometry
1 1
[geometry]
-O -H
6-311G*
****
-Cu
SDD
****
-Cu
SDD
-35 5
There are some key points of this template file:
- Because parameters of DFT-D3(BJ) are not built-in for TPSSh in Gaussian 16 C.02 (the version used in the present tutorial), therefore according to Section 2.5, the parameters are supplemented using IOp(3/174=1000000,3/175=2238200,3/177=452900,3/178=4655000).
- -35 in the last line indicates TPSSh is used, as mentioned in Section 2.5.
- Net charge of 1 and spin multiplicity of 1 defined in this file correspond to the setting used in the calculation of the entire system (while the values for fragment calculations are defined in fragment.txt).
- Customized basis sets are used in this example, the same basis set setting is applied to calculations of both entire system and every fragment. Therefore, - must be added before the element names, requesting that Gaussian will not stop if one or more of these elements do not appear in the present calculation.
- Pseudopotential “SDD” is used, due to the same reason as above, - is added prior to Cu.
Run sobEDA.sh, you will obtain the following sobEDA result:
Total interaction energy: -51.03 kcal/mol
Physical components of interaction energy derived by sobEDA:
Electrostatic (E_els): -59.76 kcal/mol
Exchange (E_x): -20.29 kcal/mol
Pauli repulsion (E_rep): 65.08 kcal/mol
Exchange-repulsion (E_xrep = E_x + E_rep): 44.79 kcal/mol
Orbital (E_orb): -28.75 kcal/mol
DFT correlation (E_DFTc): -6.52 kcal/mol
Dispersion correction (E_dc): -0.78 kcal/mol
Coulomb correlation (E_c = E_DFTc + E_dc): -7.31 kcal/mol4.4 sobEDA analysis for CH3-NH2 at B3LYP-D3(BJ)/def2-TZVP level¶
This section illustrates sobEDA analysis for open-shell interfragment interaction. Structure of the example system CH3-NH2 is shown below:
All files of this example have been provided in “CH3NH2” folder, please check them to understand how to set up this analysis. The fragment.txt in this example is given below, it can be seen that both CH3 and NH2 fragments are calculated as doublet radical, however, the CH3 has unpaired electron of α spin, while NH2 has that of β spin (indicated as negative spin multiplicity):
The result of running sobEDA.sh is:
Total interaction energy: -97.50 kcal/mol
Physical components of interaction energy derived by sobEDA:
Electrostatic (E_els): -171.02 kcal/mol
Exchange (E_x): -115.55 kcal/mol
Pauli repulsion (E_rep): 469.61 kcal/mol
Exchange-repulsion (E_xrep = E_x + E_rep): 354.06 kcal/mol
Orbital (E_orb): -260.78 kcal/mol
DFT correlation (E_DFTc): -18.38 kcal/mol
Dispersion correction (E_dc): -1.39 kcal/mol
Coulomb correlation (E_c = E_DFTc + E_dc): -19.76 kcal/molBy the way, after finishing running sobEDA.sh, the intermediate files fragment1.fch and fragment2.fch can be found in current folder, you can use Multiwfn to plot spin density to understand distribution of unpaired electrons of the fragments. For example, boot up Multiwfn, load fragment2.fch, and then input:
Then you will see the following spin density isosurface map, from which we can find unpaired electron distributes above and below the NH2 plane, which fully meets our expectations:
4.5 sobEDAw analysis between multiple fragments: Water tetramer at BHandHLYP-D3(BJ)/6-31+G(d,p) w.CP level¶
This section illustrates sobEDA analysis for interaction between more than two fragments. Structure of example system water tetramer is shown below:
All files of this example have been provided in “water_tetramer” folder, please check them to understand how to set up this calculation. The fragment.txt in this example is shown below, it can be seen that four fragments are defined, each one corresponds to a water molecule:
The template.gjf used in this example is shown below. It can be seen that DFT-D3(BJ) parameters for BHandHLYP are supplemented by IOp options because they are not built-in (at least as of G16 C.02):
#p BHandHLYP/6-31+G(d,p) em=gd3bj IOp(3/174=1000000,3/175=1035400,3/177=279300,3/178=4961500) ExtraLinks=L608 nosymm
opted geometry
0 1
[geometry]
-3 5
Note that in the sobEDA.sh in this example folder:
- iCP has been set to 1.
- sobEDAw has been set to 1.
- The line defining c, a, r parameters for sobEDAw corresponding to the BHandHLYP-D3(BJ)/6-31+G(d,p) w.CP level has been uncommented to make them take effect.
After running sobEDA.sh, you will obtain the following result:
Total interaction energy: -33.15 kcal/mol
Physical components of interaction energy derived by sobEDA:
Electrostatic (E_els): -57.56 kcal/mol
Exchange (E_x): -46.49 kcal/mol
Pauli repulsion (E_rep): 116.98 kcal/mol
Exchange-repulsion (E_xrep = E_x + E_rep): 70.48 kcal/mol
Orbital (E_orb): -28.90 kcal/mol
DFT correlation (E_DFTc): -14.62 kcal/mol
Dispersion correction (E_dc): -2.56 kcal/mol
Coulomb correlation (E_c = E_DFTc + E_dc): -17.17 kcal/mol
Variant of sobEDA for weak interaction (sobEDAw):
Note: 97.29% DFT correlation is combined with dispersion correction to yield a SAPT-like dispersion term
Electrostatic (E_els): -57.56 kcal/mol
Exchange-repulsion (including scaled DFT correlation): 70.09 kcal/mol
Orbital (E_orb): -28.90 kcal/mol
Dispersion (E_disp): -16.78 kcal/mol4.6 sobEDA analysis for a series interfragment distances of CH3-CN¶
This section illustrates sobEDA analysis for a series interfragment distances of CH3-CN. Equilibrium structure of this system is shown below:
All files of this example have been provided in “CH3-CN_scan” folder. From the folder you can find scan.sh, in which variable i is looped from -5 to 40 by 1 each time, corresponding to displacing CH3 fragment from -0.5 to 4.0 by 0.1 Å each time towards positive direction of Z axis. As can be seen from the script, the displacement of the fragment is realized by automatically invoking Multiwfn, whose option 1 of subfunction 7 of main function 300 is able to make displacement for specific atoms according to a given vector, and then the updated geometry is exported to system.xyz in current folder. For the system.xyz corresponding to each interfragment distance, scan.sh invokes sobEDA.sh to perform sobEDA analysis, and the sobEDA terms are extracted from the outputted information to scan.txt in current folder. Once running scan.sh is fully finished, the scan.txt contains displacement values and various sobEDA terms as different columns. Then you can use the data to plot various physical components of interaction energy with respect to displacements.
Content of the resulting scan.txt (the first line labelled meaning of each column):
Index E_tot E_els E_x E_rep E_xrep E_orb E_DFTc E_dc E_c
-5 137.62 -288.26 -256.09 1129.33 873.24 -414.49 -31.17 -1.70 -32.87
-4 1.60 -278.48 -205.90 867.43 661.53 -351.29 -28.48 -1.68 -30.16
-3 -74.50 -246.18 -164.36 665.15 500.79 -301.60 -25.86 -1.65 -27.51
-2 -114.27 -207.01 -130.35 508.85 378.50 -260.81 -23.34 -1.61 -24.95
-1 -132.06 -168.66 -102.73 388.11 285.39 -226.27 -20.95 -1.57 -22.51
0 -136.65 -134.51 -80.46 295.00 214.54 -196.47 -18.69 -1.52 -20.21
1 -133.44 -105.64 -62.63 223.39 160.76 -170.50 -16.60 -1.46 -18.06
2 -125.79 -82.02 -48.45 168.54 120.09 -147.78 -14.66 -1.41 -16.07
3 -115.80 -63.14 -37.26 126.72 89.47 -127.86 -12.90 -1.36 -14.25
4 -104.76 -48.31 -28.47 94.99 66.52 -110.35 -11.29 -1.33 -12.62
...ignoredNote
It is worth noting that the current scan task does not consider geometry variation of fragments during the change of interfragment distance. To account for this effect, you need to perform relaxed scan, and then perform sobEDA analysis for every resulting geometry.
4.7 sobEDAw analysis for a series interfragment distances of benzene···H2O¶
This section illustrates sobEDAw analysis for a series intermolecular distances of benzene···H2O. The optimized equilibrium structure at B3LYP-D3(BJ)/6-311G** level is shown below:
All files of this example have been provided in “benzene-water_scan” folder. From the folder you can find scan.sh, in which variable i is looped from -5 to 50 by 1 each time, corresponding to displacing water from -0.5 to 5.0 by 0.1 Å each time towards positive direction of Z axis. equil.xyz contains optimized dimer geometry corresponding to above map. As can be seen from template.gjf and sobEDA.sh, the current sobEDAw analysis will be conducted at B3LYP-D3(BJ)/6-31+G(d,p) w.CB level.
Now run scan.sh. Once the calculation is fully finished, the scan.txt in current folder will contain displacement values and various sobEDAw terms as different columns. At the same time, you can find scan.xyz in current folder, which contains geometry of each calculated structure, you can use VMD program to load it and visualize the structures as trajectory.
Content of the scan.txt:
Index E_tot E_els E_xrep E_orb E_disp
-5 3.83 -11.75 30.09 -5.78 -8.73
-4 1.00 -9.43 22.58 -4.47 -7.69
-3 -0.95 -7.63 16.93 -3.47 -6.78
-2 -2.24 -6.24 12.69 -2.72 -5.98
-1 -3.04 -5.16 9.53 -2.14 -5.27
0 -3.49 -4.31 7.17 -1.69 -4.65
1 -3.69 -3.65 5.41 -1.35 -4.10
2 -3.71 -3.12 4.11 -1.08 -3.62
3 -3.61 -2.70 3.14 -0.87 -3.18
...ignoredThe distance between center of benzene and oxygen of water in equil.xyz is 3.17637 Å, therefore the distance of the structure with index i is 3.17637 + 0.1 * i. Plot the various sobEDAw terms with respect to the distance as curve map:
4.8 Obtaining full energy variation terms of interfragment interaction¶
This example illustrates how to obtain full energy variation terms of interfragment interaction, as exhibited in Fig. 1 of original paper of sobEDA. OC-BH3 is taken as example.
All files of this example have been provided in “OC-BH3_full” folder. As you can see from the sobEDA.sh in this folder, idetail has been set to 1. After running this script, you will see the following content from the outputted information. It can be seen that variation of various energy components during the three steps of formation of final wavefunction from isolated fragment wavefunctions are shown clearly, they sum up to total interaction energy (-43.40 kcal/mol). Via the detailed information behind sobEDA, you can better understand how sobEDA terms are derived:
[...ignored]
Details of energy variations:
Total kin els exch corr disp
isolate->promol: -147.06 0.00 -69.50 -65.15 -10.89 -1.52 kcal/mol
promol ->frozen: 221.03 427.12 -193.69 -12.46 0.05 0.00 kcal/mol
frozen ->final: -117.36 -380.34 227.00 37.26 -1.28 0.00 kcal/mol
Whole process: -43.40 46.78 -36.19 -40.35 -12.12 -1.52 kcal/mol
*************************
***** Final results *****
*************************
Total interaction energy: -43.40 kcal/mol
[...ignored]4.9 More examples¶
Files for additional examples are also provided, please check them yourself:
C18_O2folder: sobEDAw for C18···O2 interaction at BHandHLYP-D3(BJ)/6-311+G(2d,p) w.CP level.(CO)5CrCH2folder: sobEDA for (CO)5Cr=CH2 interaction at TPSSh/6-311G*&SDD level.C6F6_Cl-folder: sobEDA for C6F6···Cl- interaction at BHandHLYP-D3(BJ)/6-311+G(2d,p) level.HC-CHfolder: sobEDA for interaction between two quartet CH fragments at B3LYP-D3(BJ)/def2-TZVP level. The files for analyzing interaction between two doublet CH fragments are given indoublet-fragmentsubfolder.