# AIQM1

New: we highly recommend you to check out UAIQM that is **our ultimate solution to the universal ML models** and considered as successor to AIQM1. The only advantages of AIQM1 at the moment are: support of excited states and calibrated, accurate heats of formation, and possibility to install locally. On the XACS cloud AIQM1 has many numerical problems due to unavailability of the analytical derivatives on the cloud version.

AIQM1 (artificial intelligence–quantum mechanical method 1) is a general-purpose method approaching the gold-standard coupled cluster quantum mechanical method with high computational speed of the approximate low-level semiempirical quantum mechanical methods for the ground-state, closed-shell species, but also transferable for calculation of charged and radical species as well as for excited-state calculations with a good accuracy. See AIQM1 paper for more details. Please cite this paper alongside other required citations:

Peikun Zheng, Roman Zubatyuk, Wei Wu, Olexandr Isayev, Pavlo O. Dral. Artificial Intelligence-Enhanced Quantum Chemical Method with Broad Applicability.

*Nat. Commun.***2021**,*12*, 7022, DOI: 10.1038/s41467-021-27340-2.

**Strengths:** AIQM1 is especially good for energy calculations and geometry optimizations of closed-shell molecules in their ground-state.

**Limitations:** This method is currently limited to compounds only containing H, C, N, and O elements.

Below we provide important notes on installation of AIQM1 and non-exhaustive examples of its uses via command line/input files and Python API for:

calculation of thermochemical properties

special input for semi-empirical part

calculations of charged species and radicals

## Installation

Refer to MLatom installation instructions. AIQM1 requires TorchANI and dftd4 as well as a semi-empirical program. The semi-empirical program can be MNDO which supports analytical gradients and excited-state calculations. On MLatom@XACS cloud computing service we cannot use MNDO but use development version of SCINE Sparrow (see also a paper on Sparrow). This version does not support analytical gradients and no excited-state calculations with AIQM1. Hence, Sparrow can be used with numerical gradients for geometry optimizations and MD but might lead to numerical instabilities in, e.g., frequency calculations. Some features described below might require interfaces to the Gaussian or ASE or PySCF packages. Installation instructions, a usage manual, and proper citations are given below.

## Single-point calculations

Single-point calculations of energies (and gradients if needed) of closed-shell molecules in electronic ground state is the simplest job, which can be run with 3-4 line MLatom input file, e.g., `sp.inp`

:

```
AIQM1 # or AIQM1@DFT or AIQM1@DFT* if you want to use these methods
xyzfile=sp.xyz
yestfile=enest.dat
ygradxyzestfile=gradest.dat
```

This input requires a `sp.xyz`

file with XYZ geometries of molecules (you can provide many molecules as usual for MLatom), e.g., for hydrogen and methane `sp.xyz`

file can look like (geometries in Å):

```
2
H 0.000000 0.000000 0.363008
H 0.000000 0.000000 -0.363008
5
C 0.000000 0.000000 0.000000
H 0.627580 0.627580 0.627580
H -0.627580 -0.627580 0.627580
H 0.627580 -0.627580 -0.627580
H -0.627580 0.627580 -0.627580
```

After you prepared your input files `sp.inp`

and `sp.xyz`

, you can run MLatom as usual:

```
mlatom sp.inp > sp.out
```

After the calculations finish, MLatom output `sp.out`

will contain the standard deviation of NN prediction and components of AIQM1 energies:

```
Standard deviation of NN contribution : 0.00892407 Hartree 5.59994 kcal/mol
NN contribution : -0.00210898 Hartree
Sum of atomic self energies : -0.08587317 Hartree
ODM2* contribution : -1.09094119 Hartree
D4 contribution : -0.00000889 Hartree
Total energy : -1.17893224 Hartree
Standard deviation of NN contribution : 0.00025608 Hartree 0.16069 kcal/mol
NN contribution : 0.00958812 Hartree
Sum of atomic self energies : -33.60470494 Hartree
ODM2* contribution : -6.86968756 Hartree
D4 contribution : -0.00010193 Hartree
Total energy : -40.46490632 Hartree
```

