{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solution\n", "\n", "Here is the solution to the [first exercise](https://mmoriniere.gitlab.io/MyBigDFT/notebooks/Exercise_01_Polarizability_Tensor_jobs.html) regarding the polarizability tensor." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the ground state calculation\n", "\n", "You will first have to define the initial geometry and the input parameters given the following indications.\n", "\n", "We want to study the N$_2$ molecule, starting from the optimized geometry ($d_{N-N} = 1.0935 \\unicode[serif]{xC5}$) found using the LDA exchange-correlation function (which is the BigDFT default). We will use free boundary conditions.\n", "\n", "Regarding the wavelet grid, we want to use a fine grid multiplying radius of 7 and a coarse grid multiplying radius of 9 (`rmult = [7, 9]`) and grid spacings of 0.35 (`hgrids = 0.35`). Hint: the above parameters are stored under the `\"dft\"` key." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from mybigdft import Posinp, Atom, InputParams\n", "\n", "# Let us start with the initial geometry\n", "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0935])]\n", "gs_pos = Posinp(atoms, units=\"angstroem\", boundary_conditions=\"free\")\n", "\n", "# Let us continue with the input parameters\n", "gs_inp = InputParams({\"dft\": {\"rmult\": [7, 9], \"hgrids\": 0.35}})" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Let us ensure that you followed the instructions:\n", "# - The system is made of two nitrogen atoms\n", "assert len(gs_pos) == 2 and all([atom.type == \"N\" for atom in gs_pos])\n", "# - correct distance between the atoms\n", "assert gs_pos.distance(0, 1) == 1.0935 and gs_pos.units == \"angstroem\"\n", "# - correct boundary conditions\n", "assert gs_pos.boundary_conditions == \"free\"\n", "# - correct rmult\n", "assert gs_inp[\"dft\"][\"hgrids\"] == 0.35\n", "# - correct hgrids\n", "assert gs_inp[\"dft\"][\"rmult\"] == [7, 9]\n", "# - no external electric field is applied\n", "assert \"elecfield\" not in gs_inp[\"dft\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can now define a `Job` instance and run it. It has to be run from within a context manager (that is, using the `with` statement).\n", "\n", "In the cell below, you must use the above-defined input parameters and positions to initialize a `Job` instance. The job must also be named `\"N2\"` and run in the `\"N2/pol_tensor/\"` directory.\n", "\n", "You might also want to run the calculations in parallel using mpirun (set the number of processors via the `nmpi` optional argument of the `run` method) or use multiple threads via OpenMP (set the number of threads via the `nomp` optional argument of the `run` method)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "import os\n", "from mybigdft import Job\n", "\n", "# Define the job and run it, using the above indications\n", "base_run_dir = os.path.relpath(\"N2/pol_tensor\")\n", "with Job(posinp=gs_pos, inputparams=gs_inp, name=\"N2\",\n", " run_dir=base_run_dir) as gs_job:\n", " gs_job.run(nmpi=6, nomp=3)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Let us ensure that you followed the instructions:\n", "# - correct posinp\n", "assert gs_job.posinp == gs_pos == gs_job.logfile.posinp\n", "# - correct input parameters\n", "assert gs_job.inputparams == gs_inp == gs_job.logfile.inputparams\n", "# - correct job name\n", "assert gs_job.name == \"N2\"\n", "# - correct run directory\n", "assert gs_job.run_dir.endswith(\"N2/pol_tensor\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run a job with an electric field in the $x$ direction\n", "\n", "Using the same scheme, you can now define another another job, where an external electric field (of amplitude `1.e-4` Ha.bohr$^{-1}$) is applied along the $x$ direction. You may keep the same name for the job, but run it in another directory, namely `\"N2/pol_tensor/EF_along_x+/\"`.\n", "\n", "Keep the same input geometry and add the key `\"elecfield\"` to the input parameters parameters. The associated value has to be a list of length 3, containing the amplitude of the electric field along the three directions of space." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "# Define the new input parameters\n", "EF_x_inp = InputParams({\"dft\": {\"rmult\": [7, 9], \"hgrids\": 0.35,\n", " \"elecfield\": [1.e-4, 0.0, 0.0]}})\n", "\n", "# Define the job and run it, using the above indications\n", "with Job(posinp=gs_pos, inputparams=EF_x_inp, name=\"N2\",\n", " run_dir=os.path.join(base_run_dir, \"EF_along_x+\")) as EF_x_job:\n", " EF_x_job.run(nmpi=6, nomp=3)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Let us ensure that you followed the instructions:\n", "# - same basic input parameters\n", "assert EF_x_inp[\"dft\"][\"rmult\"] == gs_inp[\"dft\"][\"rmult\"]\n", "assert EF_x_inp[\"dft\"][\"hgrids\"] == gs_inp[\"dft\"][\"hgrids\"]\n", "# - correct electric field\n", "assert EF_x_job.logfile.inputparams[\"dft\"][\"elecfield\"] == [1.e-4, 0.0, 0.0]\n", "# - correct job name\n", "assert EF_x_job.name == \"N2\"\n", "# - correct run directory\n", "assert EF_x_job.run_dir.endswith(\"N2/pol_tensor/EF_along_x+\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Start building the polarizability tensor\n", "\n", "Given the formula in the beginning of this notebook, you can readily compute three elements of the polarizability tensor: $\\alpha_{xx}$, $\\alpha_{yx}$ and $\\alpha_{zx}$. They can be found computing the difference of the dipole vector found with and without an applied electric field and divide by the amplitude of the electric field. Note that you can access the dipole of a from the logfile of a job via its `dipole` attribute (*i.e.*, via `job.logfile.dipole`). You might want to use the [numpy](https://docs.scipy.org/doc/numpy/user/quickstart.html) module to build the polarizability tensor." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.075753e+01 0.000000e+00 0.000000e+00]\n", " [ 5.400000e-04 0.000000e+00 0.000000e+00]\n", " [-1.740000e-03 0.000000e+00 0.000000e+00]]\n" ] } ], "source": [ "import numpy as np\n", "\n", "pol_tensor = np.zeros((3, 3))\n", "# Update the polarizability tensor\n", "d_0 = np.array(gs_job.logfile.dipole)\n", "d_x = np.array(EF_x_job.logfile.dipole)\n", "delta_ef = 1.e-4\n", "pol_x_column = (d_x - d_0) / delta_ef\n", "pol_tensor[:, 0] = pol_x_column\n", "print(pol_tensor)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "expected_pol_tensor = [[ 1.075753e+01, 0.0, 0.0],\n", " [ 5.400000e-04, 0.0, 0.0],\n", " [-1.740000e-03, 0.0, 0.0]]\n", "assert np.allclose(pol_tensor, expected_pol_tensor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Apply an electric field in the $y$ and $z$ direction\n", "\n", "The first column of the polarizability tensor now being computed, it should be easy to define and run two more jobs: one with the electric field along the $y$ axis, the other with an electric field along the $z$ axis. Do not forget to change the name of the run directory for each job." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "# Define and run a job with an electric field along the the y direction\n", "EF_y_inp = InputParams({\"dft\": {\"rmult\": [7, 9], \"hgrids\": 0.35,\n", " \"elecfield\": [0.0, 1.e-4, 0.0]}})\n", "with Job(posinp=gs_pos, inputparams=EF_y_inp, name=\"N2\",\n", " run_dir=os.path.join(base_run_dir, \"EF_along_y+\")) as EF_y_job:\n", " EF_y_job.run(nmpi=6, nomp=3)\n", "\n", "# Define and run a job with an electric field along the the z direction\n", "EF_z_inp = InputParams({\"dft\": {\"rmult\": [7, 9], \"hgrids\": 0.35,\n", " \"elecfield\": [0.0, 0.0, 1.e-4]}})\n", "with Job(posinp=gs_pos, inputparams=EF_z_inp, name=\"N2\",\n", " run_dir=os.path.join(base_run_dir, \"EF_along_z+\")) as EF_z_job:\n", " EF_z_job.run(nmpi=6, nomp=3)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Let us ensure that you followed the instructions:\n", "# - correct electric fields applied\n", "assert EF_y_job.logfile.inputparams[\"dft\"][\"elecfield\"] == [0.0, 1.e-4, 0.0]\n", "assert EF_z_job.logfile.inputparams[\"dft\"][\"elecfield\"] == [0.0, 0.0, 1.e-4]\n", "# - correct run directories\n", "assert EF_y_job.run_dir.endswith(\"N2/pol_tensor/EF_along_y+\")\n", "assert EF_z_job.run_dir.endswith(\"N2/pol_tensor/EF_along_z+\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Complete the computation of the polarizability tensor\n", "\n", "You now are able to compute the rest of the elements of the polarizability tensor, using the same approach as for the first one. You may use a `for` loop running over the last two jobs (`EF_y_job` and `EF_z_job`) to do so." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1.075753e+01 5.400000e-04 -9.900000e-04]\n", " [ 5.400000e-04 1.075753e+01 -9.900000e-04]\n", " [-1.740000e-03 -1.740000e-03 1.525630e+01]]\n" ] } ], "source": [ "# Define the rest of the polarizability tensor from the logfiles\n", "for i, job in enumerate([EF_y_job, EF_z_job]):\n", " d_i = np.array(job.logfile.dipole)\n", " pol_tensor[:, i+1] = (d_i - d_0) / delta_ef\n", "print(pol_tensor)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Let us ensure that you obtained the expected result\n", "expected_pol_tensor = [[ 1.075753e+01, 5.400000e-04, -9.900000e-04],\n", " [ 5.400000e-04, 1.075753e+01, -9.900000e-04],\n", " [-1.740000e-03, -1.740000e-03, 1.525630e+01]]\n", "assert np.allclose(pol_tensor, expected_pol_tensor)" ] } ], "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.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }