{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulate an Interchange with LAMMPS\n", "\n", "In this example, we'll quickly construct an `Interchange` and then run a simulation in LAMMPS. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need an `Interchange` to get started, so let's put that together quickly. For more explanation on this process, take a look at the [packed_box] and [protein_ligand] examples.\n", "\n", "[packed_box]: https://github.com/openforcefield/openff-interchange/tree/main/examples/packed_box\n", "[protein_ligand]: https://github.com/openforcefield/openff-interchange/tree/main/examples/protein_ligand" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mdtraj as md\n", "import nglview\n", "from lammps import lammps\n", "from openff.toolkit.topology import Molecule, Topology\n", "from openff.toolkit.typing.engines.smirnoff import ForceField\n", "from openff.toolkit.utils import get_data_file_path\n", "from openmm.app import PDBFile\n", "\n", "from openff.interchange import Interchange\n", "\n", "# Read a structure from the Toolkit's test suite into a Topology\n", "pdbfile = PDBFile(get_data_file_path(\"systems/packmol_boxes/propane_methane_butanol_0.2_0.3_0.5.pdb\"))\n", "molecules = [Molecule.from_smiles(smi) for smi in [\"CCC\", \"C\", \"CCCCO\"]]\n", "off_topology = Topology.from_openmm(pdbfile.topology, unique_molecules=molecules)\n", "\n", "# Construct the Interchange with the OpenFF \"Sage\" force field\n", "interchange = Interchange.from_smirnoff(\n", " force_field=ForceField(\"openff-2.0.0.offxml\"),\n", " topology=off_topology,\n", ")\n", "interchange.positions = pdbfile.positions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tada! A beautiful solvent system:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interchange.visualize(\"nglview\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Run a simulation\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we export a `.lmp` file that can be read by LAMMPS' `read_data` command:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interchange.to_lammps_datafile(\"interchange.lmp\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need to write an input file for LAMMPS. Parts of these input files depend on force field parameters, so we should use a sample input file written for our interchange as a starting point. We can generate such a sample file with the `to_lammps_input()` method:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interchange.to_lammps_input(\"interchange_pointenergy.in\", data_file=\"interchange.lmp\")\n", "with open(\"interchange_pointenergy.in\") as f:\n", " print(f.read())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the `read_data` line must match the data file name - in this case, `interchange.lmp`. That's why `to_lammps_input` needs the data file name. We could alternatively use the `to_lammps` method to write both files in a consistent way:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interchange.to_lammps(\"interchange\") # writes interchange.lmp and interchange_pointenergy.in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That sample file will only perform a single point energy calculation; here's a more complete file that includes the above parameters but will run an actual MD simulation. \n", "\n", "
\n", " ⚠️ Don't use example input files in production
\n", " Note that this is still just an example! You should carefully write and inspect input files for all production simulations.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lammps_in = \"\"\" # These commands may not be appropriate for all systems\n", "units real\n", "atom_style full\n", "\n", "# PBC in 3 dimensions\n", "dimension 3\n", "boundary p p p\n", "\n", "# Bonded interactions in Sage force field\n", "bond_style hybrid harmonic\n", "angle_style hybrid harmonic\n", "dihedral_style hybrid fourier\n", "improper_style cvff\n", "special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333333333\n", "\n", "# Non-bonded interactions in Sage force field\n", "pair_style lj/cut/coul/long 9.0 9.0\n", "pair_modify mix arithmetic tail yes\n", "\n", "# Load the parameterized system\n", "read_data interchange.lmp\n", "\n", "# Thermostat and velocity generation\n", "fix 3 all nvt temp 300.0 300.0 500\n", "velocity all create 300.0 29348 mom yes rot yes\n", "\n", "# Output control\n", "dump traj all dcd 10 traj.dcd\n", "thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe\n", "\n", "# PME electrostatics in Sage force field\n", "kspace_style pppm 1e-6\n", "\n", "# Run for 1000 steps at 2 fs δt\n", "timestep 2\n", "run 1000\n", "\n", "\"\"\"\n", "\n", "with open(\"lammps.in\", \"w\") as f:\n", " f.write(lammps_in)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll use LAMMPS' Python interface to run the simulation. `lmp.file(\"lammps.in\")` is equivalent to `lammps -in lammps.in` on the command line. This will run in serial, which is fine for so few steps and such a small system." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lmp = lammps()\n", "\n", "lmp.file(\"lammps.in\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we visualize! We can construct an MDTraj `Topology` from our Interchange by using OpenMM as a lingua franca. LAMMPS produces coordinates that are in the central unit cell, so for a simple system like this we just need to make molecules whole to visualize:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj = md.load(\n", " \"traj.dcd\",\n", " top=md.Topology.from_openmm(interchange.topology.to_openmm()),\n", ")\n", "traj.make_molecules_whole()\n", "nglview.show_mdtraj(traj)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "category": "parametrization_evaluation", "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.11.0" } }, "nbformat": 4, "nbformat_minor": 4 }