You can also get AIQM1 energies saved in file `enest.dat`

looking like this (energies in Hartree):

```
-1.178932238420
-40.464906315250
```

and gradients saved in file `gradest.dat`

looking like this (gradients in Hartree/Å):

```
2
0.000000000000 0.000000000000 0.000032023551
0.000000000000 0.000000000000 -0.000032023551
5
-0.000000000000 -0.000000000000 0.000000000000
0.000490470799 0.000490470714 0.000490470881
-0.000490470799 -0.000490470714 0.000490470881
0.000490470799 -0.000490470714 -0.000490470881
-0.000490470799 0.000490470714 -0.000490470881
```

You can also perform the same calculation in Python API. The code snippet looks like:

```
import mlatom as ml
# read molecule from .xyz file
molDB = ml.data.molecular_database.from_xyz_file('sp.xyz')
# define method
model = ml.models.methods(method='AIQM1', qm_program='MNDO')
# perform the calculations
model.predict(molecular_database=molDB, calculate_energy_gradients=True, calculate_hessian=True)
# print the energy. You can print extract properties from each molecule object, such as energy gradients and Hessian
print(molDB[0].energy)
print(molDB[1].energy)
print(molDB[1].energy-molDB[0].energy)
```

The examples above were tested with the interface to the MNDO program. If you use the MLatom@XACS cloud computing service, it runs calculations using the interface to the SCINE Sparrow package, thus, the numbers may be slightly different. See tutorial for more details on single point energy calculations.

## Geometry optimizations

Geometry optimization of closed-shell molecules in electronic ground state is as simple as running single point calculations and 4-line MLatom input file, e.g., `opt.inp`

, looks like this:

```
AIQM1
xyzfile=opt.xyz
optxyz=final_geo.xyz
geomopt
optprog=gaussian # or optprog=ase if you choose ASE
```

This input requires `opt.xyz`

file with initial XYZ geometries of molecules to be optimized (you can provide many molecules as usual for MLatom), e.g., for methane `opt.xyz`

file can look like (geometries in Å):

```
5
C 0.0000000000 0.0000000000 0.0000000000
H 1.0870000000 0.0000000000 0.0000000000
H -0.3623333220 -1.0248334322 -0.0000000000
H -0.3623333220 0.5124167161 -0.8875317869
H -0.3623333220 0.5124167161 0.8875317869
```

After you prepared your input files `opt.inp`

and `opt.xyz`

, you can run MLatom as usual:

```
mlatom opt.inp
```

After the calculations finish, the optimized geometries are saved in either Gaussian or ASE-style output files.

If you use Gaussian, the Gaussian output files of geometry optimizations are saved in `mol_1.log`

, `mol_2.log`

, … files for each molecule, where you can find the optimized geometries and visualize them by GausianView program as in any other typical Gaussian job.

If you use ASE, you can directly obtain optimized XYZ geometries in a single file `final_geo.xyz`

with MLatom format, which for our example looks like (geometries in Å):

```
5
C 0.00000000 0.00000000 0.00000000
H 1.08666998 -0.00000000 0.00000000
H -0.36222332 -1.02452229 -0.00000000
H -0.36222332 0.51226114 -0.88726233
H -0.36222332 0.51226114 0.88726233
```

In Python API, the code snippet looks like:

```
import mlatom as ml
# Get the initial guess for the molecules to optimize
initmol = ml.data.molecule.from_xyz_file('opt.xyz')
# Choose one of the predifined (automatically recognized) methods
mymodel = ml.models.methods(method='AIQM1')
# Optimize the geometry with the choosen optimizer:
geomopt = ml.optimize_geometry(model=mymodel,
initial_molecule=initmol,
program='ASE',
maximum_number_of_steps=10000)
# Get the final geometry
final_mol = geomopt.optimized_molecule
# and its XYZ coordinates
final_mol.xyz_coordinates
# Check and save the final geometry
print('Optimized coordinates:')
print(final_mol.get_xyz_string())
final_mol.write_file_with_xyz_coordinates(filename='final_geo.xyz')
```

See tutorial for more details on geometry optimization.

## TS optimization

The geometry optimization of transition states structures (TSs) is almost the same as above with the only difference that we replace `geomopt`

with `ts`

keyword. These calculations require calculations of Hessians and thus are more expensive. Currently, they can only be performed with Gaussian interface, i.e., use `optprog=gaussian`

in the input file. Example of input file `ts.inp`

:

```
AIQM1
xyzfile=ts.xyz
optxyz=final_geo.xyz
ts
optprog=gaussian
```

The xyz file `ts.xyz`

required by the input file should be the initial guess for TS. For a pericyclic reaction:

the file `ts.xyz`

can look like:

```
14
C -0.12282076 -1.09054053 1.14065341
C 0.17199635 0.23301894 1.48248953
C -0.12282076 1.35549534 0.70296425
C -0.12282076 1.35549534 -0.70296425
C 0.17199635 0.23301894 -1.48248953
C -0.12282076 -1.09054053 -1.14065341
H 0.82454258 0.40945200 2.34048883
H -0.03331673 2.33317669 1.17510642
H -0.03331673 2.33317669 -1.17510642
H 0.82454258 0.40945200 -2.34048883
H 0.43031244 -1.88602881 -1.64114758
H -1.14789312 -1.35457363 -0.93625533
H 0.43031244 -1.88602881 1.64114758
H -1.14789312 -1.35457363 0.93625533
```

When the `ts.inp`

and `ts.xyz`

are both ready, the TS geometry optimization could be conducted via running MLatom as usual:

```
mlatom ts.inp
```

After the calculations finish, the optimized TS geometries will be saved in the XYZ output file `final_geo.xyz`

. You can check the TS optimization by investigating frequencies: the number of imaginary frequencies for successful TS optimization should be one and analyze the direction of nuclear motions by visualizing Gaussian output `.log`

files with your favorite software (e.g., GaussView).

```
14
6 -0.14824052 -1.18232681 1.09011331
6 0.27697295 0.08303553 1.40719277
6 -0.14824052 1.26264088 0.68226088
6 -0.14824052 1.26264088 -0.68226088
6 0.27697295 0.08303553 -1.40719277
6 -0.14824052 -1.18232681 -1.09011331
1 1.16206956 0.21862475 2.03298866
1 -0.22714441 2.21025461 1.20566513
1 -0.22714441 2.21025461 -1.20566513
1 1.16206956 0.21862475 -2.03298866
1 0.37749585 -2.04342303 -1.49163954
1 -1.19537249 -1.36555399 -0.92796582
1 0.37749585 -2.04342303 1.49163954
1 -1.19537249 -1.36555399 0.92796582
```

In Python API, the code snippet looks like:

```
import mlatom as ml
# Get the initial guess for the molecules to optimize
initmol = ml.data.molecule.from_xyz_file('ts.xyz')
# Choose one of the predifined (automatically recognized) methods
mymodel = ml.models.methods(method='AIQM1')
# Optimize the geometry with the choosen optimizer:
geomopt = ml.optimize_geometry(model=mymodel,
initial_molecule=initmol,
ts=True,
program='ASE',
maximum_number_of_steps=10000)
# Get the final geometry
final_mol = geomopt.optimized_molecule
# and its XYZ coordinates
final_mol.xyz_coordinates
# Check and save the final geometry
print('Optimized coordinates:')
print(final_mol.get_xyz_string())
final_mol.write_file_with_xyz_coordinates(filename='final_geo.xyz')
```

See tutorial for more details on geometry optimization of transition states structures.

## Thermochemical calculations

Thermochemical properties of closed-shell molecules in electronic ground state can be calculated at AIQM1 level by adding an option freq to the MLatom input file, e.g., `freq.inp`

and they are typically run on AIQM1-optimized geometries. An example of MLatom input file:

```
AIQM1
xyzfile=freq.xyz
freq
optprog=gaussian # or optprog=ase if you choose ASE
```

This input requires `freq.xyz`

file with initial XYZ geometries of molecules (you can provide many molecules as usual for MLatom), e.g., for methane `freq.xyz`

file can look like (geometries in Å):

```
5
C 0.000000 0.000000 0.000000
H 0.627580 0.627580 0.627580
H -0.627580 -0.627580 0.627580
H 0.627580 -0.627580 -0.627580
H -0.627580 0.627580 -0.627580
```

After you prepared your input files `freq.inp`

and `freq.xyz`

, you can run MLatom as usual:

```
mlatom freq.inp > freq.out
```

After the calculations finish, MLatom output `freq.out`

will contain the summary with atomization enthalpy at 0 K, ZPVE-exclusive atomization energy at 0 K, and heat of formation at 298 K for the molecule.

If you use Gaussian, the Gaussian output files of frequency calculations are saved in `mol_1.log`

, `mol_2.log`

, … files for each molecule; these files contain ZPVE energy and lots of thermochemical data such as entropy and the Gibbs free energy. MLatom output will contain additional information as described above with following lines for methane example:

```
Standard deviation of NN contribution : 0.00025608 Hartree 0.16069 kcal/mol
NN contribution : 0.00958812 Hartree
Sum of atomic self energies : -33.60470494 Hartree
ODM2* contribution : -6.86968756 Hartree
D4 contribution : -0.00010193 Hartree
Total energy : -40.46490632 Hartree
Atomization enthalpy at 0 K : 391.58894 kcal/mol
ZPE exclusive atomization energy at 0 K : 419.90907 kcal/mol
Heat of formation at 298.15 K : -17.30543 kcal/mol
```

When you use ASE for the calculation of thermochemical properties, you should specify `ase.linear`

and `ase.symmetrynumber`

this two keywrods. `ase.linear`

is 0 for nonlinear molecule, 1 for linear molecule, and `ase.symmetrynumber`

is the rotational symmetry number for each molecule, (see Table 10.1 and Appendix B of C. Cramer “Essentials of Computational Chemistry”, 2nd Ed.). For example, for hydrogen and methane this two molecules, you should set `ase.linear=1,0`

and `ase.symmetrynumber=2,12`

.

The Python API provides more flexibility to use more complicated models and build custom workflows. Below is a workflow of optimizing geometry and calculating frequency.

```
import mlatom as ml
# Get the initial guess for the molecules to optimize
initmol = ml.data.molecule.from_xyz_file('freq.xyz')
# Choose one of the predifined (automatically recognized) methods
mymodel = ml.models.methods(method='AIQM1')
# Optimize the geometry with the choosen optimizer:
geomopt = ml.optimize_geometry(model=mymodel,
initial_molecule=initmol,
program='ASE',
maximum_number_of_steps=10000)
# Get the optimized molecule
final_mol = geomopt.optimized_molecule
# Do vibration analysis and thermochemistry calculation
freq = ml.thermochemistry(model=mymodel,
molecule=final_mol,
program='ASE')
# or vibration analysis only
# freq = ml.freq(model=mymodel, molecule=final_mol, program='')
# Save the molecule with vibration analysis and thermochemistry results
final_mol.dump(filename='final_mol.json',format='json')
# Check vibration analysis
print("Mode Frequencies Reduced masses Force Constants")
print(" (cm^-1) (AMU) (mDyne/A)")
for ii in range(len(final_mol.frequencies)):
print("%d %13.4f %13.4f %13.4f"%(ii,final_mol.frequencies[ii],final_mol.reduced_masses[ii],final_mol.force_constants[ii]))
# Check thermochemistry results
print(f"Zero-point vibrational energy: {final_mol.ZPE} Hartree")
print(f"Enthalpy at 298 K: {final_mol.H} Hartree")
print(f"Gibbs Free energy at 298 K: {final_mol.G} Hartree")
print(f"Heat of formation at 298 K: {final_mol.DeltaHf298} Hartree")
```

