ase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Executing DFT simulations using ASE, and post-processing the output via McStas\n",
    "\n",
    "### Authors:\n",
    " - Mads Bertelsen (ESS)\n",
    " - Mousumi Upadhyay Kahaly (ELI-ALPS)\n",
    " - Shervin Nourbakhsh (ILL)\n",
    " - Gergely Nagy (ELI-ALPS)\n",
    "\n",
    "This notebook demonstrates a very simple, yet powerful, integrated workflow of executing a DFT-based crystal structure relaxation, and then using the result in a simulated neutron-scattering experiment. The DFT calculation is conducted via the open-source *[Quantum Espresso](https://www.quantum-espresso.org)* software, and the neurton experiment is simulated via *[McStas](https://www.mcstas.org)*. \n",
    "\n",
    "Two Python module is used to seamlessly interface between these softwares and the notebook. The *[Atomistic Simulation Environment](https://wiki.fysik.dtu.dk/ase/index.html)* (***ASE***), a very powerful tool for atomistic simulations in general, is used to obtain the initial structure from a public database, [*COD*](https://www.crystallography.net/cod/) and communicate with Quantum Espresso; and the [***McStasScript*** module](https://github.com/PaNOSC-ViNYL/McStasScript) is used to set up and execute the McStas simulation. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your Quantum-Espresso binaries are in a local folder, add it to `$PATH` here. If they are already there, you can skip the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QE_bin_path = os.environ[\"HOME\"]+\"/PANOSC/bin\"\n",
    "os.environ['PATH']=os.environ['PATH']+\":\"+QE_bin_path\n",
    "\n",
    "mcstas_outdir = \"mcstas_output\"\n",
    "os.environ['PATH']=os.environ['PATH']+\":/usr/lib64/mpich/bin:\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(os.environ['PATH'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set the path for the temporary files\n",
    "\n",
    "This folder will be used for temporary files created while running the simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmpdir='/tmp/jupiter/'\n",
    "print('Create temporary directory: '+tmpdir)\n",
    "os.makedirs(tmpdir,exist_ok=True)\n",
    "os.chdir(tmpdir)\n",
    "os.makedirs(mcstas_outdir,exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the input files for ASE and Quantum Espresso\n",
    "ASE can take in input several file formats.\n",
    "In this demo we will check that the conversion is done properly and that the output of the Quantum-Espresso (QE) simulation can be carried on with them.\n",
    "\n",
    "In the following only CIF files will be considered as inputs for the simulation workflow\n",
    "\n",
    "\n",
    "### Convert CIF to QE input file\n",
    "ASE is able to convert from different formats. If you plan to run QE as a standalone package you need to use files in its input format, so you need to convert for example a CIF file into QE format. This can be done in the following way\n",
    "```ase convert myfile.cif myfile.pwi```\n",
    "\n",
    "If you run the simulation using ASE, this step is not needed since conversions are done internally and transparently \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the structure \n",
    "\n",
    "First we download a selected cif file from the Crystallography Open Database. (Of course, this could also be done manually.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CIF_file = '1527603.cif'\n",
    "print('Downloading CIF file '+CIF_file+' from crystallography.net')\n",
    "os.system(\"wget -c https://www.crystallography.net/cod/\"+CIF_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the pseudo potential for Nitrogen\n",
    "\n",
    "Quantum Espresso needs a suitable pseudopotential. We get it from their collection at "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pseudopotfile = 'N.pbe-n-kjpaw_psl.1.0.0.UPF'\n",
    "pseudo_dir = tmpdir+\"/pseudo/\"\n",
    "os.makedirs(pseudo_dir,exist_ok=True)\n",
    "os.chdir(pseudo_dir)\n",
    "os.system(\"wget -c https://www.quantum-espresso.org/upf_files/\"+pseudopotfile)\n",
    "os.system(\"wget -c https://raw.githubusercontent.com/PaNOSC-ViNYL/workshop2020/team2/demo/team2/N.pbe-n-radius_5.UPF\")\n",
    "os.chdir(tmpdir)\n",
    "pseudopotfile = 'N.pbe-n-radius_5.UPF'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the list of files in the current working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.listdir()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup the simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we read the CIF file and display the (initial) structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase import io, Atom, Atoms\n",
    "atomCIF = io.read(CIF_file)\n",
    "\n",
    "print(atomCIF)\n",
    "print(atomCIF.get_positions())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.visualize import view\n",
    "view(atomCIF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we setup the Quantum Espresso calculation. The parameters correspond to the ones found in the input files of pw.x, for which the documentation is available [here](https://www.quantum-espresso.org/Doc/INPUT_PW.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from ase.calculators.espresso import Espresso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pseudopotentials={'N': pseudopotfile}\n",
    "\n",
    "calc = Espresso(pseudopotentials=pseudopotentials,\n",
    "                tstress=True, \n",
    "                tprnfor=True, \n",
    "                kpts=(6, 6, 6),\n",
    "                ecutrho=480,\n",
    "                ecutwfc=60,\n",
    "                ibrav=0,\n",
    "                nat=8,\n",
    "                ntyp=1,\n",
    "                calculation='relax',\n",
    "                occupations='smearing',\n",
    "                smearing='cold',\n",
    "                degauss=0.001,\n",
    "                outdir=tmpdir,\n",
    "                pseudo_dir=pseudo_dir,\n",
    "                conv_thr=1e-7,\n",
    "                mixing_mode='plain',\n",
    "                electron_maxstep=80,\n",
    "                mixing_beta=0.5,\n",
    "                ion_dynamics='bfgs',\n",
    ")\n",
    "\n",
    "atom = atomCIF\n",
    "atom.calc = calc\n",
    "\n",
    "#atom.set_calculator(calc)\n",
    "#atom.get_potential_energy()\n",
    "#fermi_level = calc.get_fermi_level()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the potential energy\n",
    "\n",
    "This takes some time, about 15 minutes.\n",
    "When the calculation starts, two additional files are created:\n",
    " - espresso.pwi: QE input file with atomic structure and parameters\n",
    " - espresso.pwo: QE output \n",
    " \n",
    "ASE will parse the espresso.pwo file, and update the `atoms` object accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_energy = atom.get_potential_energy()\n",
    "\n",
    "print(\"Total energy: {0} eV\".format(potential_energy))\n",
    "print(\"Total energy: {0} eV\".format(atom.get_total_energy()))\n",
    "fermi_level = calc.get_fermi_level()\n",
    "print(\"Fermi energy: {0} eV\".format(fermi_level))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare atom positions before and after calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first read the output of the QE calculation, index=-1 allow to read only the last set of positions (those at convergence)\n",
    "atomsOUT = io.read('espresso.pwo',index=-1)\n",
    "atomsOUT.get_positions() - atom.get_positions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### output the result to CIF format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ase_outfile = 'output.cif'\n",
    "hklfile=ase_outfile+'.hkl'\n",
    "io.write(ase_outfile, atom)\n",
    "os.listdir()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system('cif2hkl '+ase_outfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.listdir()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neutron scattering experiment simulation\n",
    "\n",
    "The optimized structure is now ready. From this point, the data can be used in whatever way we desire. \n",
    "\n",
    "As an example, we use McStas to simulate a neutron-scattering experiment using the optimiized structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get McStasScript\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(\"git clone git@github.com:PaNOSC-ViNYL/McStas_ViNYL_concept.git\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir('McStas_ViNYL_concept')\n",
    "#os.chdir(tmpdir)\n",
    "os.makedirs(mcstas_outdir,exist_ok=True)\n",
    "os.listdir()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your McStas is installed in a different location, update the paths below accordingly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import McStasCalculator\n",
    "import McStasParameters\n",
    "import math\n",
    "from mcstasscript.interface import instr, plotter\n",
    "from mcstasscript.interface import functions\n",
    "my_configurator = functions.Configurator()\n",
    "my_configurator.set_mcrun_path(\"/usr/local/bin/\")\n",
    "my_configurator.set_mcstas_path(\"/usr/local/mcstas/2.6/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Instr = instr.McStas_instr(\"powder_diffractometer\")\n",
    "\n",
    "Instr.add_parameter(\"wavelength\", value=1.2, comment=\"[AA]\")\n",
    "\n",
    "src = Instr.add_component(\"Source\", \"Source_Maxwell_3\")\n",
    "src.xwidth = 0.12\n",
    "src.yheight = 0.12\n",
    "src.Lmin = \"wavelength*0.94\" # Simulate wavelengths in small band around requested wavelength\n",
    "src.Lmax = \"wavelength*1.06\"\n",
    "src.dist = 3.0\n",
    "src.focus_xw = guide_width = 0.04\n",
    "src.focus_yh = guide_height = 0.08\n",
    "\n",
    "# Set source spectrum to ILL\n",
    "src.T1 = 683.7\n",
    "src.I1 = 0.5874E13\n",
    "src.T2 = 257.7\n",
    "src.I2 = 2.5099E13\n",
    "src.T3 = 16.7\n",
    "src.I3 = 1.0343E12\n",
    "\n",
    "guide = Instr.add_component(\"guide\", \"Guide_gravity\", AT=[0,0,3.0], RELATIVE=\"Source\")\n",
    "guide.w1 = guide_width\n",
    "guide.h1 = guide_height\n",
    "guide.l = guide_length = 10 # 10 m long guide\n",
    "guide.m = 3.0\n",
    "guide.G = -9.82 # Gravity\n",
    "\n",
    "Instr.add_component(\"guide_end\", \"Arm\", AT=[0, 0, guide_length], RELATIVE=\"guide\")\n",
    "Instr.add_component(\"mono_pos\", \"Arm\", AT=[0, 0, 0.2], RELATIVE=\"guide_end\")\n",
    "\n",
    "Instr.add_parameter(\"mono_Q\", value=3.355, comment=\"Monochromator scattering vector length (PG) [AA^-1]\")\n",
    "Instr.add_declare_var(\"double\", \"mono_theta\")\n",
    "Instr.add_declare_var(\"double\", \"wavevector\")\n",
    "\n",
    "# Calculate wavevector and find theta from Q = 2k sin(theta)\n",
    "Instr.append_initialize(\"wavevector = 2*PI/wavelength;\")\n",
    "Instr.append_initialize(\"mono_theta = RAD2DEG*asin(0.5*mono_Q/wavevector);\")\n",
    "\n",
    "mono = Instr.add_component(\"mono\", \"Monochromator_curved\", AT=[0,0,0], RELATIVE=\"mono_pos\")\n",
    "mono.Q = \"mono_Q\"\n",
    "mono.height = 0.1\n",
    "mono.zwidth = 0.03\n",
    "mono.NH = 3 \n",
    "mono.NV = 11\n",
    "mono.RV = 1.5 # Focusing \n",
    "mono.set_ROTATED([0, \"mono_theta\", 0], RELATIVE=\"mono_pos\")\n",
    "\n",
    "Instr.add_component(\"mono_out\", \"Arm\", AT=[0,0,0], ROTATED=[0, \"mono_theta\", 0], RELATIVE=\"mono\")\n",
    "\n",
    "L_mon = Instr.add_component(\"L_mon\", \"L_monitor\", AT=[0, 0, 1.0], RELATIVE=\"mono_out\")\n",
    "L_mon.Lmin = \"wavelength*0.94\"\n",
    "L_mon.Lmax = \"wavelength*1.06\"\n",
    "L_mon.filename = '\"L_mon.dat\"'\n",
    "L_mon.xwidth = 0.1\n",
    "L_mon.yheight = 0.1\n",
    "L_mon.nL = 150\n",
    "\n",
    "sample = Instr.add_component(\"sample\", \"PowderN\", AT=[0, 0, 1.5], RELATIVE=\"mono_out\")\n",
    "sample.radius = 0.008\n",
    "sample.yheight = 0.03\n",
    "sample.reflections = '\"../output.cif.hkl\"'\n",
    "sample.barns = 0 # output.cif.hkl cross section read as fm^2 \n",
    "\n",
    "# Wish to focus on the detector, specify height and radius for use in focusing.\n",
    "detector_height = 0.2\n",
    "detector_radius = 1.0\n",
    "sample.d_phi = math.atan(detector_height/detector_radius)*180/3.14159\n",
    "sample.tth_sign = 1.0\n",
    "sample.set_SPLIT(1000)\n",
    "\n",
    "# Set up a banana monitor to measure scattering pattern\n",
    "monitor = Instr.add_component(\"monitor\", \"Monitor_nD\", RELATIVE=\"sample\")\n",
    "monitor.xwidth = 2.0*detector_radius\n",
    "monitor.yheight = detector_height\n",
    "monitor.options = '\"banana, theta limits=[10,170], bins=320\"'\n",
    "monitor.filename = '\"banana.dat\"'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture # the output of McStas is verbose and unneccesary for us.\n",
    "\n",
    "data = Instr.run_full_instrument(ncount=5E6, mpi=4,\n",
    "                                 foldername=mcstas_outdir, increment_folder_name=True)\n",
    "\n",
    "# Instr.show_instrument() Uncomment to view instrum\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotter.make_sub_plot(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}