{ "cells": [ { "cell_type": "markdown", "id": "b6817d5d-0866-4a1e-b620-58b29dd920f2", "metadata": {}, "source": [ "## Initial imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "ef6710ec-780d-45b6-ad77-b9f50025cadc", "metadata": {}, "outputs": [], "source": [ "from scm.plams import *" ] }, { "cell_type": "markdown", "id": "9d2387fb-904c-4ec4-96e9-60fb7a2bb603", "metadata": {}, "source": [ "## Elements, coordinates, lattice vectors, and charge" ] }, { "cell_type": "markdown", "id": "88e8c36e-f59b-4ab4-8884-ed17ca1253a2", "metadata": {}, "source": [ "### Manual molecule definition" ] }, { "cell_type": "code", "execution_count": 2, "id": "5e91e4e1-b795-405c-a832-18b29a0af3fd", "metadata": {}, "outputs": [], "source": [ "molecule = Molecule()\n", "molecule.add_atom(Atom(symbol='O', coords=(0,0,0)))\n", "molecule.add_atom(Atom(symbol='H', coords=(1,0,0)))\n", "molecule.add_atom(Atom(symbol='H', coords=(0,1,0)))" ] }, { "cell_type": "markdown", "id": "6c493a49-eb6a-42d2-b726-033cb344fcf0", "metadata": {}, "source": [ "To see the input that will be passed to AMS, create an AMSJob and print the input:" ] }, { "cell_type": "code", "execution_count": 3, "id": "6b29c9d5-8358-4dcf-908f-c23b62c13ead", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "ec4268e2-b503-440b-8d07-a10f8d5d9fb9", "metadata": {}, "source": [ "### Lattice vectors: 1D-periodic\n", "\n", "For periodic systems in 1 dimension, the lattice vector must be along the x direction (with 0 components along y and z)" ] }, { "cell_type": "code", "execution_count": 4, "id": "3732eb50-66f0-479e-8be3-f3094a0c777e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", " Lattice\n", " 10.0000000000 0.0000000000 0.0000000000\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.lattice = [\n", " [10, 0, 0]\n", "]\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "f7b5fa8d-9be5-466b-8a15-0ec05d732ee4", "metadata": {}, "source": [ "### Lattice vectors: 2D-periodic\n", "\n", "For 2 dimensions, the two lattice vectors must lie in the xy plane (with 0 component along z)." ] }, { "cell_type": "code", "execution_count": 5, "id": "cef68617-74b1-4a53-abf9-718e3fb9ccd5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", " Lattice\n", " 10.0000000000 0.0000000000 0.0000000000\n", " 0.0000000000 11.0000000000 0.0000000000\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.lattice = [\n", " [10, 0, 0],\n", " [0, 11, 0],\n", "]\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "d8f5e01c-cf15-4808-b11c-ab70e514aef3", "metadata": {}, "source": [ "### Lattice vectors: 3D-periodic" ] }, { "cell_type": "code", "execution_count": 6, "id": "fbf61735-c2b2-497c-a64b-7c97e91ac6eb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", " Lattice\n", " 10.0000000000 0.0000000000 0.0000000000\n", " 0.0000000000 11.0000000000 0.0000000000\n", " -1.0000000000 0.0000000000 12.0000000000\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.lattice = [\n", " [10, 0, 0],\n", " [0, 11, 0],\n", " [-1, 0, 12]\n", "]\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "4d33f54b-6088-44e4-a3c6-9c918350e85f", "metadata": {}, "source": [ "### Delete lattice vectors" ] }, { "cell_type": "code", "execution_count": 7, "id": "b63866ef-6d97-4c52-9d2c-a48de74632fa", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.lattice = []\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "15234e4f-78f5-4d33-8c7e-44842667e427", "metadata": {}, "source": [ "### Charge" ] }, { "cell_type": "code", "execution_count": 8, "id": "6d3764e5-7f9d-47b3-b3f2-6843935eb5f3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", " Charge -1\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.properties.charge = -1\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "c02683b1-9193-41b0-a874-9d39c4d31b8c", "metadata": {}, "source": [ "To get the charge of a molecule, use ``molecule.properties.get(\"charge\", 0)``. If the charge is not defined you will then get 0 as the charge." ] }, { "cell_type": "code", "execution_count": 9, "id": "bea500ce-dab7-493c-b708-f1cd3c0de244", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The charge is -1\n" ] } ], "source": [ "my_charge = molecule.properties.get(\"charge\", 0)\n", "print(f\"The charge is {my_charge}\")" ] }, { "cell_type": "markdown", "id": "2ff4b6ec-e7fc-44bd-8563-d77ca994ca6a", "metadata": {}, "source": [ "Unset the charge:" ] }, { "cell_type": "code", "execution_count": 10, "id": "964da05c-be27-4ca0-997a-a099807d0ec3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The charge is 0\n" ] } ], "source": [ "if 'charge' in molecule.properties:\n", " del molecule.properties.charge\n", " \n", "my_charge = molecule.properties.get(\"charge\", 0)\n", "print(f\"The charge is {my_charge}\")" ] }, { "cell_type": "markdown", "id": "cdc5ce93-8dd8-4502-b827-6485e050a8d6", "metadata": {}, "source": [ "## Atomic properties: masses, regions, force field types ...\n", "\n", "In the AMS system block most atomic properties are given as a suffix at the end of the line.\n", "\n", "To access an individual atom, use for example ``molecule[1]``, which corresponds to the first atom. **Note that the indexing starts with 1**, unlike normal Python lists that start with 0!" ] }, { "cell_type": "markdown", "id": "490de15c-9919-4747-8215-5baaa3459a75", "metadata": {}, "source": [ "### Isotopes (atomic masses)" ] }, { "cell_type": "code", "execution_count": 11, "id": "1329c211-26e2-4b68-b6df-ea78c106b131", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 mass=2.014\n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule[2].properties.mass = 2.014\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "d72cf84b-c8fe-436c-aa75-f3b66a645ff0", "metadata": {}, "source": [ "### Regions\n", "\n", "Regions are used for example to \n", "\n", "* set special basis sets on a subset of atoms, or \n", "* apply a thermostat in molecular dynamics to only a subset of atoms, \n", "* visualize atoms easily in the AMS GUI,\n", "* and much more!\n", "\n", "Use Python sets to specify regions. In this way, one atom can belong to multiple regions." ] }, { "cell_type": "code", "execution_count": 12, "id": "dd08e942-d9f9-4e58-81ed-bd4e7828ff91", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 region=region1\n", " H 1.0000000000 0.0000000000 0.0000000000 mass=2.014 region=region1\n", " H 0.0000000000 1.0000000000 0.0000000000 region=region1,region2\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule[1].properties.region = {\"region1\"}\n", "molecule[2].properties.region = {\"region1\"}\n", "molecule[3].properties.region = {\"region1\", \"region2\"}\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "05cdb58c-51f4-4367-93d0-8d8d31496b53", "metadata": {}, "source": [ "### Force field types\n", "\n", "Some force fields need to know the specific atom type and not just the chemical element. Use ``ForceField.Type`` for this when you use the ForceField engine:" ] }, { "cell_type": "code", "execution_count": 13, "id": "75c951cf-d221-4fb9-adb4-2c894b56c427", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 ForceField.Type=OW region=region1\n", " H 1.0000000000 0.0000000000 0.0000000000 ForceField.Type=HW mass=2.014 region=region1\n", " H 0.0000000000 1.0000000000 0.0000000000 ForceField.Type=HW region=region1,region2\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule[1].properties.ForceField.Type = \"OW\" # these types would depend on what type of force field you use!\n", "molecule[2].properties.ForceField.Type = \"HW\"\n", "molecule[3].properties.ForceField.Type = \"HW\"\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "00460a3e-be07-43a1-94c2-0003271970ba", "metadata": {}, "source": [ "### Delete all atom-specific options\n", "Loop over the atoms and set ``atom.properties`` to an empty ``Settings()``:" ] }, { "cell_type": "code", "execution_count": 14, "id": "2b6c82d6-b5fd-4f06-8cfa-498b9b763821", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "for at in molecule:\n", " at.properties = Settings()\n", " \n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "80ffc47a-e4f1-4b7b-8d9e-8b74f6656ca3", "metadata": {}, "source": [ "## Bonds\n", "\n", "Most methods (DFT, DFTB, ML Potential, ReaxFF) ignore any specified bonds. \n", "\n", "When using force fields, you sometimes need to specify the bonds that connect atoms. Some force fields (UFF, GAFF) can automatically guess the correct types of bonds. \n", "\n", "So **most of the time you do not manually need to specify bonds**.\n", "\n", "If you **need** to specify bonds, it is easiest \n", "\n", "* to handle in the AMS GUI: use File -> Export Coordinates -> .in, and then load the file with ``molecule = Molecule(\"my_file.in\")``\n", "* to use the ``from_smiles`` function to generate a molecule from SMILES code, for example ``molecule = from_smiles(\"O\")`` for water.\n", "\n", "If you need to add bonds manually in PLAMS you can do it as follows:" ] }, { "cell_type": "code", "execution_count": 15, "id": "a2e36cfa-8419-418e-950b-71ee1cc2325a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", " BondOrders\n", " 1 2 1.0\n", " 1 3 1.0\n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "molecule.add_bond(molecule[1], molecule[2], order=1.0)\n", "molecule.add_bond(molecule[1], molecule[3], order=1.0)\n", "print(AMSJob(molecule=molecule).get_input())" ] }, { "cell_type": "markdown", "id": "b524e108-ee87-47d4-b37d-aafc04314169", "metadata": {}, "source": [ "## Multiple systems\n", "\n", "Some tasks like NEB (nudged elastic band) require more than 1 system in the input file. This can be accomplished by using a Python dictionary.\n", "\n", "In AMS, \n", "\n", "* the \"main system\" has no name. It should have the key ``\"\"`` (empty string) in the dictionary.\n", "\n", "* every additional system needs to have a name, that is used as the key in the dictionary.\n", "\n", "Let's first define two ``Molecule`` in the normal way:" ] }, { "cell_type": "code", "execution_count": 16, "id": "9f056427-36a9-4d6b-9951-63d1873cc623", "metadata": {}, "outputs": [], "source": [ "molecule1 = Molecule()\n", "molecule1.add_atom(Atom(symbol='O', coords=(0,0,0)))\n", "molecule1.add_atom(Atom(symbol='H', coords=(1,0,0)))\n", "molecule1.add_atom(Atom(symbol='H', coords=(0,1,0)))\n", "\n", "molecule2 = Molecule()\n", "molecule2.add_atom(Atom(symbol='O', coords=(0,0,0)))\n", "molecule2.add_atom(Atom(symbol='H', coords=(3.33333,0,0)))\n", "molecule2.add_atom(Atom(symbol='H', coords=(0,5.555555,0)))" ] }, { "cell_type": "markdown", "id": "9bd03e4d-2623-4043-b16a-68d7a2aeff38", "metadata": {}, "source": [ "Then create the ``mol_dict`` dictionary:" ] }, { "cell_type": "code", "execution_count": 17, "id": "8296d8ae-5099-4aaa-b545-726849f76abf", "metadata": {}, "outputs": [], "source": [ "mol_dict = {\n", " \"\": molecule1, # main system, empty key (no name)\n", " \"final\": molecule2 # for NEB, use \"final\" as the name for the other endpoint\n", "}" ] }, { "cell_type": "markdown", "id": "4e0e2482-9b42-4570-9740-db1ed4bdd346", "metadata": {}, "source": [ "Pass the ``mol_dict`` as the ``molecule`` argument to ``AMSJob``:" ] }, { "cell_type": "code", "execution_count": 18, "id": "5309858b-c9f1-4c36-a6bd-9c7bf73f2d94", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "system\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 1.0000000000 0.0000000000 0.0000000000 \n", " H 0.0000000000 1.0000000000 0.0000000000 \n", " End\n", "End\n", "system final\n", " Atoms\n", " O 0.0000000000 0.0000000000 0.0000000000 \n", " H 3.3333300000 0.0000000000 0.0000000000 \n", " H 0.0000000000 5.5555550000 0.0000000000 \n", " End\n", "End\n", "\n", "\n" ] } ], "source": [ "print(AMSJob(molecule=mol_dict).get_input())" ] }, { "cell_type": "markdown", "id": "b9c36083-be5d-4173-a094-636d43af16fc", "metadata": {}, "source": [ "Above we see that the main system is printed just as before. A second system block \"system final\" is also added with ``molecule2``." ] } ], "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.8.12" } }, "nbformat": 4, "nbformat_minor": 5 }