Methods
There are different methods available in MLatom, visit http://mlatom.com/docs/api.html#methods for more details. Here we use three methods for geometry optimization and vibration analysis of hydrogen and vinyl acetylene.
Before optimizing, the initial geometries of hydrogen and vinyl acetylene (h2_vinylacetylene_init.xyz)
should be provided.
# h2_vinylacetylene_init.xyz
2
H 0.0000000000000 0.0000000000000 0.0000000000000
H 0.0000000000000 0.0000000000000 0.8000000000000
8
C 0.00000 0.00000 0.00000
C 1.33718 0.00000 0.00000
C 2.04862 1.23465 0.00000
C 2.65517 2.30065 0.00000
H 3.17326 3.20500 0.00000
H -0.56880 0.91919 0.00000
H -0.56652 -0.92292 -0.00000
H 1.95564 -0.90366 0.00000
ANI-1ccx
ANI-1ccx is a pure machine learning method, which approaches CCSD(T)/CBS
accuracy on benchmarks for reaction thermochemistry, isomerization, and drug-like molecular torsions. And ANI-1ccx is
orders of magnitude faster than common DFT methods. For command line, prepare the input file of geometry optimization
(h2_vinylacetylene_opt_ANI1ccx.inp).
# h2_vinylacetylene_opt_ANI1ccx.inp
ANI-1ccx
geomopt
xyzfile=h2_vinylacetylene_init.xyz
optxyz=final_ANI1ccx.xyz
optprog=ase
Perform geometry optimization. And the optimized geometry will be saved in final_ANI1ccx.xyz.
mlatom h2_vinylacetylene_opt_ANI1ccx.inp
There is the code of geometry optimization in Python.
import mlatom as ml
from mlatom import constants
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('h2_vinylacetylene_init.xyz')
# define the model with methods
model = ml.models.methods(method='ANI-1ccx')
for mol in molDB:
print(mol.get_xyz_string())
# run geometry optimization
ml.optimize_geometry(model=model, molecule=mol, program='ASE')
print(mol.get_xyz_string())
To ensure the optimized geometry is a minimum rather than a saddle point, we need to do vibration analysis. Now we use
the optimized geometry final_ANI1ccx.xyz in the previous step to do vibration analysis. We need to prepare the input
file: h2_vinylacetylene_freq_ANI1ccx.inp.
# h2_vinylacetylene_freq_ANI1ccx.inp
ANI-1ccx
freq
xyzfile=final_ANI1ccx.xyz
optprog=ase
ase.linear=1,0
Run it.
mlatom h2_vinylacetylene_freq_ANI1ccx.inp
In Python, we need to provide the optimized geometry final_ANI1ccx.xyz.
import mlatom as ml
from mlatom import constants
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('final_ANI1ccx.xyz')
molDB[0].shape = 'linear'
# define the model with methods
model = ml.models.methods(method='ANI-1ccx')
for mol in molDB:
ml.thermochemistry(model=model, molecule=mol, program='ASE')
fmt = ' %-41s: %15.5f Hartree'
print(fmt % ('Standard deviation of NNs', mol.ani1ccx.energy_standard_deviation), end='')
print(' %15.5f kcal/mol' % (mol.ani1ccx.energy_standard_deviation * constants.Hartree2kcalpermol))
print(fmt % ('ZPE-exclusive internal energy at 0 K', mol.energy))
print(fmt % ('Zero-point vibrational energy', mol.ZPE))
print(fmt % ('Internal energy at 0 K', mol.U0))
print(fmt % ('Enthalpy at 298 K', mol.H))
print(fmt % ('Gibbs free energy at 298 K', mol.G))
if 'DeltaHf298' in mol.__dict__:
print('')
fmt = ' %-41s: %15.5f Hartree %15.5f kcal/mol'
print(fmt % ('Atomization enthalpy at 0 K', mol.atomization_energy_0K, mol.atomization_energy_0K * constants.Hartree2kcalpermol))
print(fmt % ('ZPE-exclusive atomization energy at 0 K', mol.ZPE_exclusive_atomization_energy_0K, mol.ZPE_exclusive_atomization_energy_0K * constants.Hartree2kcalpermol))
print(fmt % ('Heat of formation at 298 K', mol.DeltaHf298, mol.DeltaHf298 * constants.Hartree2kcalpermol))
if mol.ani1ccx.energy_standard_deviation > 1.68*constants.kcalpermol2Hartree:
print(' * Warning * Heat of formation have high uncertainty!')
AIQM1
AIQM1 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. And it’s also orders of magnitude faster than common DFT methods.
To use AIQM1, extra programs are required. Please move from Colab to XACS cloud (https://xacs.xmu.edu.cn/newcloud/).
Necessary files can be found in /export/home/xacscloud/workshop2023/2023_workshop_Torun/
The same as above, we prepare the input file h2_vinylacetylene_opt_AIQM1.inp in command line.
# h2_vinylacetylene_opt_AIQM1.inp
AIQM1
geomopt
xyzfile=h2_vinylacetylene_init.xyz
optxyz=final_AIQM1.xyz
optprog=ase
Run it. And the optimized geometry will be saved in final_AIQM1.xyz.
mlatom h2_vinylacetylene_opt_AIQM1.inp
There is the code of geometry optimization in Python.
import mlatom as ml
from mlatom import constants
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('h2_vinylacetylene_init.xyz')
molDB[0].shape = 'linear'
# define the model with methods
model = ml.models.methods(method='AIQM1')
for mol in molDB:
print(mol.get_xyz_string())
# run geometry optimization
ml.optimize_geometry(model=model, molecule=mol, program='ASE')
print(mol.get_xyz_string())
Now we use the optimized geometry final_AIQM1.xyz in the previous step to do vibration analysis in command line.
We need to prepare the input file: h2_vinylacetylene_freq_AIQM1.inp.
# h2_vinylacetylene_freq_AIQM1.inp
AIQM1
freq
xyzfile=final_AIQM1.xyz
optprog=ase
ase.linear=1,0
Run it.
mlatom h2_vinylacetylene_freq_AIQM1.inp
In Python, we need to provide the optimized geometry final_AIQM1.xyz.
import mlatom as ml
from mlatom import constants
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('final_AIQM1.xyz')
molDB[0].shape = 'linear'
# define the model with methods
model = ml.models.methods(method='AIQM1')
for mol in molDB:
# run thermochemistry calculations
ml.thermochemistry(model=model, molecule=mol, program='ASE')
fmt = ' %-41s: %15.5f Hartree'
print(fmt % ('Standard deviation of NN contribution', mol.aiqm1_nn.energy_standard_deviation), end='')
print(' %15.5f kcal/mol' % (mol.aiqm1_nn.energy_standard_deviation * constants.Hartree2kcalpermol))
print(fmt % ('NN contribution', mol.aiqm1_nn.energy))
print(fmt % ('Sum of atomic self energies', mol.aiqm1_atomic_energy_shift.energy))
print(fmt % ('ODM2* contribution', mol.odm2star.energy))
print(fmt % ('D4 contribution', mol.d4_wb97x.energy))
print(fmt % ('Total energy', mol.energy))
print(fmt % ('ZPE-exclusive internal energy at 0 K', mol.energy))
print(fmt % ('Zero-point vibrational energy', mol.ZPE))
print(fmt % ('Internal energy at 0 K', mol.U0))
print(fmt % ('Enthalpy at 298 K', mol.H))
print(fmt % ('Gibbs free energy at 298 K', mol.G))
if 'DeltaHf298' in mol.__dict__:
print('')
fmt = ' %-41s: %15.5f Hartree %15.5f kcal/mol'
print(fmt % ('Atomization enthalpy at 0 K', mol.atomization_energy_0K, mol.atomization_energy_0K * constants.Hartree2kcalpermol))
print(fmt % ('ZPE-exclusive atomization energy at 0 K', mol.ZPE_exclusive_atomization_energy_0K, mol.ZPE_exclusive_atomization_energy_0K * constants.Hartree2kcalpermol))
print(fmt % ('Heat of formation at 298 K', mol.DeltaHf298, mol.DeltaHf298 * constants.Hartree2kcalpermol))
if mol.aiqm1_nn.energy_standard_deviation > 0.41*constants.kcalpermol2Hartree:
print(' * Warning * Heat of formation have high uncertainty!')
B3LYP
B3LYP is a common DFT method which is widely used. Now we use this method to optimize the geometry. Prepare the input
file fistly: h2_vinylacetylene_opt_B3LYP.inp
# h2_vinylacetylene_opt_B3LYP.inp
method=B3LYP/6-31G*
QMprog=PySCF
geomopt
xyzfile=h2_vinylacetylene_init.xyz
optxyz=final_B3LYP.xyz
optprog=ase
Run it. The optimized geometry will be saved in final_B3LYP.xyz.
mlatom h2_vinylacetylene_opt_B3LYP.inp
There is the code of geometry optimization in Python.
import mlatom as ml
import numpy as np
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('h2_vinylacetylene_init.xyz')
molDB[0].shape = 'linear'
# define the model with methods
model = ml.models.methods(method='B3LYP/6-31G*', program='PySCF')
for mol in molDB:
print(mol.get_xyz_string())
# run geometry optimization
ml.optimize_geometry(model=model, molecule=mol, program='ASE')
print(mol.get_xyz_string())
Now we use the optimized geometry final_B3LYP.xyz in the previous step to do vibration analysis in command line.
Prepare the input file: h2_vinylacetylene_freq_B3LYP.inp.
# h2_vinylacetylene_freq_B3LYP.inp
method=B3LYP/6-31G*
QMprog=PySCF
freq
xyzfile=final_B3LYP.xyz
optprog=ase
ase.linear=1,0
Run it.
mlatom h2_vinylacetylene_freq_B3LYP.inp
In Python, we need to provide the optimized geometry final_B3LYP.xyz.
import mlatom as ml
import numpy as np
# load initial geometries
molDB = ml.data.molecular_database.from_xyz_file('final_B3LYP.xyz')
molDB[0].shape = 'linear'
# define the modsel with methods
model = ml.models.methods(method='B3LYP/6-31G*', program='PySCF')
for mol in molDB:
# run thermochemistry calculations
ml.thermochemistry(model=model, molecule=mol, program='ASE')
fmt = ' %-41s: %15.5f Hartree'
print(fmt % ('ZPE-exclusive internal energy at 0 K', mol.energy))
print(fmt % ('Zero-point vibrational energy', mol.ZPE))
print(fmt % ('Internal energy at 0 K', mol.U0))
print(fmt % ('Enthalpy at 298 K', mol.H))
print(fmt % ('Gibbs free energy at 298 K', mol.G))
Note
You can compare the accuracy and time consuming of different methods.