{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "64aa60ce-422a-4237-b4f6-68b768728e7e", "metadata": {}, "outputs": [], "source": [ "import mlatom as ml" ] }, { "cell_type": "markdown", "id": "47237837-7d5f-4af5-9faa-fa8413fb4405", "metadata": {}, "source": [ "# Universal ANI-1ccx MLP\n", "ANI-1ccx is targeting CCSD(T)/CBS - let's see how well it works.\n", "\n", "This is just for demonstration - we recommend to use UAIQM in real applications (http://mlatom.com/docs/tutorial_uaiqm.html)" ] }, { "cell_type": "code", "execution_count": 2, "id": "e2aee6d6-9830-442a-841a-0fa4719ab6f1", "metadata": {}, "outputs": [], "source": [ "ani1ccx = ml.models.methods(method='ANI-1ccx')" ] }, { "cell_type": "markdown", "id": "80951cf2-474b-4be1-ae50-c8cfaac9600e", "metadata": {}, "source": [ "## Example 1: Heat of formation\n", "Let's calculate the heat of formation of vinylacetylene" ] }, { "cell_type": "code", "execution_count": 3, "id": "ea3d8332-2b29-484e-974b-0c2062c3f127", "metadata": {}, "outputs": [], "source": [ "# initial xyz geometry\n", "va_init = ml.molecule.from_xyz_string('''8\n", "\n", "C 0.00000 0.00000 0.00000\n", "C 1.33718 0.00000 0.00000\n", "C 2.04862 1.23465 0.00000\n", "C 2.65517 2.30065 0.00000\n", "H 3.17326 3.20500 0.00000\n", "H -0.56880 0.91919 0.00000\n", "H -0.56652 -0.92292 -0.00000\n", "H 1.95564 -0.90366 0.00000\n", "''')" ] }, { "cell_type": "code", "execution_count": 4, "id": "70ee2eb4-3364-4c3b-b052-7095953edbb0", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# step 1 - optimize the geometry\n", "geomopt = ml.simulations.optimize_geometry(model=ani1ccx, initial_molecule=va_init)\n", "va_opt = geomopt.optimized_molecule\n", "va_opt.view()" ] }, { "cell_type": "code", "execution_count": 5, "id": "18c97e55-6d6f-44d7-8f48-3e37c023fca1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# step 2 - run frequency and thermochemistry calculations\n", "ml.simulations.thermochemistry(model=ani1ccx, molecule=va_opt, program='pyscf')" ] }, { "cell_type": "code", "execution_count": 6, "id": "fe83cc61-3f72-40c5-a237-6941635a42cd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mode Frequencies Reduced masses Force Constants\n", " (cm^-1) (AMU) (mDyne/A)\n", "0 235.1708 3.5411 0.1154\n", "1 305.7176 2.6230 0.1444\n", "2 542.5744 1.4393 0.2496\n", "3 547.0089 1.2320 0.2172\n", "4 560.1454 2.4045 0.4445\n", "5 588.3991 1.5591 0.3180\n", "6 901.1611 1.3035 0.6237\n", "7 968.7884 2.6410 1.4604\n", "8 991.7167 1.2770 0.7400\n", "9 1164.2582 1.4018 1.1195\n", "10 1332.4180 1.4170 1.4821\n", "11 1500.0635 1.2295 1.6300\n", "12 1747.1288 4.0202 7.2302\n", "13 2210.5392 5.6998 16.4098\n", "14 3205.7453 1.0907 6.6043\n", "15 3268.7570 1.0644 6.7009\n", "16 3346.6297 1.1204 7.3933\n", "17 3392.5651 1.1767 7.9797\n", "\n", "Zero-point vibrational energy: 0.06107491155855896 Hartree\n", "Enthalpy at 298 K: -154.46479363485236 Hartree\n", "Gibbs Free energy at 298 K: -154.49652545131622 Hartree\n", "Heat of formation at 298 K: 68.55605492834317 kcal/mol\n" ] } ], "source": [ "# Check vibration analysis\n", "print(\"Mode Frequencies Reduced masses Force Constants\")\n", "print(\" (cm^-1) (AMU) (mDyne/A)\")\n", "for ii in range(len(va_opt.frequencies)):\n", " print(\"%d %13.4f %13.4f %13.4f\"%(ii,va_opt.frequencies[ii],va_opt.reduced_masses[ii],va_opt.force_constants[ii]))\n", "\n", "# Check thermochemistry results\n", "print('')\n", "print(f\"Zero-point vibrational energy: {va_opt.ZPE} Hartree\")\n", "print(f\"Enthalpy at 298 K: {va_opt.H} Hartree\")\n", "print(f\"Gibbs Free energy at 298 K: {va_opt.G} Hartree\")\n", "print(f\"Heat of formation at 298 K: {va_opt.DeltaHf298 * ml.constants.Hartree2kcalpermol} kcal/mol\")" ] }, { "cell_type": "code", "execution_count": 7, "id": "fc856a19-f7d0-418a-8e69-8f2186ccf57e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.36496787659098034" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check the uncertainty (UQ threshold is 1.68 kcal/mol)\n", "va_opt.ani1ccx.energy_standard_deviation * ml.constants.Hartree2kcalpermol" ] }, { "cell_type": "code", "execution_count": 8, "id": "696f78d1-a894-46b5-8124-ef4d7f82cba2", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# let's view the first vibrational normal mode\n", "va_opt.view(normal_mode=0)" ] }, { "cell_type": "markdown", "id": "2ba423f9-d949-4843-a151-ccadb9847d0d", "metadata": {}, "source": [ "## Example 2: MD\n", "Let's do MD of the boring hydrogen molecule H2" ] }, { "cell_type": "code", "execution_count": 9, "id": "fedc636b-40ce-43c8-a0bc-0f88ed955166", "metadata": {}, "outputs": [], "source": [ "# use random initial conditions\n", "mol = ml.data.molecule.from_xyz_string('''2\n", "\n", "H 0.0000000000000 0.0000000000000 0.0000000000000\n", "H 0.0000000000000 0.0000000000000 0.8000000000000\n", "''')\n", "init_cond_db = ml.generate_initial_conditions(molecule=mol,\n", " generation_method='random',\n", " initial_temperature=300, # set initial temperature in K,\n", " number_of_initial_conditions=1 # we just want to run a single trajectory\n", " )\n", "init_mol = init_cond_db[0]" ] }, { "cell_type": "code", "execution_count": 10, "id": "4d1c5f28-f45c-4b00-9e35-032b6b6086b7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of steps in the trajectory: 201\n" ] } ], "source": [ "# Initializing thermostat\n", "nose_hoover = ml.md.Nose_Hoover_thermostat(temperature=300, molecule=init_mol)\n", "\n", "# Run molecular dynamics\n", "dyn = ml.md(model=ani1ccx,\n", " molecule_with_initial_conditions=init_mol,\n", " thermostat=nose_hoover,\n", " ensemble='NVT',\n", " time_step=0.5, # fs\n", " maximum_propagation_time=100.0 # fs\n", " )\n", "\n", "# Dump trajectory\n", "traj = dyn.molecular_trajectory\n", "traj.dump(filename='traj', format='plain_text')\n", "traj.dump(filename='traj.h5', format='h5md')\n", "\n", "print(f\"Number of steps in the trajectory: {len(traj.steps)}\")" ] }, { "cell_type": "code", "execution_count": 11, "id": "0aed6b6e-5550-4688-a29f-2775e874f85b", "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", "text/html": [ "
\n", "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "traj.view()" ] }, { "cell_type": "code", "execution_count": 12, "id": "a0c3cdad-18d3-491b-af14-fd4ca6b7209b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 0.00 0.80 0.80\n", " 1 0.50 0.77 0.77\n", " 2 1.00 0.71 0.71\n", " 3 1.50 0.61 0.61\n", " 4 2.00 0.51 0.51\n", " 5 2.50 0.43 0.43\n", " 6 3.00 0.37 0.37\n", " 7 3.50 0.34 0.34\n", " 8 4.00 0.31 0.31\n", " 9 4.50 0.28 0.28\n", " 10 5.00 0.26 0.26\n", " 11 5.50 0.25 0.25\n", " 12 6.00 0.23 0.23\n", " 13 6.50 0.22 0.22\n", " 14 7.00 0.21 0.21\n", " 15 7.50 0.20 0.20\n", " 16 8.00 0.19 0.19\n", " 17 8.50 0.18 0.18\n", " 18 9.00 0.17 0.17\n", " 19 9.50 0.17 0.17\n", " 20 10.00 0.16 0.16\n", " 21 10.50 0.15 0.15\n", " 22 11.00 0.15 0.15\n", " 23 11.50 0.14 0.14\n", " 24 12.00 0.14 0.14\n", " 25 12.50 0.14 0.14\n", " 26 13.00 0.13 0.13\n", " 27 13.50 0.13 0.13\n", " 28 14.00 0.12 0.12\n", " 29 14.50 0.12 0.12\n", " 30 15.00 0.12 0.12\n", " 31 15.50 0.12 0.12\n", " 32 16.00 0.11 0.11\n", " 33 16.50 0.11 0.11\n", " 34 17.00 0.11 0.11\n", " 35 17.50 0.11 0.11\n", " 36 18.00 0.10 0.10\n", " 37 18.50 0.10 0.10\n", " 38 19.00 0.10 0.10\n", " 39 19.50 0.10 0.10\n", " 40 20.00 0.10 0.10\n", " 41 20.50 0.09 0.09\n", " 42 21.00 0.09 0.09\n", " 43 21.50 0.09 0.09\n", " 44 22.00 0.09 0.09\n", " 45 22.50 0.09 0.09\n", " 46 23.00 0.08 0.08\n", " 47 23.50 0.08 0.08\n", " 48 24.00 0.08 0.08\n", " 49 24.50 0.08 0.08\n", " 50 25.00 0.08 0.08\n", " 51 25.50 0.08 0.08\n", " 52 26.00 0.07 0.07\n", " 53 26.50 0.07 0.07\n", " 54 27.00 0.07 0.07\n", " 55 27.50 0.07 0.07\n", " 56 28.00 0.07 0.07\n", " 57 28.50 0.06 0.06\n", " 58 29.00 0.06 0.06\n", " 59 29.50 0.06 0.06\n", " 60 30.00 0.06 0.06\n", " 61 30.50 0.06 0.06\n", " 62 31.00 0.05 0.05\n", " 63 31.50 0.05 0.05\n", " 64 32.00 0.05 0.05\n", " 65 32.50 0.04 0.04\n", " 66 33.00 0.04 0.04\n", " 67 33.50 0.04 0.04\n", " 68 34.00 0.03 0.03\n", " 69 34.50 0.03 0.03\n", " 70 35.00 0.03 0.03\n", " 71 35.50 0.02 0.02\n", " 72 36.00 0.02 0.02\n", " 73 36.50 0.01 0.01\n", " 74 37.00 0.01 0.01\n", " 75 37.50 0.00 0.00\n", " 76 38.00 0.01 -0.01\n", " 77 38.50 0.01 -0.01\n", " 78 39.00 0.02 -0.02\n", " 79 39.50 0.03 -0.03\n", " 80 40.00 0.04 -0.04\n", " 81 40.50 0.05 -0.05\n", " 82 41.00 0.06 -0.06\n", " 83 41.50 0.07 -0.07\n", " 84 42.00 0.09 -0.09\n", " 85 42.50 0.10 -0.10\n", " 86 43.00 0.12 -0.12\n", " 87 43.50 0.13 -0.13\n", " 88 44.00 0.15 -0.15\n", " 89 44.50 0.17 -0.17\n", " 90 45.00 0.19 -0.19\n", " 91 45.50 0.21 -0.21\n", " 92 46.00 0.23 -0.23\n", " 93 46.50 0.25 -0.25\n", " 94 47.00 0.27 -0.27\n", " 95 47.50 0.29 -0.29\n", " 96 48.00 0.31 -0.31\n", " 97 48.50 0.34 -0.34\n", " 98 49.00 0.36 -0.36\n", " 99 49.50 0.39 -0.39\n", "100 50.00 0.41 -0.41\n", "101 50.50 0.44 -0.44\n", "102 51.00 0.46 -0.46\n", "103 51.50 0.48 -0.48\n", "104 52.00 0.49 -0.49\n", "105 52.50 0.50 -0.50\n", "106 53.00 0.49 -0.49\n", "107 53.50 0.48 -0.48\n", "108 54.00 0.47 -0.47\n", "109 54.50 0.45 -0.45\n", "110 55.00 0.43 -0.43\n", "111 55.50 0.42 -0.42\n", "112 56.00 0.40 -0.40\n", "113 56.50 0.39 -0.39\n", "114 57.00 0.39 -0.39\n", "115 57.50 0.38 -0.38\n", "116 58.00 0.38 -0.38\n", "117 58.50 0.37 -0.37\n", "118 59.00 0.37 -0.37\n", "119 59.50 0.37 -0.37\n", "120 60.00 0.38 -0.38\n", "121 60.50 0.38 -0.38\n", "122 61.00 0.39 -0.39\n", "123 61.50 0.39 -0.39\n", "124 62.00 0.40 -0.40\n", "125 62.50 0.41 -0.41\n", "126 63.00 0.42 -0.42\n", "127 63.50 0.43 -0.43\n", "128 64.00 0.44 -0.44\n", "129 64.50 0.45 -0.45\n", "130 65.00 0.46 -0.46\n", "131 65.50 0.47 -0.47\n", "132 66.00 0.47 -0.47\n", "133 66.50 0.47 -0.47\n", "134 67.00 0.46 -0.46\n", "135 67.50 0.46 -0.46\n", "136 68.00 0.45 -0.45\n", "137 68.50 0.44 -0.44\n", "138 69.00 0.43 -0.43\n", "139 69.50 0.42 -0.42\n", "140 70.00 0.41 -0.41\n", "141 70.50 0.41 -0.41\n", "142 71.00 0.40 -0.40\n", "143 71.50 0.40 -0.40\n", "144 72.00 0.40 -0.40\n", "145 72.50 0.40 -0.40\n", "146 73.00 0.40 -0.40\n", "147 73.50 0.40 -0.40\n", "148 74.00 0.41 -0.41\n", "149 74.50 0.41 -0.41\n", "150 75.00 0.42 -0.42\n", "151 75.50 0.43 -0.43\n", "152 76.00 0.44 -0.44\n", "153 76.50 0.45 -0.45\n", "154 77.00 0.45 -0.45\n", "155 77.50 0.46 -0.46\n", "156 78.00 0.46 -0.46\n", "157 78.50 0.47 -0.47\n", "158 79.00 0.46 -0.46\n", "159 79.50 0.46 -0.46\n", "160 80.00 0.45 -0.45\n", "161 80.50 0.44 -0.44\n", "162 81.00 0.43 -0.43\n", "163 81.50 0.42 -0.42\n", "164 82.00 0.41 -0.41\n", "165 82.50 0.40 -0.40\n", "166 83.00 0.40 -0.40\n", "167 83.50 0.39 -0.39\n", "168 84.00 0.38 -0.38\n", "169 84.50 0.38 -0.38\n", "170 85.00 0.38 -0.38\n", "171 85.50 0.38 -0.38\n", "172 86.00 0.38 -0.38\n", "173 86.50 0.38 -0.38\n", "174 87.00 0.39 -0.39\n", "175 87.50 0.39 -0.39\n", "176 88.00 0.40 -0.40\n", "177 88.50 0.41 -0.41\n", "178 89.00 0.43 -0.43\n", "179 89.50 0.44 -0.44\n", "180 90.00 0.46 -0.46\n", "181 90.50 0.47 -0.47\n", "182 91.00 0.49 -0.49\n", "183 91.50 0.49 -0.49\n", "184 92.00 0.49 -0.49\n", "185 92.50 0.48 -0.48\n", "186 93.00 0.47 -0.47\n", "187 93.50 0.45 -0.45\n", "188 94.00 0.43 -0.43\n", "189 94.50 0.41 -0.41\n", "190 95.00 0.38 -0.38\n", "191 95.50 0.36 -0.36\n", "192 96.00 0.35 -0.35\n", "193 96.50 0.33 -0.33\n", "194 97.00 0.31 -0.31\n", "195 97.50 0.30 -0.30\n", "196 98.00 0.28 -0.28\n", "197 98.50 0.27 -0.27\n", "198 99.00 0.26 -0.26\n", "199 99.50 0.25 -0.25\n", "200 100.00 0.23 -0.23\n" ] } ], "source": [ "for step in traj.steps:\n", " 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}')" ] }, { "cell_type": "code", "execution_count": null, "id": "31afe0dd-186a-4e93-8efa-739ae0c660c5", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }