diff --git a/doc/tutorials/mlip/mlip.ipynb b/doc/tutorials/mlip/mlip.ipynb index 6c5ac75781..9e991ef3da 100644 --- a/doc/tutorials/mlip/mlip.ipynb +++ b/doc/tutorials/mlip/mlip.ipynb @@ -171,12 +171,14 @@ "\n", "In this tutorial, we will use [Simplified Molecular Input Line Entry System](https://doi.org/10.1021/ci00057a005) (SMILES) representations and [RDKit](https://github.com/rdkit/rdkit) to generate starting structures.\n", "SMILES strings are text representations of molecules without the hydrogen atoms.\n", - "For example, the SMILES strings for methane, ethane, and propane are \"C\", \"CC\", and \"CCC\" respectively.\n", + "For example, the SMILES strings for methane, ethane, and propane are `C`, `CC`, and `CCC` respectively.\n", + "Other elements are represented by their atomic symbol.\n", + "Ammonia for example would be `N`.\n", "SMILES can also represent cyclic structures, double bonds, and other chemically relevant properties, but we won't get into this today.\n", "RDKit, the tool mentioned earlier, provides powerful utilities for creating 3D structures from these string representations of molecules.\n", "\n", - "In the first part of the tutorial, we will perform a geometry optimization on an ethanol molecule. \n", - "The chemical formula of ethanol is CH3CH2OH." + "In the first part of the tutorial, we will perform a geometry optimization for water. \n", + "The chemical formula of ethanol is H2O." ] }, { @@ -184,7 +186,7 @@ "metadata": {}, "source": [ "**Exercise:**\n", - "- Write down the SMILES string for ethanol in the variable `ethanol_smiles`, we will use it to create a molecular structure.\n", + "- Write down the SMILES string for water in the variable `smiles`, we will use it to create a molecular structure.\n", " Remember, you do not include hydrogens in a SMILES representation." ] }, @@ -197,14 +199,14 @@ "outputs": [], "source": [ "# SOLUTION CELL\n", - "ethanol_smiles = \"CCO\"" + "smiles = \"O\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Using the SMILES string, we can create a 2D sketch of the molecule:" + "Using the SMILES string, we can create a 2D sketch of the molecule using rdkit:" ] }, { @@ -213,7 +215,7 @@ "metadata": {}, "outputs": [], "source": [ - "m = rdkit.Chem.MolFromSmiles(ethanol_smiles)\n", + "m = rdkit.Chem.MolFromSmiles(smiles)\n", "rdkit.Chem.Draw.MolToImage(rdkit.Chem.AddHs(m))" ] }, @@ -230,7 +232,7 @@ "metadata": {}, "outputs": [], "source": [ - "ethanol_atoms = rdkit2ase.smiles2atoms(ethanol_smiles)\n", + "ethanol_atoms = rdkit2ase.smiles2atoms(smiles)\n", "ethanol_atoms" ] }, @@ -297,7 +299,7 @@ "- Set the timestep to 0.5 femtoseconds. Don't forget to convert the units to simulation units first.\n", "\n", "**Hint:**\n", - "- Use a box length of `16` Angstrom, which is also the length of the unit system.\n", + "- Use an arbitrary box length, which is also the length of the unit system and disable periodic boundary conditions.\n", "- For the conversion of the timestep you should be familiar with the [ESPResSo Units](https://espressomd.github.io/doc/introduction.html#on-units) and you can use pint for the conversion. You can get the magnitude of the timestep in the correct unit system with pints' [`m_as()`](https://pint.readthedocs.io/en/stable/api/base.html#pint.Quantity.m_as) method." ] }, @@ -308,7 +310,8 @@ "outputs": [], "source": [ "# SOLUTION CELL\n", - "system = espressomd.System(box_l=[16] * 3)\n", + "system = espressomd.System(box_l=[1, 1, 1])\n", + "system.periodicity = [False, False, False]\n", "system.cell_system.skin = 0.4\n", "system.time_step = (0.5 * ureg.fs).m_as(((1 * ureg.u * ureg.angstrom**2) / ureg.electron_volt) ** 0.5)" ] @@ -412,7 +415,9 @@ "metadata": {}, "source": [ "**Exercise:**\n", - "- Take a look at the energies and the molecule in the visualizer during the minimisation." + "- Take a look at the energies and the molecule in the visualizer during the minimisation.\n", + "- Why would you stop at $\\text{f}_{max} = 0.1~\\text{eV}\\cdot\\AA^{-1}$?\n", + "- Compare the final energy to the value of water in the materials project database [mp-697111](https://next-gen.materialsproject.org/materials/mp-697111/tasks/mp-697111?chemsys=H-O) of $\\approx -14.89$ eV. The `mp-697111` structure is a periodic ice crystal belonging to the space group Cmc2." ] }, { @@ -422,10 +427,12 @@ "outputs": [], "source": [ "# SOLUTION CELL\n", + "# Due to the PBC and interactions the energy of the ice crystal is lower than the energy of the isolated water molecule.\n", + "# Going beyond the given fmax would be outside the achievable accuracy of the MLIP.\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(energies)\n", "plt.xlabel(\"Frame\")\n", - "plt.ylabel(\"Energy (eV)\")\n", + "plt.ylabel(\"Energy / eV\")\n", "plt.show()" ] }, @@ -444,6 +451,7 @@ "source": [ "system.part.clear()\n", "system.thermostat.turn_off()\n", + "system.periodicity = [False, False, False]\n", "\n", "del vis.zndraw[:]" ] @@ -631,8 +639,7 @@ "energies = []\n", "for atoms in vis.zndraw:\n", " OH_1.append(atoms.get_distance(300, 305, mic=True))\n", - " OH_2.append(atoms.get_distance(304, 306, mic=True))\n", - " energies.append(atoms.get_potential_energy())" + " OH_2.append(atoms.get_distance(304, 306, mic=True))" ] }, { @@ -657,20 +664,7 @@ "**Exercise:**\n", "- Indentify the highlighted atoms in the visualization window, which are the two hydrogen atoms of the sulfuric acid.\n", "- Observe the moment they undergo the protonation event visually and also in the hydrogen-oxygen distance plots.\n", - "- Can you also see changes in the potential energy? Why is that?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(10, 6))\n", - "plt.plot(energies)\n", - "plt.xlabel(\"Frame\")\n", - "plt.ylabel(\"Energy (eV)\")\n", - "plt.show()" + "- Discuss the results and the potential for using MLIPs to aid in the simulation of chemical reactions." ] }, { @@ -713,7 +707,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.10.14" } }, "nbformat": 4,