4. Universal machine learning models
4.1. Slides
4.2. What are the universal ML models?
Before we saw that ML models can be very fast and accurate at the same time. However, we needed to use (and train) the model specifically for the H2 example. This is not ideal and it would be nice to have ML models that can be applied to different systems of interest without needing specific training. Such models do exist and they are referred in literature as universal ML models. Example of the universal ML-based model is the UAIQM and its earliest incarnation AIQM1:
[Figure credit: reproduced from Pavlo O. Dral. AI in computational chemistry through the lens of a decade-long journey. Chem. Commun. 2024, 60, 3240–3258. DOI: 10.1039/D4CC00010B under the CC-BY license.]
AIQM1 is approaching the accuracy of CCSD(T)/CBS method and it is an ML-improved semi-empirical QM method, i.e., it is slower than a pure ML model. An example of a pure ML model is ANI-1ccx which also targets the accuracy of the CCSD(T)/CBS method (and, in fact, was published earlier than AIQM1). If both AIQM1 and ANI-1ccx methods are approaching the same CCSD(T)/CBS, then why bother using slower AIQM1? Let’s run an experiment:
Example ML-uni-1.
Optimize geometries of H2 and H2O with ANI-1ccx and AIQM1 (you are welcome to try UAIQM too!).
Check:
Method you used.
Bond lengths (and angle for H2O) you obtained.
How much time it took to complete the calculations.
Example of an input file:
geomopt # Request geometry optimization
AIQM1 # or ANI-1ccx - no model file is needed!
XYZfile='
2
H 0.0000000000000 0.0000000000000 0.0000000000000
H 0.0000000000000 0.0000000000000 0.7000000000000
'
optXYZ=h2_opt_AIQM1.xyz # optimized geometry output (this line is optional, by default it will be in optgeoms.xyz)
After you run these simulations, you should have obtained the following result (compared to other approaches):
One thing that you have noticed is that calculations with both methods are indeed very fast, much faster than with CCSD(T) methods you tried in your homework or task in Lecture 1. The optimized geometry of water molecule is also very close to experimental structure, much closer than with a slower B3LYP/6-31G* method. ANI-1ccx is slightly worse but still provides very good agreement with experiment for H2O.
Note
AIQM1 requires the semi-empirical program. On Aitomistic Lab@XMU, only eligible users can use the required semi-empirical program MNDO.
There is no such problem with other AIQM methods, so you can, e.g., try AIQM2, but we did not fit AIQM2 to accurately represent heats of formation (it is still good for many common reaction enthalpies though).
However, you might be surprised how bad ANI-1ccx is for H2. AIQM1 is also far from CCSD(T) quality. What happened?
AIQM1 and ANI-1ccx are both approaching CCSD(T)/CBS but AIQM1 ≠ ANI-1ccx ≠ CCSD(T)/CBS! They are ML-based models and trained on data which might or might not contain enough training points for H2. How the model was trained is also important. In general, you observe that AIQM1 is more robust than ANI-1ccx for the price of being slower. We will explore why AIQM1 is better in the future lectures.
The important lesson to learn is that any ML model might have limitations (like QM models too). Hence, always pay attention and use your expertise to check the results, i.e., by following recommendations:
compare with other high-level QM calculations (which can be unavailable or prohibitively expensive),
compare to experiment (may be not reliable either!),
check uncertainty (easiest and often works, see the second part or this lecture).
Those were just few molecules we checked. Overall, the quality of AIQM1 is excellent as evidenced by extensive benchmarks on many different molecules which show that AIQM1 optimized geometries are close to experiment and definetely better than many popular DFT methods which might even be qualitatively wrong (see example below for C18 carbon):
[Figure credit: reproduced from Pavlo O. Dral. AI in computational chemistry through the lens of a decade-long journey. Chem. Commun. 2024, 60, 3240–3258. DOI: 10.1039/D4CC00010B under the CC-BY license.]
Even more, AIQM1 was used to revise experimental geometries obtained via X-ray diffraction (see figure above).
You can safely use AIQM1 for quite large systems, where it is faster and more accurate than DFT, i.e., no reason to spend more time with DFT to get worse result!
[Figure credit: reproduced from Pavlo O. Dral, et al. 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 under CC-BY 4.0 license.]
You should be, however, aware of the current limitations of the universal models in general. See up-to-date list. The models available in MLatom on July 3, 2024, are overviewed below (partly based on exerpt from the MLatom 3 paper):
Our best offering: the UAIQM models.
Universal potentials ANI-1ccx, ANI-1x, ANI-2x, ANI-1x-D4, ANI-2x-D4, AIMnet-2, ANI-1xnr.
ANI-1ccx is the most accurate and approaches gold-standard CCSD(T) accuracy. Other methods approach the density functional theory (DFT) level.
ANI-1ccx and ANI-1x are limited to CHNO elements, while ANI-2x can be used for CHNOFClS elements.
We allow the user to use D4 dispersion-corrected universal ANI potentials that might be useful for noncovalent complexes.
Hybrid QM/ML methods AIQM2, AIQM1, AIQM1@DFT, AIQM1@DFT*, and DM21 are more transferable and accurate than pretrained ML models but slower (the speed of semiempirical QM methods, which are still much faster than DFT).
AIQM1 is approaching gold-standard CCSD(T) accuracy, while AIQM1@DFT and AIQM1@DFT* target the DFT accuracy for neutral, closed-shell molecules in their ground state. All these methods are limited to the CHNO elements. AIQM1 and AIQM1@DFT include explicit D4 dispersion corrections for the ωB97X functional, while AIQM1@DFT* does not.
These methods can also be used to calculate charged species, radicals, excited states, and other QM properties such as dipole moments, charges, oscillator strengths, and nonadiabatic couplings.
AIQM2 has improved accuracy over AIQM1, particularly, for reaction energies and barriers.
4.3. Frequencies and thermochemistry
After geometry optimization, the next step is to do frequency calculations (typically within rigid-rotor harmonic oscillator approximation) which achieves several goals: 1) obtain zero-point vibrational energy (molecules are always in motion and this should be taken into account!), 2) other thermochemical properties (enthalpies, Gibbs free energies, heats of formation, etc.), 3) perform normal mode analysis to check the nature of vibrations and their frequencies. These calculations also allow you to check whether your optimized geometries correspond to a true minima on potential energy surface (PES) - the minima have all real frequencies while you will see imaginary (shown as negative) frequencies for non-minima.
In principle, AIQM1 and ANI-1ccx are able to deliver very accurate thermochemical properties, particularly heats of formation with close-to chemical accuracy (errors not much larger than 1 kcal/mol). See it for yourself:
Example ML-uni-2.
Calculate frequencies in H2, H2O, and vinylacetylene molecules with ANI-1ccx.
Analyze your calculation results:
Are your optimized geometries true minima?
What are the heats of formation?
Is any of the calculations uncertain?
Example of an input file:
freq # Request frequency and thermochemistry calculations
ANI-1ccx
XYZfile=h2_opt_ani1ccx.xyz # file with optimized geometry
Note
First, you need to optimize geometries as described in the previous tasks. Frequency calculations must be performed at the same level as geometry optimizations.
Initial geometry for vinylacetylene:
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
Note
Frequency calculations require the evaluation of the second-order energy derivatives. Many methods do not have analytical derivatives, while numerical derivatives are often noisy and may lead to high inaccuracies and lots of negative frequencies which are artificial. This is the case, e.g., if AIQM1 is used on the cloud where the default QM program is Sparrow which has no analytical derivatives for AIQM1. Check in the MLatom program output what program is used - Sparrow (not reliable yet) or MNDO (can be trusted). Calculations with the semiempirical MNDO program are possible but the reader has to install it locally as described in MLatom’s documentation.
Creation of the initial XYZ geometry of vinylacetylene takes time. We do provide you with the initial guess, but what is a faster way to get XYZ geometry for a new molecule? I have asked ChatGPT in 2023 to generate XYZ coordinates or an input file for Gaussian program and it generated some gibberish:
That might improve with time though.
Note
Indeed, above statement is already outdated, now with our Aitomia on Aitomistic Lab@XMU, you can retrieve xyz files by providing names to the molecules. And visualize them in the web browser too! I did not have time to update it this here.
One more trick is to look up the coordinates in the internet as many molecules were already either calculated or have known experimental structures. Indeed, for vinylacetylene we can find XYZ structure optimized at a quite decent CCSD(T)/6-31G* level of theory and made available on a very useful NIST website (check it out – it is a treasure trove for any chemist!):
Going back to our exercise, you can find that for water molecule, it has no negative frequencies, i.e., it is a true minimum. Its heat of formation is predicted to be -58.46 kcal/mol at AIQM1 which is quite close to the experimental value of -57.8 kcal/mol (in gas phase! [J.B.Pedley, R.D.Naylor, and S.P.Kirby, “Thermochemical Data of Organic Compounds”, 2nd ed., Chapman and Hall, London, 1986]; the error is only 0.7 kcal/mol). And there is no warning about uncertainty, see below. The output file looks like this:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
Rotational symmetry number: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 1734.1034 1.0833 1.9194
2 3809.1282 1.0819 9.2489
3 3829.5496 1.0449 9.0285
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NNs : 0.00028358 Hartree 0.17795 kcal/mol
Energy : -76.38362520 Hartree
ZPE-exclusive internal energy at 0 K: -76.38363 Hartree
Zero-point vibrational energy : 0.02135 Hartree
Internal energy at 0 K: -76.36227 Hartree
Enthalpy at 298 K: -76.35849 Hartree
Gibbs free energy at 298 K: -76.38056 Hartree
Atomization enthalpy at 0 K: 0.35083 Hartree 220.14828 kcal/mol
ZPE-exclusive atomization energy at 0 K: 0.37218 Hartree 233.54734 kcal/mol
Heat of formation at 298 K: -0.09316 Hartree -58.45749 kcal/mol
==============================================================================
Just for your reference, it is -58.33 kcal/mol at AIQM1 (the error of only 0.5 kcal/mol), see the output file:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 1654.6962 1.0829 1.7470
2 3850.4309 1.0449 9.1277
3 3931.4391 1.0817 9.8508
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NN contribution : 0.00009088 Hartree 0.05703 kcal/mol
NN contribution : -0.00110508 Hartree
Sum of atomic self energies : -63.04934947 Hartree
ODM2* contribution : -13.33326549 Hartree
D4 contribution : -0.00002289 Hartree
Total energy : -76.38374293 Hartree
ZPE-exclusive internal energy at 0 K: -76.38374 Hartree
Zero-point vibrational energy : 0.02150 Hartree
Internal energy at 0 K: -76.36225 Hartree
Enthalpy at 298 K: -76.35847 Hartree
Gibbs free energy at 298 K: -76.37988 Hartree
Atomization enthalpy at 0 K: 0.35062 Hartree 220.01876 kcal/mol
ZPE-exclusive atomization energy at 0 K: 0.37212 Hartree 233.50896 kcal/mol
Heat of formation at 298 K: -0.09295 Hartree -58.32770 kcal/mol
==============================================================================
Wall-clock time: 4.54 s (0.08 min, 0.00 hours)
MLatom terminated on 25.04.2024 at 11:15:15
==============================================================================
4.4. Uncertainty quantification
For hydrogen molecule, we see that it is also a true minimum, its heat of formation is -23.09 kcal/mol. And as a chemist, you must know what the experimental value is and realize that the error is too large. However, you will notice that the output file contains a warning message * Warning * Heat of formation have high uncertainty!:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
Rotational symmetry number: 1
This is a linear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 2866.6581 1.0080 4.8805
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NNs : 0.00705607 Hartree 4.42775 kcal/mol
Energy : -1.20973475 Hartree
ZPE-exclusive internal energy at 0 K: -1.20973 Hartree
Zero-point vibrational energy : 0.00653 Hartree
Internal energy at 0 K: -1.20320 Hartree
Enthalpy at 298 K: -1.19990 Hartree
Gibbs free energy at 298 K: -1.21434 Hartree
Atomization enthalpy at 0 K: 0.20144 Hartree 126.40693 kcal/mol
ZPE-exclusive atomization energy at 0 K: 0.20797 Hartree 130.50502 kcal/mol
Heat of formation at 298 K: -0.03680 Hartree -23.09143 kcal/mol
* Warning * Heat of formation have high uncertainty!
==============================================================================
For your reference, the heat of formation at AIQM1 is not very accurate either: it is -2.86 kcal/mol and uncertain too:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
This is a linear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 3500.7003 1.0078 7.2769
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NN contribution : 0.00891919 Hartree 5.59688 kcal/mol
NN contribution : -0.00210675 Hartree
Sum of atomic self energies : -0.08587317 Hartree
ODM2* contribution : -1.09094346 Hartree
D4 contribution : -0.00000889 Hartree
Total energy : -1.17893227 Hartree
ZPE-exclusive internal energy at 0 K: -1.17893 Hartree
Zero-point vibrational energy : 0.00797 Hartree
Internal energy at 0 K: -1.17096 Hartree
Enthalpy at 298 K: -1.16765 Hartree
Gibbs free energy at 298 K: -1.18240 Hartree
Atomization enthalpy at 0 K: 0.16920 Hartree 106.17224 kcal/mol
ZPE-exclusive atomization energy at 0 K: 0.17717 Hartree 111.17663 kcal/mol
Heat of formation at 298 K: -0.00455 Hartree -2.85652 kcal/mol
* Warning * Heat of formation have high uncertainty!
==============================================================================
Wall-clock time: 4.49 s (0.07 min, 0.00 hours)
MLatom terminated on 25.04.2024 at 11:26:42
==============================================================================
Indeed, the calculations have high uncertainty as judged by the standard deviation of neural network (NN) contribution, it is quite large: 4.43 and 5.60 kcal/mol at ANI-1ccx and AIQM1, respectively, compared to only 0.18 and 0.06 kcal/mol for a water molecule. AIQM1 and ANI-1ccx both contain many (eight) NNs whose predictions are averaged. These NNs learned the same property and, hence, if they start to deviate from each other it means they are likely to stray away too far from the region they were trained on. The deviation between NNs is often used in ML to judge how bad the predictions are (e.g., used in query by committee and active learning, see next lectures). The question is how large deviation is still acceptable and hence we calibrated them on a standard data set and MLatom flags when deviation exceeds the calibrated thresholds.
This is consistent with our previous observation that ANI-1ccx and AIQM1 are not particularly good for the H2 molecule, i.e., the optimized bond lengths were quite off from the experimental value.
In case of vinylacetylene, we observe that it is also true minimum, its heat of formation is 68.53 kcal/mol, and it is not uncertain:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
Rotational symmetry number: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 234.3120 3.5401 0.1145
2 305.3683 2.6251 0.1442
3 542.0444 1.4458 0.2503
4 546.4789 1.2326 0.2169
5 559.3526 2.3867 0.4400
6 588.7133 1.5581 0.3182
7 900.8130 1.3032 0.6230
8 966.6489 2.6463 1.4569
9 991.1206 1.2769 0.7390
10 1163.7123 1.4019 1.1185
11 1332.3004 1.4117 1.4764
12 1499.1878 1.2277 1.6258
13 1750.8294 4.0546 7.3229
14 2206.4408 5.7716 16.5550
15 3205.7596 1.0906 6.6035
16 3262.3771 1.0654 6.6806
17 3342.5988 1.1208 7.3781
18 3388.8071 1.1735 7.9404
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NNs : 0.00057841 Hartree 0.36296 kcal/mol
Energy : -154.53147078 Hartree
ZPE-exclusive internal energy at 0 K: -154.53147 Hartree
Zero-point vibrational energy : 0.06102 Hartree
Internal energy at 0 K: -154.47045 Hartree
Enthalpy at 298 K: -154.46484 Hartree
Gibbs free energy at 298 K: -154.49658 Hartree
Atomization enthalpy at 0 K: 1.29896 Hartree 815.10992 kcal/mol
ZPE-exclusive atomization energy at 0 K: 1.35999 Hartree 853.40367 kcal/mol
Heat of formation at 298 K: 0.10920 Hartree 68.52699 kcal/mol
==============================================================================
At AIQM1, the heat of formation is 69.04 kcal/mol and is also confident:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 203.2675 3.0894 0.0752
2 293.6230 2.4319 0.1235
3 537.1646 3.5850 0.6095
4 671.2905 1.5015 0.3987
5 686.5656 1.2252 0.3403
6 720.9366 1.4999 0.4593
7 879.5200 2.4730 1.1271
8 931.0665 1.3842 0.7070
9 996.7269 1.0846 0.6348
10 1099.9357 1.6593 1.1828
11 1313.2985 1.3102 1.3315
12 1431.8700 1.2556 1.5167
13 1645.4021 3.8620 6.1603
14 2162.6121 6.3527 17.5052
15 3135.1808 1.0791 6.2496
16 3174.0848 1.0667 6.3318
17 3260.7611 1.1187 7.0083
18 3456.1498 1.1511 8.1013
==============================================================================
Thermochemistry for molecule 1
==============================================================================
Standard deviation of NN contribution : 0.00008761 Hartree 0.05498 kcal/mol
NN contribution : 0.00715111 Hartree
Sum of atomic self energies : -133.90358074 Hartree
ODM2* contribution : -20.63402319 Hartree
D4 contribution : -0.00056531 Hartree
Total energy : -154.53101813 Hartree
ZPE-exclusive internal energy at 0 K: -154.53102 Hartree
Zero-point vibrational energy : 0.06060 Hartree
Internal energy at 0 K: -154.47042 Hartree
Enthalpy at 298 K: -154.46492 Hartree
Gibbs free energy at 298 K: -154.49656 Hartree
Atomization enthalpy at 0 K: 1.29803 Hartree 814.52617 kcal/mol
ZPE-exclusive atomization energy at 0 K: 1.35863 Hartree 852.55199 kcal/mol
Heat of formation at 298 K: 0.11002 Hartree 69.04119 kcal/mol
==============================================================================
Wall-clock time: 6.03 s (0.10 min, 0.00 hours)
MLatom terminated on 25.04.2024 at 11:31:47
==============================================================================
Nevertheless, the experimental value used in one of the databases for fitting semi-empirical methods, was 73 kcal/mol. This value deviates from both ANI-1ccx and AIQM1’s predictions by 4 kcal/mol which is too much for confident calculations. Hence, we double checked (and you can too) by calculating the value with another QM method (G4 as implemented in the Gaussian program) and it also gave heat of formation of 69.1 kcal/mol, i.e., basically the same as AIQM1. Hence, I suspected that something might be wrong with that experimental value and, indeed, I found in the NIST database a value of 70.4 kcal/mol which is much closer to the ANI-1ccx and AIQM1 predictions. This is another example (among many similar) where accurate AIQM1 calculations helped to revise the experimental values.
Now let’s see what we can get if we use DFT.
Example ML-uni-3.
Calculate heat of formation of H2O with B3LYP/6-31G*.
Report:
What is the heat of formation?
Is this calculation uncertain?
Example of an input file (do not forget to optimize geometry too!):
freq # Request frequency and thermochemistry calculations
method=B3LYPG/6-31G* # Using PySCF instead of Gaussian
qmprog=pyscf # Use PySCF instead of Gaussian
freqprog=pyscf
XYZfile=h2o_opt_B3LYP.xyz # file with optimized geometry
As you can see in the output file, the heat of formation at B3LYP is -42.30 kcal/mol which is very wrong compared to the experimental -57.8 kcal/mol!
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
Rotational symmetry number: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
1 1710.8681 1.0830 1.8677
2 3720.9628 1.0452 8.5264
3 3843.9202 1.0813 9.4132
==============================================================================
Thermochemistry for molecule 1
==============================================================================
ZPE-exclusive internal energy at 0 K: -76.36989 Hartree
Zero-point vibrational energy : 0.02113 Hartree
Internal energy at 0 K: -76.34876 Hartree
Enthalpy at 298 K: -76.34498 Hartree
Gibbs free energy at 298 K: -76.36708 Hartree
Atomization enthalpy at 0 K: 0.32508 Hartree 203.98967 kcal/mol
ZPE-exclusive atomization energy at 0 K: 0.34621 Hartree 217.25003 kcal/mol
Heat of formation at 298 K: -0.06741 Hartree -42.29876 kcal/mol
==============================================================================
Wall-clock time: 3.43 s (0.06 min, 0.00 hours)
It is indeed well known that DFT methods struggle with getting heats of formation properly. We will discuss this issue in the next lecture.
Another important observation is that DFT (or any other QM method) never provides you any uncertainty of its prediction! This is a big advantage of ML methods that they do provide uncertainty quantification. That’s why when many people complain that they are not sure whether ML is good, they forget that they can never be sure whether their DFT method is good. In ML we at least get an uncertainty of predictions and can improve the ML model later with more data.
4.5. Qualitative problems with PES shape
One of the problems with ML and universal ML is that they often do not describe correctly the dissociation curves of molecules. See it for yourself:
Example ML-uni-4.
Optimize and calculate frequencies of H2 with ANI-1ccx but starting not from initial bond length of 0.7 A as in the previous task but with 0.8 A.
Questions:
Is the geometry you have obtained the same as in the previous task?
Is your optimized geometry true minimum?
What is its heats of formation?
Is this calculation confident?
4.6. Jupyter notebook for additional interactive practice
import mlatom as ml
Universal ANI-1ccx MLP¶
ANI-1ccx is targeting CCSD(T)/CBS - let's see how well it works.
This is just for demonstration - we recommend to use UAIQM in real applications (http://mlatom.com/docs/tutorial_uaiqm.html)
ani1ccx = ml.models.methods(method='ANI-1ccx')
Example 1: Heat of formation¶
Let's calculate the heat of formation of vinylacetylene
# initial xyz geometry
va_init = ml.molecule.from_xyz_string('''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
''')
# step 1 - optimize the geometry
geomopt = ml.simulations.optimize_geometry(model=ani1ccx, initial_molecule=va_init)
va_opt = geomopt.optimized_molecule
va_opt.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
# step 2 - run frequency and thermochemistry calculations
ml.simulations.thermochemistry(model=ani1ccx, molecule=va_opt, program='pyscf')
<mlatom.simulations.thermochemistry at 0x2ad6bfd14310>
# Check vibration analysis
print("Mode Frequencies Reduced masses Force Constants")
print(" (cm^-1) (AMU) (mDyne/A)")
for ii in range(len(va_opt.frequencies)):
print("%d %13.4f %13.4f %13.4f"%(ii,va_opt.frequencies[ii],va_opt.reduced_masses[ii],va_opt.force_constants[ii]))
# Check thermochemistry results
print('')
print(f"Zero-point vibrational energy: {va_opt.ZPE} Hartree")
print(f"Enthalpy at 298 K: {va_opt.H} Hartree")
print(f"Gibbs Free energy at 298 K: {va_opt.G} Hartree")
print(f"Heat of formation at 298 K: {va_opt.DeltaHf298 * ml.constants.Hartree2kcalpermol} kcal/mol")
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
0 235.1708 3.5411 0.1154
1 305.7176 2.6230 0.1444
2 542.5744 1.4393 0.2496
3 547.0089 1.2320 0.2172
4 560.1454 2.4045 0.4445
5 588.3991 1.5591 0.3180
6 901.1611 1.3035 0.6237
7 968.7884 2.6410 1.4604
8 991.7167 1.2770 0.7400
9 1164.2582 1.4018 1.1195
10 1332.4180 1.4170 1.4821
11 1500.0635 1.2295 1.6300
12 1747.1288 4.0202 7.2302
13 2210.5392 5.6998 16.4098
14 3205.7453 1.0907 6.6043
15 3268.7570 1.0644 6.7009
16 3346.6297 1.1204 7.3933
17 3392.5651 1.1767 7.9797
Zero-point vibrational energy: 0.06107491155855896 Hartree
Enthalpy at 298 K: -154.46479363485236 Hartree
Gibbs Free energy at 298 K: -154.49652545131622 Hartree
Heat of formation at 298 K: 68.55605492834317 kcal/mol
# check the uncertainty (UQ threshold is 1.68 kcal/mol)
va_opt.ani1ccx.energy_standard_deviation * ml.constants.Hartree2kcalpermol
0.36496787659098034
# let's view the first vibrational normal mode
va_opt.view(normal_mode=0)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Example 2: MD¶
Let's do MD of the boring hydrogen molecule H2
# use random initial conditions
mol = ml.data.molecule.from_xyz_string('''2
H 0.0000000000000 0.0000000000000 0.0000000000000
H 0.0000000000000 0.0000000000000 0.8000000000000
''')
init_cond_db = ml.generate_initial_conditions(molecule=mol,
generation_method='random',
initial_temperature=300, # set initial temperature in K,
number_of_initial_conditions=1 # we just want to run a single trajectory
)
init_mol = init_cond_db[0]
# Initializing thermostat
nose_hoover = ml.md.Nose_Hoover_thermostat(temperature=300, molecule=init_mol)
# Run molecular dynamics
dyn = ml.md(model=ani1ccx,
molecule_with_initial_conditions=init_mol,
thermostat=nose_hoover,
ensemble='NVT',
time_step=0.5, # fs
maximum_propagation_time=100.0 # fs
)
# Dump trajectory
traj = dyn.molecular_trajectory
traj.dump(filename='traj', format='plain_text')
traj.dump(filename='traj.h5', format='h5md')
print(f"Number of steps in the trajectory: {len(traj.steps)}")
Number of steps in the trajectory: 201
traj.view()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
for step in traj.steps:
print(f'{step.step:3d} {step.time:8.2f} {step.molecule.internuclear_distance(0, 1):8.2f} {step.molecule.atoms[1].xyz_coordinates[2] - step.molecule.atoms[0].xyz_coordinates[2]:8.2f}')
0 0.00 0.80 0.80 1 0.50 0.77 0.77 2 1.00 0.71 0.71 3 1.50 0.61 0.61 4 2.00 0.51 0.51 5 2.50 0.43 0.43 6 3.00 0.37 0.37 7 3.50 0.34 0.34 8 4.00 0.31 0.31 9 4.50 0.28 0.28 10 5.00 0.26 0.26 11 5.50 0.25 0.25 12 6.00 0.23 0.23 13 6.50 0.22 0.22 14 7.00 0.21 0.21 15 7.50 0.20 0.20 16 8.00 0.19 0.19 17 8.50 0.18 0.18 18 9.00 0.17 0.17 19 9.50 0.17 0.17 20 10.00 0.16 0.16 21 10.50 0.15 0.15 22 11.00 0.15 0.15 23 11.50 0.14 0.14 24 12.00 0.14 0.14 25 12.50 0.14 0.14 26 13.00 0.13 0.13 27 13.50 0.13 0.13 28 14.00 0.12 0.12 29 14.50 0.12 0.12 30 15.00 0.12 0.12 31 15.50 0.12 0.12 32 16.00 0.11 0.11 33 16.50 0.11 0.11 34 17.00 0.11 0.11 35 17.50 0.11 0.11 36 18.00 0.10 0.10 37 18.50 0.10 0.10 38 19.00 0.10 0.10 39 19.50 0.10 0.10 40 20.00 0.10 0.10 41 20.50 0.09 0.09 42 21.00 0.09 0.09 43 21.50 0.09 0.09 44 22.00 0.09 0.09 45 22.50 0.09 0.09 46 23.00 0.08 0.08 47 23.50 0.08 0.08 48 24.00 0.08 0.08 49 24.50 0.08 0.08 50 25.00 0.08 0.08 51 25.50 0.08 0.08 52 26.00 0.07 0.07 53 26.50 0.07 0.07 54 27.00 0.07 0.07 55 27.50 0.07 0.07 56 28.00 0.07 0.07 57 28.50 0.06 0.06 58 29.00 0.06 0.06 59 29.50 0.06 0.06 60 30.00 0.06 0.06 61 30.50 0.06 0.06 62 31.00 0.05 0.05 63 31.50 0.05 0.05 64 32.00 0.05 0.05 65 32.50 0.04 0.04 66 33.00 0.04 0.04 67 33.50 0.04 0.04 68 34.00 0.03 0.03 69 34.50 0.03 0.03 70 35.00 0.03 0.03 71 35.50 0.02 0.02 72 36.00 0.02 0.02 73 36.50 0.01 0.01 74 37.00 0.01 0.01 75 37.50 0.00 0.00 76 38.00 0.01 -0.01 77 38.50 0.01 -0.01 78 39.00 0.02 -0.02 79 39.50 0.03 -0.03 80 40.00 0.04 -0.04 81 40.50 0.05 -0.05 82 41.00 0.06 -0.06 83 41.50 0.07 -0.07 84 42.00 0.09 -0.09 85 42.50 0.10 -0.10 86 43.00 0.12 -0.12 87 43.50 0.13 -0.13 88 44.00 0.15 -0.15 89 44.50 0.17 -0.17 90 45.00 0.19 -0.19 91 45.50 0.21 -0.21 92 46.00 0.23 -0.23 93 46.50 0.25 -0.25 94 47.00 0.27 -0.27 95 47.50 0.29 -0.29 96 48.00 0.31 -0.31 97 48.50 0.34 -0.34 98 49.00 0.36 -0.36 99 49.50 0.39 -0.39 100 50.00 0.41 -0.41 101 50.50 0.44 -0.44 102 51.00 0.46 -0.46 103 51.50 0.48 -0.48 104 52.00 0.49 -0.49 105 52.50 0.50 -0.50 106 53.00 0.49 -0.49 107 53.50 0.48 -0.48 108 54.00 0.47 -0.47 109 54.50 0.45 -0.45 110 55.00 0.43 -0.43 111 55.50 0.42 -0.42 112 56.00 0.40 -0.40 113 56.50 0.39 -0.39 114 57.00 0.39 -0.39 115 57.50 0.38 -0.38 116 58.00 0.38 -0.38 117 58.50 0.37 -0.37 118 59.00 0.37 -0.37 119 59.50 0.37 -0.37 120 60.00 0.38 -0.38 121 60.50 0.38 -0.38 122 61.00 0.39 -0.39 123 61.50 0.39 -0.39 124 62.00 0.40 -0.40 125 62.50 0.41 -0.41 126 63.00 0.42 -0.42 127 63.50 0.43 -0.43 128 64.00 0.44 -0.44 129 64.50 0.45 -0.45 130 65.00 0.46 -0.46 131 65.50 0.47 -0.47 132 66.00 0.47 -0.47 133 66.50 0.47 -0.47 134 67.00 0.46 -0.46 135 67.50 0.46 -0.46 136 68.00 0.45 -0.45 137 68.50 0.44 -0.44 138 69.00 0.43 -0.43 139 69.50 0.42 -0.42 140 70.00 0.41 -0.41 141 70.50 0.41 -0.41 142 71.00 0.40 -0.40 143 71.50 0.40 -0.40 144 72.00 0.40 -0.40 145 72.50 0.40 -0.40 146 73.00 0.40 -0.40 147 73.50 0.40 -0.40 148 74.00 0.41 -0.41 149 74.50 0.41 -0.41 150 75.00 0.42 -0.42 151 75.50 0.43 -0.43 152 76.00 0.44 -0.44 153 76.50 0.45 -0.45 154 77.00 0.45 -0.45 155 77.50 0.46 -0.46 156 78.00 0.46 -0.46 157 78.50 0.47 -0.47 158 79.00 0.46 -0.46 159 79.50 0.46 -0.46 160 80.00 0.45 -0.45 161 80.50 0.44 -0.44 162 81.00 0.43 -0.43 163 81.50 0.42 -0.42 164 82.00 0.41 -0.41 165 82.50 0.40 -0.40 166 83.00 0.40 -0.40 167 83.50 0.39 -0.39 168 84.00 0.38 -0.38 169 84.50 0.38 -0.38 170 85.00 0.38 -0.38 171 85.50 0.38 -0.38 172 86.00 0.38 -0.38 173 86.50 0.38 -0.38 174 87.00 0.39 -0.39 175 87.50 0.39 -0.39 176 88.00 0.40 -0.40 177 88.50 0.41 -0.41 178 89.00 0.43 -0.43 179 89.50 0.44 -0.44 180 90.00 0.46 -0.46 181 90.50 0.47 -0.47 182 91.00 0.49 -0.49 183 91.50 0.49 -0.49 184 92.00 0.49 -0.49 185 92.50 0.48 -0.48 186 93.00 0.47 -0.47 187 93.50 0.45 -0.45 188 94.00 0.43 -0.43 189 94.50 0.41 -0.41 190 95.00 0.38 -0.38 191 95.50 0.36 -0.36 192 96.00 0.35 -0.35 193 96.50 0.33 -0.33 194 97.00 0.31 -0.31 195 97.50 0.30 -0.30 196 98.00 0.28 -0.28 197 98.50 0.27 -0.27 198 99.00 0.26 -0.26 199 99.50 0.25 -0.25 200 100.00 0.23 -0.23