4.2. Example: (H2O)6

Tip

The sample input and output files can be found in testfiles/rigidmol/1-h2o6.

4.2.1. Global Optimization

The system here is \(\left(\mathrm{H}_2\mathrm{O}\right)_{6}\). Its parameter file is misc/charmm36/tip4p.xyz.

Step 1: Prepare an input file named h2o6.inp with the following content:

h2o6.inp
1h2o6.cluster # cluster file name
220           # population size
320           # maximal generations
43            # scout limit
54.0          # amplitude
6h2o6         # save optimized configuration
730           # number of LMs to be saved

The texts after # are comments and can be arbitrary. I will explain the input line by line:

  1. The name of the file that contains cluster components and geometry. It will be explained in the next step.

  2. The population size \(SN\). See Theoretical Background for details.

  3. The maximum cycle number \(g_{\mathrm{max}}\). See Theoretical Background for details.

  4. The scout limit \(g_{\mathrm{limit}}\). See Theoretical Background for details.

  5. The estimated size of the cluster in Angstrom. See Theoretical Background for details.

  6. The folder name to save local minima.

  7. The number of local minima to be saved.

Step 2: Copy misc/charmm36/tip4p.xyz to the current path. Then prepare the cluster file named h2o6.cluster with the following content:

h2o6.cluster
11
2tip4p.xyz 6
3* 4.0000

  1. The kind of molecules in the cluster. For \(\left(\mathrm{H}_2\mathrm{O}\right)_{6}\), this is 1; For \(\mathrm{Na}^{+}\left(\mathrm{H}_2\mathrm{O}\right)_{3}\left(\mathrm{CH}_3\mathrm{OH}\right)_2\), this is 3.

  2. The parameter file name (with path, if necessary), and the number of this molecules in the cluster. If there are more than one kind of molecules in the cluster, then each molecule should be given by a line like this.

  3. Here, * 4.0000 means generating a random initial guess with \(L\) = 4.0000.

Tip

In cluster files, the parameter files can be given with relative or absolute path. For example, if the water parameter file is at /home/you/abcluster/misc/tip4p.xyz and you do not want to copy it, you can simply give the absolute path of this file in h2o6.cluster like /home/you/abcluster/misc/tip4p.xyz 6.

Note that, rigidmol will also generate *.cluster files, but it does not end with something like * 4.00 but with some real numbers, for example:

A cluster file generated by rigidmol
11
2tip4p.xyz 6
3Water hexamer
4 3.5  7.0  -0.2  1.3  2.9  0.6
5 4.0  0.9   3.1  4.3  3.2  3.8
6-1.8  4.4   1.0  0.8  0.2  0.6
7 4.2  3.4   3.9  1.9  0.2  3.2
8 3.6  4.7   2.9  3.6  1.7  1.7
9 5.1  5.2  -0.6  1.6  2.8  3.5

These numbers are rigid coordinates of each molecule in the cluster. You do not need to care about them.

Step 3: Run the global optimization:

$ rigidmol h2o6.inp > h2o6.out

After a few seconds, you will find several new files:

  • h2o6.out The main output file.

  • h2o6-OPT.cluster The global minimum in ABCluster format. It can only be read by rigidmol. Also, it can be used as the initial guess of a new optimization.

  • h2o6-OPT.xyz The global minimum in XYZ format. It can be read by, for example, VMD, CYLView, or VESTA.

  • h2o6-OPT.gjf The global minimum in Gaussian input format. It can be read by, for example, GaussView.

  • h2o6-LM A folder containing all local minima, each one having two files in XYZ, Gaussian input, and ABCluster format, respectively. They are sorted in energy-increasing order, e.g. 0.xyz is lower in energy than 13.xyz.

  • abcluster*.cluster/xyz/gjf/cluster The file containing the currently found global minimum during the running of rigidmol. You can check the current stable structure before rigidmol terminates. If it crashes, one can use this abcluster*.cluster to start a new optimization.

You will find the global minimum in h2o6-OPT.xyz and local minima in h2o6-LM. Below, the global minimum h2o6-LM/0.xyz (cage isomer) and an important local minimum h2o6-LM/2.xyz (prism isomer).

alternate text alternate text

4.2.2. Many-Body Energy Decomposition Analysis

Many-body energy decomposition analysis (MB-EDA) is an important method to analyze the cluster energy. Its detail can be seen in Theoretical Background and For Molecular Clusters: Many-Body Energy Decomposition Analysis. Here, we just show how to perform MB-EDA for h2o6-OPT.xyz.

  1. Get the quantum chemical software Qbics. You can download the executable or source code and compile by yourself from: https://zhjun-sci.com/qbics/. The details of compilation can be found here: https://zhjun-sci.com/qbics/doc/compilation.html.

  2. When Qbics is available, prepare the following input file:

h2o6-mbedas.inp
 1basis
 2   def2-svp
 3end
 4
 5scf
 6   charge     0 # Total charge.
 7   spin2p1    1
 8   type       U
 9end
