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. - :ref:`ANI-1ccx ` - :ref:`AIQM1 ` - :ref:`B3LYP ` Before optimizing, the initial geometries of hydrogen and vinyl acetylene (:download:`h2_vinylacetylene_init.xyz `) should be provided. .. code-block:: # 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 ~~~~~~~~ `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 (:download:`h2_vinylacetylene_opt_ANI1ccx.inp `). .. code-block:: # 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``. .. code-block:: mlatom h2_vinylacetylene_opt_ANI1ccx.inp There is the code of geometry optimization in Python. .. code-block:: 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: :download:`h2_vinylacetylene_freq_ANI1ccx.inp `. .. code-block:: # h2_vinylacetylene_freq_ANI1ccx.inp ANI-1ccx freq xyzfile=final_ANI1ccx.xyz optprog=ase ase.linear=1,0 Run it. .. code-block:: mlatom h2_vinylacetylene_freq_ANI1ccx.inp In Python, we need to provide the optimized geometry ``final_ANI1ccx.xyz``. .. code-block:: 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 ~~~~~ `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 :download:`h2_vinylacetylene_opt_AIQM1.inp ` in command line. .. code-block:: # 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``. .. code-block:: mlatom h2_vinylacetylene_opt_AIQM1.inp There is the code of geometry optimization in Python. .. code-block:: 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: :download:`h2_vinylacetylene_freq_AIQM1.inp `. .. code-block:: # h2_vinylacetylene_freq_AIQM1.inp AIQM1 freq xyzfile=final_AIQM1.xyz optprog=ase ase.linear=1,0 Run it. .. code-block:: mlatom h2_vinylacetylene_freq_AIQM1.inp In Python, we need to provide the optimized geometry ``final_AIQM1.xyz``. .. code-block:: 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 ~~~~~~ B3LYP is a common DFT method which is widely used. Now we use this method to optimize the geometry. Prepare the input file fistly: :download:`h2_vinylacetylene_opt_B3LYP.inp ` .. code-block:: # 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``. .. code-block:: mlatom h2_vinylacetylene_opt_B3LYP.inp There is the code of geometry optimization in Python. .. code-block:: 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: :download:`h2_vinylacetylene_freq_B3LYP.inp `. .. code-block:: # h2_vinylacetylene_freq_B3LYP.inp method=B3LYP/6-31G* QMprog=PySCF freq xyzfile=final_B3LYP.xyz optprog=ase ase.linear=1,0 Run it. .. code-block:: mlatom h2_vinylacetylene_freq_B3LYP.inp In Python, we need to provide the optimized geometry ``final_B3LYP.xyz``. .. code-block:: 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.