See tutorial for more details on frequency and thermochemistry calculations.

## Special input for semi-empirical programs

You need to use `QMprogramKeywords`

Sometimes, you may need to change the semiempirical quantum mechanical part of AIQM1 by informing the MNDO or Sparrow program. MLatom implementation therefore supports the request to read user-defined MNDO or Sparrow keywords via option `QMprogramKeywords=[file name with MNDO or Sparrow keywords]`

, which are passed to the MNDO program. For historical reasons, `mndokeywords`

is also supported for the MNDO program. Please consult the MNDO program manual for the available keywords.

Note

Whenever you request special keywords for MNDO, keywords `iop=-22 immdp=-1 nsav15=3 igeom=1 iform=1`

are always required and AIQM1 calculations can be run only for a single molecule in xyz file.

## Calculations of charged species and radicals

For example, if we want to optimize geometry of protonated water H_{3}O^{+}, then you can use the following MLatom input `optcs.inp`

:

```
AIQM1
xyzfile=optcs.xyz
optxyz=final_geo.xyz
geomopt
optprog=ase
charges=1
```

Similarly, for radicals, you can provide the keyword `multiplicities`

.

This input requires `optcs.xyz`

file with initial XYZ geometries of molecules to be optimized, e.g. (geometries in Å):

```
4
O 0.00000000 0.00000000 0.08727273
H 0.00000000 0.90509668 -0.23272727
H -0.78383672 -0.45254834 -0.23272727
H 0.78383672 -0.45254834 -0.23272727
```

After you prepared your input files `optcs.inp`

and `optcs.xyz`

, you can run MLatom as usual:

```
mlatom optcs.inp
```

The output geometry is saved in `final_geo.xyz`

with MLatom format, which for our example looks like (geometries in Å):

```
4
O 0.00000000 -0.00000000 0.12401796
H 0.00000000 0.90385318 -0.24497568
H -0.78275982 -0.45192659 -0.24497568
H 0.78275982 -0.45192659 -0.24497568
```

When you want to calculate the heat of formation of H_{3}O^{+} with the above optimized geometry, then you can use the following MLatom input `freqcs.inp`

:

```
AIQM1
xyzfile=final_geo.xyz
freq
optprog=ase
charges=1
```

After you run MLatom by `mlatom freqcs.inp > freqcs.out`

, you can get the following heat of formation from the output:

```
Heat of formation at 298.15 K : 150.84947 kcal/mol
* Warning * Heat of formation have high uncertainty!
```

## Vertical excitation energies

For example, if we want to calculate the vertical excitation energies of ethene in ^{1}B_{1u} state, then you can use the following MLatom input `vee.inp`

:

```
AIQM1
xyzfile=vee.xyz
mndokeywords=mndokw
```

This input requires a `vee.xyz`

file with XYZ geometries of molecules, e.g. (geometries in Å):

```
6
Ethene 1B1u
C 0.000000 0.000000 0.668188
C 0.000000 0.000000 -0.668188
H 0.000000 0.923274 1.238289
H 0.000000 -0.923274 1.238289
H 0.000000 0.923274 -1.238289
H 0.000000 -0.923274 -1.238289
```

and `mndokw`

file with MNDO keywords:

```
iop=-22 immdp=-1 +
igeom=1 iform=1 +
jop=-2 nsav15=3 +
kci=5 ici1=1 ici2=1 iroot=2 ioutci=2 +
movo=-1 nciref=1 mciref=3 levexc=2
```

