Skip to content

2. Use of sobEDA.sh

2.1 Running sobEDA.sh

After preparing the following files in current folder, and properly setting sobEDA.sh, then you can input path of sobEDA.sh script to start sobEDA/sobEDAw analysis:

  • system.xyz: Standard .xyz format file containing coordinates of the entire system. See Section 2.3 for details.
  • fragment.txt: Definition of net charge, spin multiplicity, atom indices in each fragment. See Section 2.4 for details.
  • template.gjf: Template input file of Gaussian, determining calculation level, net charge, spin multiplicity of the entire system, and computational resources. See Section 2.5 for details.

After running of sobEDA.sh, Gaussian and Multiwfn will be automatically invoked. Details of calculation procedure will be shown on screen so that you can monitor the progress, and finally sobEDA and sobEDAw terms will be displayed.

2.2 Setting sobEDA.sh

Before running sobEDA.sh, you should edit the beginning part of this script to properly set some options:

  • gau: Command of running Gaussian program. Assume that Gaussian 16 has been properly installed, just set it to g16.
  • mwfn: Command of running Multiwfn program. Assume that Multiwfn has been properly installed, just set it to Multiwfn.
  • iCP: If set to 0, then monomer basis functions will be used during fragment calculations. If set to 1, then complex basis functions will be used during fragment calculations. The latter is more expensive, but highly recommended for sobEDAw analysis to yield a better result. When iCP=1, it is equivalent to employing Counterpoise correction to account for basis set superposition error (BSSE) problem.
  • sobEDAw: If set to 0, only sobEDA result will be finally given. If set to 1, both sobEDA and sobEDAw results will be finally shown. In the latter case, you should also define fitted parameters c, a, r involved in sobEDAw method by parm_c, parm_a and parm_r variables, respectively. In sobEDA.sh you can find that many c, a, r parameters for commonly used calculation levels have already been given as comments, so you can just uncomment the corresponding line to use them. Parameters for more levels can be found in Tables 1 and 2 of the original paper of sobEDA method.
  • irun: This is for debug use. Commonly it should be set to 1. If set to 0, then the script will only process existing files generated by previous calculation, no additional calculations by Gaussian will be actually conducted.
  • idetail: Commonly it should be set to 0. If set to 1, full energy variation details will be finally printed, and additional computational cost is needed.

2.3 system.xyz

This file contains geometry of the entire system in .xyz format used for sobEDA/sobEDAw analysis. The first line is the number of atoms, the second line is a comment and does not affect analysis, and the rest of this file contains element name and X, Y, Z coordinates in Å for each atom. For example:

7
optimized
 O                 -0.03836200    1.54613400    0.00000000
 H                  0.06356500    0.57566500    0.00000000
 H                  0.85354600    1.90026300    0.00000000
 N                 -0.03836200   -1.37627100    0.00000000
 H                 -1.03728700   -1.54652400    0.00000000
 H                  0.34780000   -1.83229000    0.81779500
 H                  0.34780000   -1.83229000   -0.81779500

Tip: Generating .xyz via Multiwfn

You can easily generate the .xyz file via Multiwfn. Multiwfn can load a wide variety of files and export the structure and wavefunction to rich formats. See Section 2.5 of Multiwfn manual on the formats that can be loaded. For example, Multiwfn can directly load the output file of an opt or opt freq task of Gaussian; after loading it, in the main menu if you input command xyz, you will enter the interface of exporting .xyz file. After inputting the file path to be exported (for example, /sob/system.xyz), the system.xyz file will be generated. You can also visualize the loaded structure in Multiwfn via main function 0 (noGUI version of Multiwfn does not support this).

2.4 fragment.txt

The first line of this file is the number of fragments. The rest defines net charge, spin multiplicity and atom indices in each fragment. The following example defines two fragments: the first fragment has net charge of 0 and spin multiplicity of 1, corresponding to atoms 1-3, 7, 10-12; the second fragment has net charge of 1 and spin multiplicity of 2, corresponding to atoms 4-6, 8, 9. Evidently, the union of all fragments must equal the entire system.

2
0 1
1-3,7,10-12
1 2
4-6,8,9

You can use negative spin multiplicity to request Multiwfn to flip spin of the fragment when constructing promolecule wavefunction from corresponding fragment wavefunctions. In the following example, fragment 1 and fragment 2 will be viewed as doublet radical with α and β unpaired electrons during this process, respectively:

2
0 2
1-4
0 -2
5-8

