{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# QM/MM Simulations with ASE\n\nQM/MM Simulations couple two (or, in principle, more)\ndescriptions to get total energy\nand forces for the entire system in an efficient manner.\nASE has a native Explicit Interaction calculator,\n:class:`~ase.calculators.qmmm.EIQMMM`, that uses an electrostatic embedding\nmodel to couple the subsystems explicitly. See\n:doi:`the method paper for more info <10.1021/acs.jctc.7b00621>`.\n\nExamples of what this code has been used for can be seen\n:doi:`here <10.1021/jz500850s>`,\nand :doi:`here <10.1021/acs.inorgchem.6b01840>`.\n\nThis section will show you how to setup up various QM/MM simulations.\nWe will be using GPAW_ for the QM part. Other QM calculators should\nbe straightforwardly compatible with the subtractive-scheme SimpleQMMM\ncalculator, but for the Excplicit Interaction EIQMMM calculator, you\nwould need to be able to put an electrostatic external potential into\nthe calculator for the QM subsystem. This is often simply a matter of:\n\n1. Making the ASE-calculator write out the positions and charge-values\n to a format that your QM calculator can parse.\n2. Read in the forces on the point charges from the QM density.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ASE-calculators that currently support EIQMM:\n\n1. GPAW_\n2. DFTBplus_\n3. CRYSTAL_\n4. TURBOMOLE_\n\nTo see examples of how to make point charge potentials for EIQMMM,\nhave a look at the :class:`~ase.calculators.dftb.PointChargePotential`\nclasses in any of the calculators above.\n\n\nYou might also be interested in the solvent MM potentials included in ASE.\nThe tutorial on `tipnp water box equilibration` could be relevant to\nhave a look at. For acetonitrile, have a look at `acn_md_tutorial`.\n\nSome MD codes have more advanced solvators, such as AMBER_, and stand-alone\nprograms such as PACKMOL_ might also come in handy.\n\n\n\n## General Workflow\nThe implementation was developed with the focus of\nmodelling ions and complexes\nin solutions, we're working on expanding its functionality to encompass\nsurfaces.\n\nIn broad strokes, the steps to performing QM/MM MD simulations for thermal\nsampling or dynamics studies, these are the steps:\n\nQM/MM MD General Strategy for A QM complex in an MM solvent:\n\n1. Equilibrate an MM solvent box using one of the MM potentials built into\n ASE (see `tipnp water box equilibration` for water potentials), one\n of the compatible external MM codes, or write your own potential\n (see `Adding new calculators`)\n2. Optimize the gas-phase structure of your QM complex in GPAW, analyze what\n level of accuracy you will need for your task.\n3. Place the relaxed structure of the QM molecule in your MM solvent box,\n deleting overlapping MM molecules.\n4. Re-equillibrate the QM/MM system.\n5. Run production runs.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Electrostatic Embedding QM/MM\nTheoretical Background\n^^^^^^^^^^^^^^^^^^^^^^\n\nThe total energy expression for the full QM/MM system is:\n\n\\begin{align}E_\\mathrm{TOT} = E_\\mathrm{QM} + E_\\mathrm{I} + E_\\mathrm{MM}.\\end{align}\n\nThe MM region is modelled using point charge force fields, with charges\n$q_i$ and $\\tau_i$ denoting their spatial coordinates, so the\nQM/MM coupling term $E_\\mathrm{I}$ will be\n\n\\begin{align}E_\\mathrm{I} = \\sum_{i=1}^C q_i \\int \\frac{n({\\bf r})}{\\mid\\!{\\bf r}\n - \\tau_i\\!\\mid}\\mathrm{d}{\\bf r} +\n \\sum_{i=1}^C\\sum_{\\alpha=1}^A\n \\frac{q_i Z_{\\alpha}}{\\mid\\!{\\bf R}_\\alpha\n - \\tau_i\\!\\mid} + E_\\mathrm{RD}\\end{align}\n\nwhere $n({\\bf r})$ is the spatial electronic density of the quantum\nregion, $Z_\\alpha$ and ${\\bf R}_\\alpha$ are the charge and\ncoordinates of the nuclei in the QM region, respectively, and\n$E_\\mathrm{RD}$ is the term describing the remaining, non-Coulomb\ninteractions between the two subsystems.\n\nFor the MM point-charge external potential in GPAW, we use the total pseudo-\ncharge density $\\tilde{\\rho}({\\bf r})$ for the coupling, and since the\nColoumb integral is evaluated numerically on the real space grid, thus the\ncoupling term ends up like this:\n\n\\begin{align}E_\\mathrm{I} = \\sum_{i=1}^C q_i \\sum_{g}\n \\frac{\\tilde{\\rho}({\\bf r})}{\\mid\\!{\\bf r}_g - \\tau_i\\!\\mid} v_g\n + E_\\mathrm{RD}\\end{align}\n\nCurrently, the term for $E_{\\mathrm{RD}}$ implemented is a Lennard-\nJones-type potential:\n\n\\begin{align}E_\\mathrm{RD} = \\sum_i^C \\sum_\\alpha^A\n 4\\epsilon\\left[\n \\left(\\frac{\\sigma}{\\mid\\!{\\bf R}_\\alpha\n - \\tau_i\\!\\mid}\\right)^{12}\n - \\left(\\frac{\\sigma}{\\mid\\!{\\bf R}_\\alpha\n - \\tau_i\\!\\mid}\\right)^{6} \\right]\\end{align}\n\n### Example: QM/MM for a Water Dimer\nLet's first do a very simple electrostatic embedding QM/MM single point\nenergy calculation on the water dimer. The necessary inputs are described in\nthe :class:`ase.calculators.qmmm.EIQMMM` class.\n\n\nThe following script will calculate the QM/MM single point energy of the\nwater dimer from the `s22`, using LDA and TIP3P,\nfor illustration purposes.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\nfrom gpaw import GPAW\n\nfrom ase.calculators.qmmm import EIQMMM, Embedding, LJInteractions\nfrom ase.calculators.tip3p import TIP3P, epsilon0, sigma0\nfrom ase.data import s22\nfrom ase.visualize.plot import plot_atoms\n\n# Create system\natoms = s22.create_s22_system('Water_dimer')\natoms.center(vacuum=4.0)\n\n# visualization\nfig, ax = plt.subplots()\nplot_atoms(atoms, ax=ax)\nax.set_axis_off()\n\n# Make QM atoms selection of first water molecule:\nqm_idx = range(3)\n\n# Set up interaction & embedding object\ninteraction = LJInteractions({('O', 'O'): (epsilon0, sigma0)})\nembedding = Embedding(rc=0.02) # Short range analytical potential cutoff\n\n# Set up calculator\natoms.calc = EIQMMM(\n qm_idx,\n GPAW(txt='qm.out', mode='fd', xc='LDA'),\n TIP3P(),\n interaction,\n embedding=embedding,\n vacuum=None, # if None, QM cell = MM cell\n output='qmmm.log',\n)\n\nprint(f'The potential energy of the system is {atoms.get_potential_energy()}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we have just used the TIP3P LJ parameters for the QM part as well. If\nthis is a good idea or not isn't trivial. The LJInteractions module needs\ncombined parameters for all possible permutations of atom types in your\nsystem, that have LJ parameters. A list of combination rules can be found\n[here](http://www.sklogwiki.org/SklogWiki/index.php/Combining_rules).\nAn overview over different TIP3P LJ parameters can be found\ne.g., on the [LAMMPS website](https://docs.lammps.org/Howto_tip3p.html).\nHere's a code snippet of how to combine LJ parameters of atom types A and B\nvia the Lorentz-Berthelot rules:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import itertools as it\n\nepsHH = 0.0\nsigHH = 1.0\nepsOO = epsilon0\nsigOO = sigma0\n\nparameters = {'H': (epsHH, sigHH), 'O': (epsOO, sigOO)}\n\n\ndef lorenz_berthelot(p):\n combined = {}\n for comb in it.product(p.keys(), repeat=2):\n combined[comb] = (\n (p[comb[0]][0] * p[comb[1]][0]) ** 0.5,\n (p[comb[0]][1] + p[comb[1]][1]) / 2,\n )\n return combined\n\n\ncombined = lorenz_berthelot(parameters)\ninteraction = LJInteractions(combined)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will (somewhat redundantly) yield:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "print(combined)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to run structural relaxations and molecular dynamics\nusing the electrostatic embedding scheme:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ase.constraints import FixBondLengths\nfrom ase.optimize import LBFGS\n\nmm_bonds = [(3, 4), (4, 5), (5, 3)]\natoms.constraints = FixBondLengths(mm_bonds)\ndyn = LBFGS(atoms=atoms, trajectory='dimer.traj')\n# Relaxation (we only relax to very loose fmax here,\n# lower values should be used for serious calculations)\ndyn.run(fmax=0.5) # e.g. fmax=0.05" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since TIP3P is a rigid potential, we constrain all interatomic distances.\nQM bond lengths can be constrained too, in the same manner.\n\n### Periodicity\nFor these types of simulations with GPAW, you'd probably want two cells:\na QM (non-periodic) and and MM cell (periodic):\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "atoms.set_pbc(True)\n# Set up calculator\natoms.calc = EIQMMM(\n qm_idx,\n GPAW(txt='qm.out', mode='fd', xc='LDA'),\n TIP3P(),\n interaction,\n embedding=embedding,\n vacuum=4.0, # Now QM cell has walls min. 4 \u00c5 from QM atoms\n output='qmmm.log',\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will center the QM subsystem in the MM cell. For QM codes with no single\nreal-space grid like GPAW, you can still use this to center your QM subsystem,\nand simply disregard the QM cell, or manually center your QM subsystem,\nand leave vacuum as ``None``.\n\n## LJInteractionsGeneral - For More Intricate Systems\nIt often happens that you will have different 'atom types' (an element in a\nspecific environment) per element in your system, i.e. you want to assign\ndifferent LJ-parameters to the oxygens of your solute molecule and the oxygens\nof water. This can be done using\n:class:`~ase.calculators.qmmm.LJInteractionsGeneral`, which takes in NumPy\narrays with sigma and epsilon values for each individual QM and MM atom,\nrespectively, and combines them itself, with Lorentz-Berthelot.\nI.e., for our water dimer from before:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ase.calculators.qmmm import LJInteractionsGeneral\nfrom ase.calculators.tip3p import epsilon0, sigma0\n\n# General LJ interaction object for the 'OHHOHH' water dimer\nsigma_mm = np.array([sigma0, 0, 0]) # Hydrogens have 0 LJ parameters\nepsilon_mm = np.array([epsilon0, 0, 0])\nsigma_qm = np.array([sigma0, 0, 0])\nepsilon_qm = np.array([epsilon0, 0, 0])\ninteraction = LJInteractionsGeneral(\n sigma_qm,\n epsilon_qm,\n sigma_mm,\n epsilon_mm,\n qm_molecule_size=3,\n mm_molecule_size=3,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ``qm_molecule_size`` and ``mm_molecule_size`` should\nbe the number of atoms\nper molecule. Often the ``qm_molecule_size`` will\nsimply be the total number of\natoms in your QM subsystem. Here, our MM subsystem is comprised of a single\nwater molecule, but say we had N water molecules in the MM subsystem,\nwe wouldn't need to repeat e.g. the ``sigma_mm`` array N times,\nwe can simply keep it as it is\nwritten in the above.\n\n\n## EIQMMM and Charged Systems - Counterions\nIf your QM subsystem is charged, it is good to charge-neutrialize the entire\nsystem. This can be done in ASE by adding MM 'counterions', which are simple,\nsingle-atomic particles that carry a charge,\nand interact with the solvent and solute\nthrough a Coulomb and an LJ term.\nThe implementation is rather simplified and should only\nserve to neutralize the total system. It might be a good idea to restrain them\nso they do not diffuse too close to the QM subsystem,\nas the effective concentration\nin your simulation cell might be a lot higher than in real life.\n\n### Example: QM/MM for a NaCl in Water\nAs simple example we add here two NaCl to our system,\nwhere we treat the 2 Na as part of the QM region and the Cl as part\nof the MM region.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ase import Atoms\n\nions = Atoms('Na2Cl2', positions=[[1, 4, 3], [2, 8, 6], [10, 4, 2], [7, 8, 6]])\natoms += ions\n\n# resort so that the regions are together\natoms = Atoms([atoms[i] for i in [0, 1, 2, 6, 7, 3, 4, 5, 8, 9]])\n\nna_idx = [i.index for i in atoms if i.number == 11]\nqm_idx = list(qm_idx) + na_idx\n\n# visualization\nfig, ax = plt.subplots()\nplot_atoms(atoms, ax=ax)\nax.set_axis_off()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use the implementation, you need to 'Combine'\ntwo MM calculators, one for the\ncounterions, and one for your solvent. This is an example of combining two Cl-\nions with TIP3P water, using the TIP3P LJ-arrays from the previous section:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ase import units\nfrom ase.calculators.combine_mm import CombineMM\nfrom ase.calculators.counterions import AtomicCounterIon as ACI\n\n# Cl-: 10.1021/ct600252r\nsigCl = 4.02\nepsCl = 0.71 * units.kcal / units.mol\n\n# in this sub-atoms object, CombineMM only sees Cl and Water,\n# and Cl is here atom 3 and 4\nmmcalc = CombineMM(\n [3, 4], # indices of the counterion atoms in the MM region\n apm1=1,\n apm2=3, # atoms per 'molecule' of each subgroup\n calc1=ACI(-1, epsCl, sigCl), # Counterion calculator\n calc2=TIP3P(), # Water calculator\n sig1=[sigCl],\n eps1=[epsCl], # LJ Params for subgroup1\n sig2=sigma_mm,\n eps2=epsilon_mm,\n) # LJ params for subgroup2\n\natoms_mm = Atoms([i for i in atoms if i.index > 2 and i.number != 11])\nprint(atoms_mm)\natoms_mm.calc = mmcalc\nprint('MM', atoms_mm.get_potential_energy())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The charge of the counterions is defined as -1 in the first\ninput in ``ACI``, which\nThen also takes LJ-parameters for interactions with other ions\nbeloning to this calculator.\n\nThis ``mmcalc`` object is then used in the initialization\nof the ``EIQMMM`` calculator.\nBut before we can do that, the QM/MM Lennard-Jones potential\nneeds to understand that\nThe total MM subsystem is now comprised of two subgroups,\nthe counterions and the water.\nThat is done by initializing the ``interaction`` object with a\ntuple of NumPy arrays for the MM\nPart. So if your QM subsystem has 5 atoms, you'd do:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Na+: DOI 10.1021/ct600252r\nsigNa = 4.07\nepsNa = 0.0005 * units.kcal / units.mol\n\nsigma_qm = np.array([sigma0, 0, 0, sigNa, sigNa])\nepsilon_qm = np.array([epsilon0, 0, 0, epsNa, epsNa])\nsigma_mm = (np.array([sigma0, 0, 0]), np.array([sigCl]))\nepsilon_mm = (np.array([epsilon0, 0, 0]), np.array([epsCl]))\n\ninteraction = LJInteractionsGeneral(\n sigma_qm,\n epsilon_qm,\n sigma_mm,\n epsilon_mm,\n qm_molecule_size=5,\n mm_molecule_size=5,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Current limitations:\n\n* No QM PBCs\n* There is currently no automated way of doing QM/MM over covalent bonds\n (inserting cap-atoms, redistributing forces ...)\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Other tips:\n\nIf you are using GPAW and water, consider having a look at the\n[much faster RATTLE constraints for water here](https://gitlab.com/gpaw/gpaw/blob/master/gpaw/test/rattle.py)_\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.13.7" } }, "nbformat": 4, "nbformat_minor": 0 }