{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise 4: Thermodynamic Properties" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Authors: Jan Janssen, Tilmann Hickel, Jörg Neugebauer ([Max-Planck-Institut für Eisenforschung](https://www.mpie.de))**\n", "\n", "The scope of this first exercise is to become familar with:\n", "* with the phonopy interface in pyiron to calculate free energies,\n", "* the harmonic and quasi-harmonic approximation and\n", "* how to combine multiple pyiron jobs in one workflow. \n", "\n", "With this last section we have all the necessary tools to calculate phase diagrams. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Update Notice \n", "Since the workshop in 2020, there have been some updates to pyiron commands, and here in the exercises new commands are used. \n", "\n", "The changes include:\n", "- Now pyiron has various modules, like pyiron_atomistics, pyiron_continuum, pyiron_experimental. So we can import `Project` from `pyiron_atomistics` directly.\n", "- using `pr.create` to create instances of various objects, such as structures and jobs. For example:\n", "```python\n", "pr.create.structure.ase.bulk() ## this replaces pr.create_ase_bulk()\n", "```\n", ", or\n", "```python\n", "pr.create.job.Lammps('job_name') ## this replaces pr.pr.create_job(job_type=pr.job_type.Lammps,'job_name')\n", "```\n", "\n", "Similary, one can create pyiron tables via `pr.create.table()`.\n", "\n", "By this approach, the user easily gets the available options via autocompeletion.\n", "\n", "**Please note that in the video tutorials, the old commands are still presented.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reminder\n", "In the first session we learned how to create a pyiron project object and then use this pyiron project object to create atomistic structure objects. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import the Project object\n", "from pyiron_atomistics import ____" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a Project object instance for a project named phonons\n", "pr = ____(\"phonons\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a cubic aluminum fcc structure\n", "al_fcc = pr.create.structure.ase.bulk(____, _____=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Confirm the final structure has 4 atoms by calculating \n", "# the length of the structure object\n", "____(al_fcc) == 4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Harmonic Approximation \n", "Calculate the phonons at the 0K equilibrium volume and then use the model of the harmonic oscillator to calculate the free energy, heat capacity and entropy contribution. These calculation can again be either executed with density functional theory (DFT) or interatomic potentials. In this workshop we use interatomic potentials to accelerate the process, but the same simulation protocol can also be executed with DFT code as quantum engine. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Minimization" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a LAMMPS job to relax the structure - here we the regular minimization\n", "# for DFT calculation it is recommended to calculate the energy volume curve \n", "# to identify the equilibrium volume as discussed in the previous section. \n", "# Name the minimization job \"lmp_mini\".\n", "job_mini = pr.create.job.Lammps(\n", " job_name=________\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "potential = \"2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1\"\n", "job_mini.potential = potential" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the the cubic fcc aluminium structure to the LAMMPS quantum engine job\n", "job_mini._________ = _________" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set the pressure during the minimization to zero \n", "job_mini.calc_minimize(_________=_______)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the minimization\n", "job_mini.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Phonon Calculation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a interatomic template job with the LAMMPS quantum engine named lmp\n", "job_lammp_template = pr.create.job.Lammps(\n", " job_name=________\n", ")\n", "job_lammp_template.potential = potential" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the the final structure of the minimization calculation job as input structure\n", "job_lammp_template._________ = _________.get_structure()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# We calculate the bulk energy using the template job in a separate calculation named \"lmp_bulk\"\n", "job_bulk = job_lammp_template.copy_to(\n", " new_job_name=___________, \n", " new_database_entry=False\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a Phonopy job from the LAMMPS quantum engine named phono\n", "phono = pr.create.job.PhonopyJob(\n", " job_name=________\n", ")\n", "phono.ref_job = job_lammp_template" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the Phonopy calculation and the bulk calculation\n", "job_bulk.run()\n", "phono.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Density of States" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the density of states over energy \n", "phono.plot_dos();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Free energy calculation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# To calculate the thermodynamic properties we define a temperature range\n", "# Starting a 0K up to the 800K using 41 steps. We import numpy and define\n", "# a linear space starting at 0 to 800 with 41 steps. \n", "import numpy as np\n", "temperatures = np.linspace(______, _______, _____)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate the thermal properties for a the defined temperature range \n", "bulk_thermal_properties = phono.get_thermal_properties(temperatures=_________) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import the matplotlib library for plotting. \n", "import matplotlib.pyplot as plt\n", "\n", "# Plot the free energy over temperature by adding the inner bulk energy and the free vibrational energy \n", "plt.plot(__________, job_bulk.output.energy_pot[-1] + bulk_thermal_properties.free_energies)\n", "plt.xlabel(\"Temperature [K]\")\n", "plt.ylabel(\"Free energy ($U+F_{vib}$) [eV]\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Heat capacity" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the heat capacity over temperature\n", "plt.plot(___________, bulk_thermal_properties.cv)\n", "plt.xlabel(\"Temperature [K]\")\n", "plt.ylabel(\"Heat capacity\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Entropy" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the entropy over temperature\n", "plt.plot(___________, bulk_thermal_properties.entropy)\n", "plt.xlabel(\"Temperature [K]\")\n", "plt.ylabel(\"Entropy\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quasi-harmonic Approximation\n", "The quasi-harmonic approximation includes the thermal volume expansion in contrast to the harmonic approximation which only considers the 0K equilibrium volume. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a interatomic template job with the LAMMPS quantum engine named lmp_strain\n", "job_strain_template = pr.create.job.Lammps(\n", " job_name=________\n", ")\n", "job_strain_template.potential = potential" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the the final structure of the minimization calculation job as input structure\n", "job_strain_template.structure = ________.get_structure()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a secondary template job to calculate the bulk energy\n", "# The second template job we name lmp_bulk_strain \n", "job_strain_bulk = job_strain_template.copy_to(\n", " new_job_name=______________, \n", " new_database_entry=False\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Use this second template job to create a Murnaghan job named murn_strain\n", "murn_strain = pr.create.job.Murnaghan(\n", " job_name=____________\n", ")\n", "murn_strain.ref_job = job_strain_bulk" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a Phonopy job from the original LAMMPS quantum engine named phono_strain\n", "phono_strain = pr.create.job.PhonopyJob(\n", " job_name=____________\n", ")\n", "phono_strain.ref_job = job_strain_template" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign this Phonopy job to a QuasiHarmonicJob to calculate the phonons at multiple volumes\n", "# This QuasiHarmonicJob is named quasi_strain\n", "quasi_strain = pr.create.job.QuasiHarmonicJob(\n", " job_name=_______________\n", ")\n", "quasi_strain.ref_job = phono_strain" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set the end temperature to be 800K and set the number of temperature steps to 41 \n", "quasi_strain.input[\"temperature_end\"] = ________\n", "quasi_strain.input[\"temperature_steps\"] = _______" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the QuasiHarmonicJob calculation and the bulk calculation\n", "murn_strain.run()\n", "quasi_strain.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Free Energy\n", "To calculate the heat capacity and entropy at constant pressure, the free energy at constant temperature is plotted over volume. The minimum at constant temperature defines the volume of constant pressure. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Visualise the temperature dependent free energy over volume using a matplotlib color map\n", "import matplotlib\n", "cmap = matplotlib.cm.get_cmap('coolwarm')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Iterate over the the output of the QuasiHarmonicJob to plot the temperature dependent free energy over volume \n", "for i, [t, free_energy, v] in enumerate(\n", " zip(____________[\"output/temperatures\"].T, \n", " ____________[\"output/free_energy\"].T, \n", " ____________[\"output/volumes\"].T)):\n", " color = cmap(i/len(quasi_strain[\"output/temperatures\"].T))\n", " # Add the energy of the Murnaghan Job to the vibrational free energy\n", " plt.plot(v, free_energy + _________[\"output/energy\"], color=color)\n", "\n", "# Add labels to the plot \n", "plt.xlabel(\"Volume\")\n", "plt.ylabel(\"Free Energy\")\n", "\n", "# Add a color bar to visualise the temperature dependence \n", "temperatures = quasi_strain[\"output/temperatures\"]\n", "normalize = matplotlib.colors.Normalize(vmin=temperatures.min(), vmax=temperatures.max())\n", "scalarmappaple = matplotlib.cm.ScalarMappable(norm=normalize, cmap=cmap)\n", "scalarmappaple.set_array(temperatures)\n", "cbar = plt.colorbar(scalarmappaple)\n", "cbar.set_label(\"Temperature\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pressure optimised thermodynamic properties\n", "After calculating the volume expansion from the free energy surface over temperature and volume, the entropy and heat capacity are calculated at these optimised volumes." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Use the optimise_volume() function of the QuasiHarmonicJob \n", "v0_lst, free_eng_lst, entropy_lst, cv_lst = quasi_strain.optimise_volume(\n", " # It requires the output energy of the energy volume curve as additional input \n", " bulk_eng=murn_strain[\"output/energy\"]\n", ")\n", "temperature_output_lst = quasi_strain[\"output/temperatures\"][0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the finite temperature volume over temperature \n", "plt.plot(temperature_output_lst, __________)\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Volume\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the pressure optimised free energy over temperature \n", "plt.plot(temperature_output_lst, ___________)\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Free Energy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the pressure optimised entropy over temperature \n", "plt.plot(temperature_output_lst, ____________)\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Entropy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the pressure optimised heat capacity over temperature \n", "plt.plot(temperature_output_lst, ____________)\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Heat Capacity\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Concentration dependence \n", "With the temperature dependence of a unary discussed above the next step is the concentration dependence. For the concentration dependence we use the special quasi-random structures introduced in the first section to calculate mixed structures. Here we select the Ni-Cr phase diagram, which can be calculated using the `2018--Howells-C-A--Cr-Ni--LAMMPS--ipr1` potential. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create Endmembers\n", "Start by creating the endmembers of both phases the nickel fcc endmember and the chromium bcc endmember." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create the Nickel fcc endmember \n", "ni_fcc = pr.create.structure.ase.bulk(____, cubic=True)\n", "\n", "# Create the Chromium bcc endmember \n", "cr_bcc_small = pr.create.structure.ase.bulk(____, cubic=True)\n", "\n", "# Compare the length of both structures\n", "len(ni_fcc), len(cr_bcc_small)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Repeat the bcc cell to have the same number of atoms as the fcc cell\n", "cr_bcc = cr_bcc_small.repeat([___, ___, ____])\n", "\n", "# Compare the length of both structure\n", "len(cr_bcc) == len(ni_fcc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Select concentrations \n", "One advantage of atomistic simulation is the ability to calculate free energies of unstable phases. In this case we calculate both the fcc phase and the bcc phase for all temperature ranges. Starting at 0% Cr up to 100%. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create 3 mixed concentrations for a 4 atom cell within the range 0.0 to 1.0 \n", "# for example these could be 0.25, 0.5 and 0.75\n", "concentration_lst = [____, ____, ___]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an FCC Cr lattice by replacing all elements of the FCC Ni lattice with Cr\n", "cr_fcc = ______.copy()\n", "cr_fcc[:] = \"Cr\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an BCC Ni lattice by replacing all elements of the BCC Cr lattice with Ni\n", "ni_bcc = ________.copy()\n", "ni_bcc[:] = \"Ni\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a list of all concentrations \n", "concentration_total_lst = [0.0] + concentration_lst + [1.0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SQS Structures\n", "To calculate SQS structures for multiple concentrations the SQSMaster job is used, it takes an SQS job as an input in addtion to a list of concentrations. Here it is used to calculate SQS structures for both the FCC phase and the BCC phase." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### FCC" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an SQS job named sqs_job_ni\n", "sqs_job_ni = pr.create.job.SQSJob(\n", " job_name=__________\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the cubic nickel fcc structure to the SQS job\n", "sqs_job_ni.structure = ______" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Limit the number of iterations to 1000\n", "sqs_job_ni.input['iterations'] = _____" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an SQSMaster named sqs_master_ni\n", "master_ni = pr.create.job.SQSMaster( \n", " job_name=__________\n", ")\n", "master_ni.ref_job = sqs_job_ni" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the mixed concentration list defined above - not the total concentration list which includes the end members\n", "master_ni.input[\"fraction_lst\"] = __________\n", "master_ni.input[\"species_one\"] = \"Ni\"\n", "master_ni.input[\"species_two\"] = \"Cr\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the calculation\n", "master_ni.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### BCC " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an SQS job named sqs_job_cr\n", "sqs_job_cr = pr.create.job.SQSJob(\n", " job_name=_________\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the cubic chromium bcc structure to the SQS job\n", "sqs_job_cr.structure = __________" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Limit the number of iterations to 1000\n", "sqs_job_cr.input['iterations'] = _______" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create an SQSMaster named sqs_master_cr\n", "master_cr = pr.create.job.SQSMaster( \n", " job_name=_________\n", ")\n", "master_cr.ref_job = sqs_job_cr" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Assign the mixed concentration list defined above - not the total concentration list which includes the end members\n", "master_cr.input[\"fraction_lst\"] = _________\n", "master_cr.input[\"species_one\"] = \"Ni\"\n", "master_cr.input[\"species_two\"] = \"Cr\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Execute the calculation\n", "master_cr.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Collect structures\n", "At this point ten structures have been generated five for each phase. Now each of these structures has to be minimized to apply the quasi-harmonic approximation for each of them. In the interest of time this step is skipped here and we only compile the list of structures." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Combine the end members with the structures from the SQS masters \n", "# to have a total of five structures for each phase\n", "structure_fcc_lst = [______] + master_ni.list_structures() + [______]\n", "structure_bcc_lst = [______] + master_cr.list_structures() + [______]\n", "\n", "len(structure_fcc_lst) == len(structure_bcc_lst)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "In this section you learned: \n", "* how to calculate free energies using the harmonic and the quasi-harmonic approximation. \n", "* how hierachical workflows can be constructed in pyiron by combining multiple master jobs. \n", "* and how to generate mixed phases with special quasi random structures. \n", "\n", "Suggestions: \n", "* Calculate the concentration dependent volume expansion for both the FCC and the BCC phase. \n", "* Calculate the harmonic approximation for a single SQS structure. As the mixed phase breaks the symmetry additional displacments are necessary. Phonopy automatically identifies the necessary displacements, but these calculation take more time. \n", "* Use the option `job.server.run_mode.interactive=True` for the LAMMPS jobs before they are assigned to the phonopy jobs to accelerate the calculations by using LAMMPS as a C-library rather than calling the executable by writing output and input files. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "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.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }