Skip to content

Commit

Permalink
update MLIP tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
@PythonFZ
PythonFZ committed Sep 30, 2024
1 parent 2cb91a6 commit bdb0d71
Showing 1 changed file with 23 additions and 29 deletions.
52 changes: 23 additions & 29 deletions doc/tutorials/mlip/mlip.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,20 +171,22 @@
"\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 CH<sub>3</sub>CH<sub>2</sub>OH."
"In the first part of the tutorial, we will perform a geometry optimization for water. \n",
"The chemical formula of ethanol is H<sub>2</sub>O."
]
},
{
"cell_type": "markdown",
"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."
]
},
Expand All @@ -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:"
]
},
{
Expand All @@ -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))"
]
},
Expand All @@ -230,7 +232,7 @@
"metadata": {},
"outputs": [],
"source": [
"ethanol_atoms = rdkit2ase.smiles2atoms(ethanol_smiles)\n",
"ethanol_atoms = rdkit2ase.smiles2atoms(smiles)\n",
"ethanol_atoms"
]
},
Expand Down Expand Up @@ -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."
]
},
Expand All @@ -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)"
]
Expand Down Expand Up @@ -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."
]
},
{
Expand All @@ -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()"
]
},
Expand All @@ -444,6 +451,7 @@
"source": [
"system.part.clear()\n",
"system.thermostat.turn_off()\n",
"system.periodicity = [False, False, False]\n",
"\n",
"del vis.zndraw[:]"
]
Expand Down Expand Up @@ -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))"
]
},
{
Expand All @@ -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."
]
},
{
Expand Down Expand Up @@ -713,7 +707,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
"version": "3.10.14"
}
},
"nbformat": 4,
Expand Down

0 comments on commit bdb0d71

Please sign in to comment.