# 6.4. geom with Gaussian

Tip

Gaussian is suitable for **any clusters in gas phase**, as long as you choose the **correct model** (pure, hybrid, or double-hybrid DFT, CCSD(T), and basis sets), but it is also expensive. So, this should only be used when **there is no better ways**.

In this Section, we just want to demonstrate how to use Gaussian with `geom`

. Actually, it is possible that xTB is used first to get reliable structures, then Gaussian is used to refine, as discussed in Theoretical Background.

## 6.4.1. Example: (H2O)4(-)

Tip

The sample input and output files can be found in `testfiles/geom/2-h2o4e-g16`

.

In this Section we will see how to use `geom`

and Gaussian to do global optimization. Actually, this is exactly **the same** as we do in isomer with Gaussian: you just need to copy the commands to `commands`

block.

The cluster we consider here is \(\left(\text{H}_2\text{O}\right)_4^{-}\), an electron solvated in 4 water molecules. We will use the following input:

```
1lm_dir h2o4e # Save the local minima to this folder.
2num_calcs 50 # Total number of calculations.
3do_coarse_opt yes # no: Do NOT the coarse optimization.
4min_energy_gap 1.E-4 # When two energies differ smaller than
5 # this value, they are treated as identical.
6 # A negative number means do not remove
7 # energetically degenerated ones.
8max_geom_iters 3000 # The maximum number of iterations for local optimization.
9 # If it is less or equal than zero, then the number is unlimited.
10
11components
12 h2o.xyz 4
13 random 0 0 0 6 6 6
14 ****
15end
16
17savegjf
18 %chk=h2o4-$index$.chk
19 %nproc=48
20 %mem=30GB
21 $hash$MP2/aug-cc-pVTZ Opt Freq Output=wfn
22 $blank$
23 Initial guess generated from ABCluster, energy = $energy$
24 $blank$
25 -1 2
26 >>>>
27 h2o6-$index$.wfn
28end
29
30commands
31 xyz2gaussian optfile $inp$ > $xxx$.gjf
32 g16 < $xxx$.gjf > $xxx$.log 2>/dev/null
33 gaussian2xyz $xxx$.log > $out$
34end
```

You can see that, the commands to call Gaussian in `commands`

is the same to the ones in isomer with Gaussian. We will calculate `50`

structures and save them to `heo4e`

. We also add a `savegjf`

block to let `geom`

save LMs in Gaussian job input format.

Also, we use `random`

instead of `box`

because **we hope to get more flexible initial guess for global optimization.**

The structure of \(\text{H}_2\text{O}\) is:

```
13
2water
3O 0.00000000 0.00000000 -0.11081188
4H 0.00000000 -0.78397589 0.44324751
5H -0.00000000 0.78397589 0.44324751
```

The Gaussian input template is:

```
1%nprocs=48
2%mem=20GB
3%chk=h2o4e.chk
4#N B3LYP/6-31++g(d) SCF(XQC) Geom(NoCrowd) Opt
5
6opt
7
8-1 2
9>>>>
10>>>>
```

Now you can run the global optimization:

```
$ geom h2o4e.inp > h2o4e.out
```

After the optimization, the end of `h2o4e.out`

is

```
-- Result Report --
Results are energy-increasingly reordered.
Structures of energies within 1.000E-04 are treated as degenerate.
All minima are saved to "h2o4e".
-------------------------------------------------------------------
# index Energy NaiveRMSD
-------------------------------------------------------------------
0 32 -305.73627036 0.00000000
1 8 -305.73593799 1.92341824
2 1 -305.73517045 1.95589313
3 47 -305.73433736 2.14224379
4 16 -305.73402914 1.52860669
5 33 -305.73371713 1.31278433
6 34 -305.73312443 0.94350149
7 25 -305.73243751 1.49950503
8 23 -305.73225990 1.78428545
9 38 -305.73083814 2.34486039
10 29 -305.73055667 1.64935620
```

So the global minimum is `h2o4e/32.xyz`

, which is shown below. Interestingly, the anionic cluster optimized starting with the GM of neutral \(\left(\text{H}_2\text{O}\right)_4\) is higher in energy than `h2o4e/32.xyz`

by 1.37 kcal/mol at B3LYP/6-31g++(d).

Tip

We use `B3LYP/6-31++g(d)`

to study the system. However, since the addition electron in this system may be very delocalized, very diffuse basis functions and high level of theory (beyond MP2) are needed to get reasonable results. That is why we let `geom`

to generate GJF files for LMs. You can run more accurate calculations for the top 5 or 10 LMs to further refine the true GM.

## 6.4.2. Example: (CF3)4@C20

Tip

The sample input and output files can be found in `testfiles/geom/3-c20cf34-g16`

.

It is argued that carbon cage like \(\text{C}_{20}\) is an electron acceptor. So what will happen if some electron-withdrawing groups are attached? Let’s try to build a cluster of 4 \(\text{-CF}_3\) on \(\text{C}_{20}\).

The structures of \(\text{C}_{20}\) and \(\text{-CF}_{3}\) are:

```
120
2C20
3C -1.20472114 0.39143763 1.65815580
4C -0.74455861 -1.02479702 1.65815580
5C 1.20472114 0.39143763 1.65815580
6C -0.00000000 1.26671877 1.65815580
7C 0.00000065 2.04959354 0.39143762
8C -1.20471985 1.65815687 -0.39143704
9C 1.20472125 1.65815590 -0.39143685
10C 0.74456049 1.02479624 -1.65815544
11C -0.74455796 1.02479756 -1.65815575
12C 1.20472204 -0.39143886 -1.65815486
13C 1.94927928 -0.63336012 -0.39143623
14C 1.94927917 0.63335933 0.39143808
15C 0.74455861 -1.02479702 1.65815580
16C 1.20472070 -1.65815601 0.39143808
17C 0.00000038 -2.04959369 -0.39143685
18C -1.20472033 -1.65815639 0.39143762
19C -1.94927899 -0.63336053 -0.39143704
20C -1.20472093 -0.39143845 -1.65815575
21C 0.00000048 -1.26672037 -1.65815458
22C -1.94927868 0.63336060 0.39143848
```

```
14
2CF3
3C 0.00000000 0.00000000 -0.36818188
4F -1.10227044 -0.63639614 0.08181821
5F 0.00000000 1.27279228 0.08181821
6F 1.10227044 -0.63639614 0.08181821
```

The input file is:

```
1lm_dir c20cf34 # Save the local minima to this folder.
2num_calcs 5 # Total number of calculations.
3do_coarse_opt yes # no: Do NOT the coarse optimization.
4min_energy_gap 1.E-4 # When two energies differ smaller than
5 # this value, they are treated as identical.
6 # A negative number means do not remove
7 # energetically degenerated ones.
8max_geom_iters 3000 # The maximum number of iterations for local optimization.
9 # If it is less or equal than zero, then the number is unlimited.
10
11components
12 c20.xyz 1
13 fix 0 0 0 0 0 0
14 ****
15 cf3.xyz 4
16 shell 1 1 1 2 3 4 0 0 0 3 3 3
17 ****
18end
19
20commands
21 xyz2gaussian optfile $inp$ > $xxx$.gjf
22 g16 < $xxx$.gjf > $xxx$.log 2>/dev/null
23 gaussian2xyz $xxx$.log > $out$
24end
```

From `components`

, we can see that \(\text{C}_{20}\) is fixed at (`0`

, `0`

, `0`

) without rotation. Since its diameter is about 6 Å, we let the 4 \(\text{-CF}_{3}\) groups form a shell on an ellipsoid centered at (`0`

, `0`

, `0`

) with radii `3`

, `3`

, and `3`

. They point to the shell with an axis defined by the average of atom `1`

, `1`

, `1`

and `2`

, `3`

, `4`

. This system is very expensive, so we do only `5`

calculations.

The Gaussian template is:

```
1%nprocs=48
2%mem=20GB
3%chk=c20cf34.chk
4#N wB97XD/gen Opt
5
6opt
7
80 1
9>>>>
1021 25 29 33
116-31g(d)
12****
13F 0
146-31g(d)
15****
161 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
173-21g
18****
19>>>>
```

Here, we choose `wB97XD`

instead of `B3LYP`

because the latter usually performs poorly for aromatic systems. After the first `>>>>`

, we define the basis sets for different atoms to save the computational cost

Tip

This basis set canNOT be used in real scientific research. A better one is `def2-SVP`

. For property calculations, triple-zeta and diffusion functions may be required.

Now you can run the global optimization:

```
$ geom c20cf34.inp > c20cf34.out
```

After the optimization, the end of `c20cf34.out`

is

```
-- Result Report --
Results are energy-increasingly reordered.
Structures of energies within 1.000E-04 are treated as degenerate.
All minima are saved to "c20cf34".
-------------------------------------------------------------------
# index Energy NaiveRMSD
-------------------------------------------------------------------
0 2 -2107.38065472 0.00000000
1 0 -2107.38019685 1.83056800
2 1 -2107.37839045 3.12425180
3 3 -2107.37619597 3.38430408
4 4 -2106.56831104 2.73267963
-------------------------------------------------------------------
```

So the global minimum is `c20cf34/2.xyz`

, which is shown below.