After you prepared your input files `vee.inp`

, `vee.xyz`

, and `mndokw`

you can run MLatom as usual:

```
mlatom vee.inp
```

The output of vertical excitation energy is saved in the MNDO program outfile `mndo.out`

, where you can see:

```
State 2, Mult. 1, B1u (4), E-E(1)= 7.910791 eV, E= -304.054934 eV
```

See also tutorial.

## Excited-state geometry optimization

For example, if we want to optimize geometry of formaldehyde in ^{1}nπ^{*} excited state (^{1}A”), then you can use the following MLatom input `optes.inp`

:

```
AIQM1
xyzfile=optes.xyz
optxyz=final_geo.xyz
geomopt
optprog=ase
mndokeywords=mndokw_es
```

This input requires `optes.xyz`

file with initial XYZ geometries of molecules to be optimized, e.g. (geometries in Å):

```
4
formaldehyde 1npi*
C -0.02667227 0.64339915 -0.00000000
O -0.02667227 -0.71674790 0.00000000
H 0.18670592 0.93679416 1.04362466
H 0.18670592 0.93679416 -1.04362466
```

and `mndokw_es`

file with MNDO keywords:

```
iop=-22 immdp=-1 +
igeom=1 iform=1 +
jop=-2 nsav15=3 +
kharge=0 imult=0 +
kci=5 ici1=6 ici2=4 jci1=1 jci2=1 ncisym=2 +
movo=-1 nciref=1 mciref=3 levexc=2 iroot=1 lroot=1
```

The first line in `mndokw_es`

file requests the use of the ODM2* Hamiltonian (you do not need to modify it), the second line specifies that geometry is given in a free XYZ format, the third line requests calculation of gradients and saving energies and gradients into `fort.15`

file (required for geometry optimization), the fourth line sets charge (kharge=0) and multiplicity (imult=0), and the following lines setup the type of excited-state calculations. `kci=5`

requests GUGA-CI approach for excited state calculation; `ici1`

and `ici2`

define the number of active occupied orbitals and unoccupied orbitals respectively; `jci1`

and `jci2`

define the number of occupied and unoccupied pi-MOs included in the active space; `ncisym`

define the state symmetry. Here the initial geometry has C^{s}point group, which has A’ and A” this two irreducible representation, and we use `ncisym=2`

to calculate the ^{1}A” state. `nciref`

and `mciref`

define the number of reference occupations and definition of reference occupations respectively; `levexc`

define the maximum excitation level; `iroot`

and `lroot`

is the number of lowest CI states computed and the number of the CI state of interest.

After you prepared your input files `optes.inp`

, `optes.xyz`

, and `mndokw_es`

you can run MLatom as usual:

```
mlatom optes.inp
```

The output geometry is saved in `final_geo.xyz`

with MLatom format, which for our example looks like (geometries in Å):

```
4
C -0.10906735 0.55934608 -0.00000000
O -0.02784720 -0.77782274 0.00000000
H 0.22849093 1.00935811 0.92370071
H 0.22849093 1.00935811 -0.92370071
```

## Visualizing molecular orbitals

If you use Sparrow for providing the ODM2* part: `vmo_sparrow.inp`

, `sp.xyz`

```
AIQM1
qmprog=sparrow
xyzfile=sp.xyz
yestfile=enest.dat
QMprogramKeywords=SparrowKW
```

You need to use `QMprogramKeywords`

for printing orbitals in Molden wavefunction format. The content of the file `SparrowKW`

should be `-W`

. After calculation, the `wavefunction.molden.input`

file will be generated. You can use Jmol or MOLDEN program to open this file and see the orbitals. For Jmol program, to get the pictures of orbitals. You should right-click in the panel and select `Surface`

and `MO`

. For MOLDEN program, to get the pictures of orbitals, you need to choose Dens. Mode from the control panel and select the `Orbital`

button, then select the `Space`