Note

Note that spin flip should be properly set, so that the sum of α (β) electrons of all fragments is equal to α (β) electrons of the entire system.

2.5 template.gjf

This is template input file of Gaussian, a typical example is shown below. Note that there are two blank lines at the end of the file.

#p b3lyp/6-311+G(2d,p) em=gd3bj ExtraLinks=L608 nosymm

title line

0 1
[geometry]

-5 5

With the template file shown above, sobEDA.sh will invoke Gaussian to perform calculation for fragments and the entire system at B3LYP-D3(BJ)/6-311+G(2d,p) level. The net charge and spin multiplicity for the entire system is specified in this file, namely 0 and 1, respectively. DO NOT change #P, ExtraLinks=L608 and nosymm keywords, all of them are needed in the analysis. The [geometry] part will be automatically replaced with the actual coordinates of fragments and entire system by sobEDA.sh when invoking Gaussian to calculate them.

The ExtraLinks=L608 keyword is used to ask Gaussian to print detailed components of energy by link 608 module, while the -5 at the last line tells link 608 that the current DFT functional is B3LYP. If you use other DFT functional, you should correspondingly change this index. See IOp(3/74) in Gaussian IOps reference for details, the index for commonly used functionals is given here: -73 (MN15), -55 (M06-2X), -54 (M06), -53 (M06L), -35 (TPSSh), -13 (PBE0), -5 (B3LYP), -6 (B3PW91), -3 (BHandHLYP), 402 (BLYP), 1009 (PBE), 2523 (TPSS).

Range-separated Functionals Not Supported

Notice that because Gaussian (at least as of G16 C.02) is unable to print energy components for range-separated functionals, functionals such as CAM-B3LYP, ωB97X, and HSE06 cannot be used in sobEDA/sobEDAw analysis.

The 5 at last line in above example requests link 608 module to use ultrafine integration grid, which is used by Gaussian 16 by default during SCF process.

Only DFT-D3(BJ) parameters of limited DFT functionals are built-in in Gaussian code. As of G16 C.02, to use functionals such as TPSSh-D3(BJ) and BHandHLYP-D3(BJ), you have to manually define DFT-D3(BJ) parameters. Since G16 C.01, the parameters should be defined using IOps; for example:

  • Adding IOp(3/174=1000000,3/175=1035400,3/177=279300,3/178=4961500) for BHandHLYP-D3(BJ)

  • Adding IOp(3/174=1000000,3/175=2238200,3/177=452900,3/178=4655000) for TPSSh-D3(BJ)

These parameters come from supplemental information of GMTKN55 benchmark paper Phys. Chem. Chem. Phys., 19, 32184 (2017). If you are using Gaussian older than G16 C.01, you need to use environment variables to define the DFT-D3(BJ) parameters, see http://sobereva.com/210 for details.

Template files for some popular DFT functionals have been provided in templates folder in this package, you can directly use them after modifying the file name.

Although sobEDA/sobEDAw is defined based on dispersion-corrected DFT, it is also possible to ignore the dispersion correction, namely just removing the em=gd3bj keyword. In this case, the dispersion correction term reported by sobEDA.sh will be exactly zero.

To define the number of CPU cores and memory used in the sobEDA calculation, you can either set -P- and -M- in Default.Rou file of Gaussian, or directly add %nproc and %mem in the template.gjf.

Solvation Effect

It is important to note that solvation effect is not directly taken into account in sobEDA method, and implicit solvation model (scrf keyword) should not be directly specified in the template file. Please search “solvation” in original paper of sobEDA on how to correctly consider the solvation effect. Briefly speaking, contribution of solvation effect to interaction energy should be manually calculated via Gaussian, it is never difficult.

2.6 Files produced by sobEDA.sh

After this script finishes running, you will find following files in current folder:

  • fragment[index].gjf, out, fch: They are input file, output file and .fch file of corresponding fragment.
  • promol.gjf, out, chk, fch: They are input file, output file, .chk and .fch file of calculation of promolecule state. Note that promol.fch is generated by combining all fragment .fch files via Multiwfn.
  • final.gjf, out, chk: They are input file, output file and .chk file of calculation of final state.

You can manually check input and output files to examine if the calculations were conducted in a correct way. Also, you can use Multiwfn to load .fch file (final.chk can be manually converted to .fch via formchk utility) to visualize orbitals via main function 0 or check real space functions such as spin density to verify whether reasonable wavefunctions have been obtained.