01-Structure-Representation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Structure Representations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [DiffPy-CMI](http://www.diffpy.org/products/diffpycmi/) framework provides several options for representing atomic structures of materials.  Here we demonstrate basic operations on structure models using the [diffpy.structure](https://github.com/diffpy/diffpy.structure) and [pyobjcryst](https://github.com/diffpy/pyobjcryst) Python packages included in the DiffPy-CMI software.\n",
    "\n",
    "**Contents**\n",
    "\n",
    "> [diffpy.structure](#diffpy.structure)<br>\n",
    "> [pyobjcryst](#pyobjcryst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## diffpy.structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In [diffpy.structure](https://github.com/diffpy.structure) the atom arrangements are represented as a collaboration of *Structure*, *Atom* and *Lattice* classes.  The *Structure* class is an enhanced Python list of *Atom* objects, where each atom stores fractional coordinates, element or ion symbol, a matrix of displacement parameters and other attributes.  The *Lattice* class defines fractional coordinates with respect to the absolute Cartesian coordinate system and provides functions for conversion between fractional and Cartesian coordinates and other geometric operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lattice=Lattice(a=4.123, b=4.123, c=4.123, alpha=90, beta=90, gamma=90)\n",
      "Cs   0.000000 0.000000 0.000000 1.0000\n",
      "Cl   0.500000 0.500000 0.500000 1.0000\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from diffpy.structure import Structure, Atom, Lattice\n",
    "\n",
    "acs = Atom('Cs', [0, 0, 0])\n",
    "acl = Atom('Cl', [0.5, 0.5, 0.5])\n",
    "cscl = Structure(atoms=[acs, acl],\n",
    "                 lattice=Lattice(4.123, 4.123, 4.123, 90, 90, 90))\n",
    "\n",
    "print(cscl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The atoms in the `cscl` object can be accessed by integer indices starting from zero.  A subset of atoms can be selected using an index range, but also using a list of indices or a mask array of boolean flags."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cs   0.000000 0.000000 0.000000 1.0000 \n",
      "---\n",
      "lattice=Lattice(a=4.123, b=4.123, c=4.123, alpha=90, beta=90, gamma=90)\n",
      "Cl   0.500000 0.500000 0.500000 1.0000\n",
      "Cs   0.000000 0.000000 0.000000 1.0000 \n",
      "---\n",
      "lattice=Lattice(a=4.123, b=4.123, c=4.123, alpha=90, beta=90, gamma=90)\n",
      "Cl   0.500000 0.500000 0.500000 1.0000 \n",
      "---\n",
      "lattice=Lattice(a=4.123, b=4.123, c=4.123, alpha=90, beta=90, gamma=90)\n",
      "Cs   0.000000 0.000000 0.000000 1.0000\n"
     ]
    }
   ],
   "source": [
    "print(cscl[0], '\\n---')     # first atom in the structure\n",
    "print(cscl[::-1], '\\n---')  # reversed order of atoms\n",
    "print(cscl[[1,]], '\\n---')  # substructure containing only the 2nd atom\n",
    "print(cscl[cscl.x < 0.2])   # substructure of atoms with x < 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Atom positions in the structure are specified in fractional coordinates; therefore, their Cartesian positions and relative distances change with a change in lattice parameters.  Here we show this behavior using the *distance* function and the integer indices of the relevant atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cl Cartesian coordinates: [2.0615 2.0615 2.0615]\n",
      "Cs-Cl distance: 3.570622739803241\n",
      "# change cell parameter to a = 3\n",
      "Cl Cartesian coordinates: [1.5 1.5 1.5]\n",
      "Cs-Cl distance: 2.598076211353316\n"
     ]
    }
   ],
   "source": [
    "print(\"Cl Cartesian coordinates:\", cscl[1].xyz_cartn)\n",
    "print(\"Cs-Cl distance:\", cscl.distance(0, 1))\n",
    "print(\"# change cell parameter to a = 3\")\n",
    "cscl.lattice.setLatPar(a=3, b=3, c=3)\n",
    "print(\"Cl Cartesian coordinates:\", cscl[1].xyz_cartn)\n",
    "print(\"Cs-Cl distance:\", cscl.distance(0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Details about each site are stored as data attributes of the *Atom* object.  These attributes can be changed either individually per each atom or using mapped arrays of the owning *Structure*.  Thus, rather than using a for loop, a single statement can be used to set isotropic displacement parameters for all atoms or for a subset of chlorine atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.003 0.003]\n",
      "[0.003 0.004]\n"
     ]
    }
   ],
   "source": [
    "cscl.Uisoequiv = 0.003\n",
    "print(cscl.Uisoequiv)\n",
    "cscl[cscl.element == 'Cl'].Uisoequiv = 0.004\n",
    "print(cscl.Uisoequiv) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The changes in the equivalent isotropic displacement parameter *Uisoequiv* are propagated to the tensor of displacement parameters *U*.  The *anisotropy* flag specifies whether anisotropic displacements are allowed on each atom site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.003 0.    0.   ]\n",
      "  [0.    0.003 0.   ]\n",
      "  [0.    0.    0.003]]\n",
      "\n",
      " [[0.004 0.    0.   ]\n",
      "  [0.    0.004 0.   ]\n",
      "  [0.    0.    0.004]]] \n",
      "---\n",
      "[False False]\n"
     ]
    }
   ],
   "source": [
    "print(cscl.U, '\\n---')\n",
    "print(cscl.anisotropy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0.004 0.    0.   ]\n",
      "  [0.    0.003 0.   ]\n",
      "  [0.    0.    0.003]]\n",
      "\n",
      " [[0.004 0.    0.   ]\n",
      "  [0.    0.004 0.   ]\n",
      "  [0.    0.    0.004]]] \n",
      "---\n",
      "[ True False]\n"
     ]
    }
   ],
   "source": [
    "# allow anisotropy for the first atom\n",
    "cscl[0].anisotropy = True\n",
    "cscl[0].U11 = 0.004\n",
    "print(cscl.U, '\\n---')\n",
    "print(cscl.anisotropy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Structure data can be loaded from several file formats such as XYZ, PDB or CIF.  The *Structure* class has been designed as a simple list of atoms with no awareness of crystal symmetry.  Therefore, when loading from a CIF file, the asymmetric unit gets expanded to a full unit cell as if in the P1 symmetry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lattice=Lattice(a=5.62, b=5.62, c=5.62, alpha=90, beta=90, gamma=90)\n",
      "Na1+ 0.000000 0.000000 0.000000 1.0000\n",
      "Na1+ 0.000000 0.500000 0.500000 1.0000\n",
      "Na1+ 0.500000 0.000000 0.500000 1.0000\n",
      "Na1+ 0.500000 0.500000 0.000000 1.0000\n",
      "Cl1- 0.500000 0.500000 0.500000 1.0000\n",
      "Cl1- 0.500000 0.000000 0.000000 1.0000\n",
      "Cl1- 0.000000 0.500000 0.000000 1.0000\n",
      "Cl1- 0.000000 0.000000 0.500000 1.0000\n"
     ]
    }
   ],
   "source": [
    "from diffpy.structure import loadStructure\n",
    "nacl = loadStructure('NaCl.cif')\n",
    "print(nacl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Na1' 'Na1_2' 'Na1_3' 'Na1_4' 'Cl1' 'Cl1_2' 'Cl1_3' 'Cl1_4']\n"
     ]
    }
   ],
   "source": [
    "print(nacl.label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although symmetry operations are not intrinsic to the *Structure* class, the diffpy.structure package provides functions for generating symmetry equivalent positions or symmetry constraints for the coordinates and displacement parameters.  The package also provides definitions for all space groups in over 500 different symmetry settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[array([0., 0., 0.]),\n",
       "  array([0. , 0.5, 0.5]),\n",
       "  array([0.5, 0. , 0.5]),\n",
       "  array([0.5, 0.5, 0. ])]]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from diffpy.structure.spacegroups import GetSpaceGroup\n",
    "from diffpy.structure.symmetryutilities import ExpandAsymmetricUnit\n",
    "\n",
    "Fm3m = GetSpaceGroup('Fm-3m')\n",
    "eau = ExpandAsymmetricUnit(Fm3m, corepos=[[0, 0, 0],])\n",
    "eau.expandedpos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *SymmetryConstraints* class can be used to determine symmetry constraints on positions and displacement tensors at the specified unit cell sites."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "asymmetric unit and the equivalent positions: {0: [0, 1, 2, 3], 4: [4, 5, 6, 7]}\n",
      "position parameters: []\n",
      "position constraints on the 1-st site: {'x': '+0', 'y': '+0', 'z': '+0'}\n",
      "ADP parameters: [('U110', 0.0), ('U114', 0.0)]\n",
      "ADP constraints on the 1-st site: {'U11': 'U110', 'U22': 'U110', 'U33': 'U110', 'U12': '0', 'U13': '0', 'U23': '0'}\n"
     ]
    }
   ],
   "source": [
    "from diffpy.structure.symmetryutilities import SymmetryConstraints\n",
    "\n",
    "symcon = SymmetryConstraints(Fm3m, positions=nacl.xyz)\n",
    "\n",
    "print('asymmetric unit and the equivalent positions:', symcon.coremap)\n",
    "print('position parameters:', symcon.pospars)\n",
    "print('position constraints on the 1-st site:', symcon.poseqns[0])\n",
    "print('ADP parameters:', symcon.Upars)\n",
    "print('ADP constraints on the 1-st site:', symcon.Ueqns[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information about the diffpy.structure package, see the online documentation at http://www.diffpy.org/diffpy.structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## pyobjcryst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option for describing material structures with DiffPy-CMI is to use the pyobjcryst package.  [pyobjcryst](https://github.com/diffpy/pyobjcryst) is a Python interface to the [ObjCryst++](https://sourceforge.net/projects/objcryst/) crystallographic library developed by Vincent Favre-Nicolin.  pyobjcryst provides advanced features for representing crystal structures with intrinsic crystal symmetries, chemical constraints such as rigid atom groups, and restraints on bond lengths and bond and/or torsion angles.  The pyobjcryst *Crystal* class represents periodic structures as a collection of general scatterers located at fractional coordinates.  Each scatterer is a single atom site in simple case; however, it can also be a collection of atom positions grouped in a *Molecule* object or described by Z-Matrix.  The pyobjcryst package supports the standard CIF and a native XML format for loading and storing structure models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "UnitCell : Sodium chloride(F m -3 m)\n",
      "    Cell dimensions :  5.62000   5.62000   5.62000  90.00000  90.00000  90.00000\n",
      "List of scattering components (atoms): 2\n",
      "Na1              at :  0.0000 0.0000 0.0000, Occup=1.0000 * 0.0208, ScattPow:Na1             , Biso=    0.0000\n",
      "Cl1              at :  0.5000 0.5000 0.5000, Occup=1.0000 * 0.0208, ScattPow:Cl1             , Biso=    0.0000\n",
      "\n",
      "Occupancy = occ * dyn, where:\n",
      "        - occ is the 'real' occupancy\n",
      "        - dyn is the dynamical occupancy correction, indicating  either\n",
      "          an atom on a special position, or several identical atoms \n",
      "          overlapping (dyn=0.5 -> atom on a symetry plane / 2fold axis..\n",
      "                               -> OR 2 atoms strictly overlapping)\n",
      "\n",
      " Total number of components (atoms) in one unit cell : 8\n"
     ]
    }
   ],
   "source": [
    "from pyobjcryst import loadCrystal\n",
    "\n",
    "nacl = loadCrystal('NaCl.cif')\n",
    "print(nacl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *CreateCrystalFromCIF* function returned an instance of the pyobjcryst *Crystal* class, which has a variety of functions for accessing and manipulating structure data.  For example, the *GetSpaceGroup()* function provides access to the active space group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SpaceGroup:\n",
      "  Schoenflies symbol = Oh^5\n",
      "  Hermann-Maugin symbol = F m -3 m\n",
      "  Hall symbol = -F 4 2 3\n",
      "  SgNumber = 225\n",
      "  Number of Seitz Matrix = 24\n",
      "  Number of Translation Vectors = 4\n",
      "  List of Seitz Matrices : \n",
      "    x,y,z\n",
      "    -y,x,z\n",
      "    -x,-y,z\n",
      "    y,-x,z\n",
      "    x,-y,-z\n",
      "    y,x,-z\n",
      "    -x,y,-z\n",
      "    -y,-x,-z\n",
      "    z,x,y\n",
      "    -x,z,y\n",
      "    -z,-x,y\n",
      "    x,-z,y\n",
      "    z,-x,-y\n",
      "    x,z,-y\n",
      "    -z,x,-y\n",
      "    -x,-z,-y\n",
      "    y,z,x\n",
      "    y,-z,-x\n",
      "    z,y,-x\n",
      "    -y,z,-x\n",
      "    -z,-y,-x\n",
      "    -y,-z,x\n",
      "    z,-y,x\n",
      "    -z,y,x\n",
      "  There is an inversion center at 0 0 0\n",
      "  List of Translation vectors :\n",
      "     0,0,0\n",
      "     0,0.5,0.5\n",
      "     0.5,0,0.5\n",
      "     0.5,0.5,0\n",
      "Extension (origin choice, rhomboedral/hexagonal):\u0000\n"
     ]
    }
   ],
   "source": [
    "print(nacl.GetSpaceGroup())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The positions of the Na and Cl scatterers are specified in fractional coordinates; therefore, the bond distances and angles change with lattice parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ORIGINAL:\n",
      "0 0 0 1 0.0208333 0x7f8cd3936660 Na1\n",
      "0.5 0.5 0.5 1 0.0208333 0x7f8cd3936c60 Cl1\n",
      "Na-Cl distance: 2.81\n",
      "EXPANDED:\n",
      "0 0 0 1 0.0208333 0x7f8cd3936660 Na1\n",
      "0.5 0.5 0.5 1 0.0208333 0x7f8cd3936c60 Cl1\n",
      "Na-Cl distance: 3.5\n"
     ]
    }
   ],
   "source": [
    "from pyobjcryst.crystal import Crystal\n",
    "\n",
    "print(\"ORIGINAL:\")\n",
    "for sc in nacl.GetScatteringComponentList():\n",
    "    print(sc)\n",
    "print(\"Na-Cl distance:\", nacl.GetMinDistanceTable()[0, 1])\n",
    "a2 = 7\n",
    "nacl.a, nacl.b, nacl.c = a2, a2, a2\n",
    "print(\"EXPANDED:\")\n",
    "for sc in nacl.GetScatteringComponentList():\n",
    "    print(sc)\n",
    "print(\"Na-Cl distance:\", nacl.GetMinDistanceTable()[0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it is preferable to have bond distances and angles independent of lattice parameters - for example when searching for crystal packing of molecules of a known shape.  The pyobjcryst library allows this by using an entire molecule as a scattering unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ORIGINAL:\n",
      "-0.25 0 0 1 1 0x7f8cd39d1670 Carbon\n",
      "0.25 0 0 1 1 0x7f8cd39d1670 Carbon\n",
      "C1-C2 distance: 2.0\n",
      "EXPANDED:\n",
      "-0.125 0 0 1 1 0x7f8cd39d1670 Carbon\n",
      "0.125 0 0 1 1 0x7f8cd39d1670 Carbon\n",
      "C1-C2 distance: 2.0\n"
     ]
    }
   ],
   "source": [
    "from pyobjcryst.crystal import Crystal\n",
    "from pyobjcryst.molecule import Molecule\n",
    "from pyobjcryst.scatteringpower import ScatteringPowerAtom\n",
    "\n",
    "crst = Crystal(4, 4, 4, 'P1')\n",
    "# atom scattering properties are defined by the\n",
    "# ScatteringPower object.  The same ScatteringPower may\n",
    "# be shared among several atom sites.\n",
    "spC = ScatteringPowerAtom('Carbon', 'C')\n",
    "crst.AddScatteringPower(spC)\n",
    "# molecules in ObjCryst++ are owned by a Crystal object and\n",
    "# thus need to be created with a reference to their owner.\n",
    "mol = Molecule(crst, \"mol\")\n",
    "# atom positions are specified in Cartesian coordinates\n",
    "# anchored at the center of mass of the Molecule.\n",
    "mol.AddAtom(-1, 0, 0, spC, 'C1')\n",
    "mol.AddAtom(+1, 0, 0, spC, 'C2')\n",
    "# activate the Molecule as a scattering entity within the Crystal.\n",
    "crst.AddScatterer(mol)\n",
    "\n",
    "print(\"ORIGINAL:\")\n",
    "for sc in crst.GetScatteringComponentList():\n",
    "    print(sc)\n",
    "print(\"C1-C2 distance:\", crst.GetMinDistanceTable()[0, 1])\n",
    "\n",
    "# enlarge lattice parameters\n",
    "crst.a, crst.b, crst.c = 8, 8, 8\n",
    "print(\"EXPANDED:\")\n",
    "for sc in crst.GetScatteringComponentList():\n",
    "    print(sc)\n",
    "print(\"C1-C2 distance:\", crst.GetMinDistanceTable()[0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fractional coordinates of the carbon atoms were adjusted to preserve the same interatomic distances within the Molecule.  The position of a *Molecule* object within a *Crystal* is specified in fractional coordinate; changing this coordinate will affect all atoms in the molecule.  Molecules can also be rotated by changing their orientation Quaternion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MOVED:\n",
      "-0.125 0.5 0 1 1 0x7f8cd39d1670 Carbon\n",
      "0.125 0.5 0 1 1 0x7f8cd39d1670 Carbon\n",
      "ROTATED:\n",
      "-6.25e-09 0.375 0 1 1 0x7f8cd39d1670 Carbon\n",
      "6.25e-09 0.625 0 1 1 0x7f8cd39d1670 Carbon\n"
     ]
    }
   ],
   "source": [
    "from math import pi\n",
    "from pyobjcryst.molecule import Quaternion\n",
    "\n",
    "# quaternions for a null rotation and for a 90deg rotation by the z-axis\n",
    "qnone = Quaternion.RotationQuaternion(0, 0, 0, 1)\n",
    "qz90 = Quaternion.RotationQuaternion(0.5 * pi, 0, 0, 1)\n",
    "\n",
    "# (1) move the Molecule mol along the Y-axis\n",
    "mol.Y = 0.5\n",
    "mol.Q0, mol.Q1, mol.Q2, mol.Q3 = qnone.Q0, qnone.Q1, qnone.Q2, qnone.Q3\n",
    "print(\"MOVED:\")\n",
    "for sc in crst.GetScatteringComponentList():\n",
    "    print(sc)\n",
    "    \n",
    "# (2) rotate the molecule by 90 degrees around the Z-axis\n",
    "mol.Q0, mol.Q1, mol.Q2, mol.Q3 = qz90.Q0, qz90.Q1, qz90.Q2, qz90.Q3\n",
    "print(\"ROTATED:\")\n",
    "for sc in crst.GetScatteringComponentList():\n",
    "    print(sc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "For more information about the pyobjcryst package see http://www.diffpy.org/pyobjcryst/.<br>\n",
    "The underlying ObjCryst++ library is documented at http://vincefn.net/ObjCryst/."
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}