option to set the contour value. Besides, you can also use Chemcraft to visualize the orbitals, but you need to rename this file with a `.molden`

suffix.

If you use MNDO for providing the ODM2* part: `vmo_mndo.inp`

, `sp.xyz`

```
AIQM1
qmprog=mndo
xyzfile=sp.xyz
yestfile=enest.dat
mndokeywords=mndokw_vmo
```

You should provide the `mndokw_vmo`

file with general MNDO keywords and add an extra keyword `nsav13=2`

:

```
iop=-22 immdp=-1 +
igeom=1 iform=1 +
jop=-1 nsav15=3 +
kharge=0 imult=0 nsav13=2
```

Then the `molden.dat`

file will be generated. You can also use Jmol or MOLDEN program to visualize the orbitals.

## Citations

If you have used AIQM1, then the following citations are appropriate in your publication:

**AIQM1:**P. Zheng, R. Zubatyuk, W. Wu, O. Isayev, P. O. Dral. Artificial Intelligence-Enhanced Quantum Chemical Method with Broad Applicability.*Nat. Commun.***2021**,*12*, 7022. DOI: 10.1038/s41467-021-27340-2.**MLatom:**Pavlo O. Dral, Fuchun Ge, Yi-Fan Hou, Peikun Zheng, Yuxinxin Chen, Mario Barbatti, Olexandr Isayev, Cheng Wang, Bao-Xin Xue, Max Pinheiro Jr, Yuming Su, Yiheng Dai, Yangtao Chen, Lina Zhang, Shuang Zhang, Arif Ullah, Quanhao Zhang, Yanchi Ou. MLatom 3: A Platform for Machine Learning-enhanced Computational Chemistry Simulations and Workflows.*J. Chem. Theory Comput.***2024**,*20*, 1193–1213. DOI: 10.1021/acs.jctc.3c01203.**MLatom:**Pavlo O. Dral, Fuchun Ge, Yi-Fan Hou, Peikun Zheng, Yuxinxin Chen, Bao-Xin Xue, Max Pinheiro Jr, Yuming Su, Yiheng Dai, Yangtao Chen, Shuang Zhang, Lina Zhang, Arif Ullah, Quanhao Zhang, Yanchi Ou. MLatom: A Package for Atomistic Simulations with Machine Learning, version [add version number], Xiamen University, Xiamen, China, 2013–2024. MLatom.com.**ODM2* Hamiltonian:**P. O. Dral, X. Wu, W. Thiel,*J. Chem. Theory Comput.***2019**,*15*, 1743.If you use MNDO program:

**MNDO program:**W. Thiel,`MNDO`

[check your version] (Max-Planck-Institut fur Kohlenforschung, Mulheim an der Ruhr, [year])If you use MLatom@XACS cloud computing service, where SCINE Sparrow is used:

F. Bosia, T. Husch, C. H. Müller, S. Polonius, J.-G. Sobez, M. Steiner, J. P. Unsleber, A. C. Vaucher, T. Weymuth, M. Reiher,

`qcscine/sparrow`

: to-be-released alpha version with modifications by P. Zheng,**2022**.T. Husch, A. C. Vaucher, M. Reiher,

*Int. J. Quantum Chem.***2018**,*118*, e25799.

**D4:**E. Caldeweyher, C. Bannwarth, S. Grimme,*J. Chem. Phys.***2017**,*147*, 034112.**D4 program:**E. Caldeweyher, S. Ehlert, S. Grimme, DFT-D4, Version [check your version], (Mulliken Center for Theoretical Chemistry, University of Bonn, [year])**ANI model:**J. S. Smith, O. Isayev, A. E. Roitberg,*Chem. Sci.***2017**,*8*, 3192**TorchANI program:**X. Gao, F. Ramezanghorbani, O. Isayev, J. S. Smith, A. E. Roitberg,*J. Chem. Inf. Model.***2020**,*60*, 3408