10
11grimmedisp
12   type bj
13end
14
15eda
16   type      mb_tso
17   mb_level  4
18   frag 0 1 1-3
19   frag 0 1 4-6
20   frag 0 1 7-9
21   frag 0 1 10-12
22   frag 0 1 13-15
23   frag 0 1 16-18
24end
25
26mol
27   O      1.56089640      0.09521244      0.97451815
28   H      1.15030687      0.86428495      0.57932761
29   H      1.98904751     -0.35049295      0.24358370
30   O     -0.43542133     -0.58119050     -1.58243636
31   H     -0.69354569     -1.00007736     -0.76137762
32   H     -0.36958470      0.34910133     -1.36690329
33   O     -0.13166974      1.91014797     -0.32847358
34   H     -0.12458087      2.86059513     -0.44175107
35   H     -0.88872147      1.74093757      0.23230478
36   O     -2.17547367      1.03965222      1.23351847
37   H     -1.77664579      0.17396488      1.32157298
38   H     -3.11557395      0.86676020      1.18300454
39   O     -0.73246357     -1.30623952      1.12082542
40   H     -0.59931071     -2.14401038      1.56426893
41   H      0.11647565     -0.86758106      1.17661794
42   O      2.19540445     -1.19867926     -1.37339883
43   H      1.27354214     -1.04544150     -1.58057004
44   H      2.60372639     -1.38450287     -2.21896056
45end
46
47task
48  eda b3lyp
49end

I will explain the input line by line:

  1. Line 48: Means we will perform EDA with density functional theory B3LYP. You can also use other density functional like PBE if you wish.

  2. Line 2: Means we will use basis set def2-svp. You can also use 6-31g(d) or cc-pvdz or other basis sets you prefer.

  3. Line 6: The total charge is 0.

  4. Line 7: The total spin multiplicity is 1.

  5. Line 8: U is unrestricted SCF method. Only U can be used for EDA calculations.

  6. Line 12: BJ method will be used for Grimme dispersion correction.

  7. Line 16: The type of eda method. We strongly recommend to use mb_tso and do NOT change it.

  8. Line 17: This means the order to truncate for the level of many-body interaction. In most cases, 4 is enough, since higher orders are very small but computationally time-comsuming.

  9. Line 18: Define the first molecule fragment, with charge 0, spin multiplicity 1, and atoms 1-3. If atoms are not continuous, you can also write things like 1-5 9 12 34-37 which atom 1,2,3,4,5,9,12,34,35,36,37.

  10. Line 19-23: Define the other 5 molecular fragments.

  11. Line 136-153: Molecular geometry.

  1. Now, run Qbics to do the calculation:

$ qbics h2o6-mseda.inp -n 24 > h2o6-mseda.out

Here -n 24 means to parallelize the calculation with 24 cores.

  1. Now we can analyze the results. Open h2o6-mseda.out:

h2o6-mbedas.out
 1Table 5. Summary (kcal/mol).
 2---------------------------------------------------------------------------------------------------------------------------------
 3    Interactions         delE_el         delE_xc         delE_pl         delE_ct       delE_bsse       delE_disp        delE_tot
 4---------------------------------------------------------------------------------------------------------------------------------
 5   SUM of 2-body -7.42190602E+01  7.00852652E+01 -8.21424492E+00 -5.41899034E+01  3.41731354E+01 -7.86889683E+00 -4.02337048E+01
 6   SUM of 3-body -2.37203844E-09 -1.00133720E+00 -4.78274426E+00  2.27628396E+00 -5.35726511E+00  4.11748888E-03 -8.86094512E+00
 7   SUM of 4-body  3.83451522E-09  1.01702862E-01 -1.66028693E-01 -1.83396053E+00  1.65625886E+00 -3.61238203E-05 -2.42063621E-01
 8          Remain -2.81793653E-09 -1.19547650E-02 -1.83316504E-02  5.57013755E-02 -9.21206667E-02  1.89762141E-07 -6.67055197E-02
 9---------------------------------------------------------------------------------------------------------------------------------
10             SUM -7.42190602E+01  6.91736761E+01 -1.31813495E+01 -5.36918786E+01  3.03800084E+01 -7.86481528E+00 -4.94034190E+01
11---------------------------------------------------------------------------------------------------------------------------------

From this table, we can see: total interaction energy is -49.40 kcal/mol, with 2-, 3-, and 4-body energy being -40.23, -8.86, and -0.24 kcal/mol, respectively. The 3-body energy is negative, meaning that a cooperative effect exists in this cluster. The many-body interaction is stabilizing the cluster.

The total energy can be decomposed into 5 components: elecstrostatic (delE_el, -74.22 kcal/mol), Pauli-exclusion (delE_ex, 69.17 kcal/mol), porlarization (delE_pl, -13.18 kcal/mol), charge transfer (delE_ct``+ ``delE_bsse, -53.69+30.38 = -23.31 kcal/mol), and dispersion energy (delE_disp, -7.86 kcal/mol).

Each components can also be decomposed into many-body contributions. For example, the 2-, 3-, and 4-body order of polarization energy is -8.21, -4.78, and -0.17 kcal/mol. The 3-body contribution is significant. However, for electrostatic energy, the 3-body contribution is almost 0.

By carefully analyzing this MB-EDA results, one can study cluster formation mechanism, like in J. Comput. Chem. 2025, 46, e70114.