{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solution\n", "\n", "Here is the solution to the [third exercise](https://mmoriniere.gitlab.io/MyBigDFT/notebooks/Exercise_03_Polarizability_Tensor_Workflow_2.html) regarding the polarizability tensor." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a class\n", "\n", "The use of the `Workflow` class in the previous notebook showed you that having a clear separation of the initialization and post-processing of a workflow is important, notably to reduce the amount of code to copy and paste each time you want to change the system or the input parameters. We will go even further in that direction by writing a class deriving from the `AbstractWorkflow` class. You will see that the code you implemented in the previous exercise can be quickly adapted to meet our needs here.\n", "\n", "Part of the body of the class is already written, you must only code the following:\n", "\n", "- add the string `\"pol_tensor\"` to the POST_PROCESSING_ATTRIBUTES list,\n", "- initialize the queue of jobs in the `__init__` method,\n", "- define the `post_proc` method to compute the value of the `pol_tensor` attribute.\n", "\n", "The `__init__` method has `ground_state` as argument, which is meant to be a `Job` instance of the reference ground state. This job must be the first one in the queue. The three other jobs must be added by taking this one job as reference, only modifying the run directory (in the same manner as in the previous exercise) and the input parameters, so as to apply an electric field in each direction ($x$, $y$ and $z$). To that end, you must add an optional argument `ef_amplitude` to the `__init__` method to state the value of the electric field amplitude to be used in all the jobs (default value: `1.e-4`).\n", "\n", "The definition of the `pol_tensor` method is given so as to give you an example of how to actually code such methods. Note the use of a `_poltensor` attribute in that method and in the `post_proc` method; this is what we meant by \"semi-private attributes\" in the previous section: defining a property returning such a semi-private attribute ensures that you can access the `pol_tensor` attribute here, but it cannot be set. The only way of modifying the value of this attribute is to change the value of the `_poltensor`. This is actually what is done in the `post_proc` method!" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "from copy import deepcopy\n", "import numpy as np\n", "from mybigdft import Job\n", "from mybigdft.workflows.workflow import AbstractWorkflow\n", "\n", "class PolTensor(AbstractWorkflow):\n", " \"\"\"\n", " This workflow allows to compute the (electronic) polarizability\n", " tensor of a given system.\n", "\n", " The polarizability tensor represents the response of the charges of\n", " a system (its dipole) to the application of an external electric\n", " field.\n", "\n", " To compute this polarizability tensor, some BigDFT calculations are\n", " performed, where the system is subject to an external electric\n", " field along each direction of space (:math:`x`, :math:`y` and\n", " :math:`z`). The elements of the polarizability tensor are then\n", " defined by the ratio of the delta of the dipole in one direction and\n", " the delta of the electric field amplitudes:\n", "\n", " .. math::\n", "\n", " \\alpha_{ij} = \\frac{\\Delta D_i}{\\Delta E_j}\n", "\n", " where :math:`i, j \\in \\{x, y, z\\}`, :math:`\\Delta D_i` is the\n", " variation of the dipole along the :math:`i` direction and\n", " :math:`\\Delta E_j` is the variation of the electric field amplitude\n", " along the :math:`j` direction.\n", " \"\"\"\n", " \n", " POST_PROCESSING_ATTRIBUTES = [\"pol_tensor\"]\n", " \n", " def __init__(self, ground_state, ef_amplitude=1.e-4):\n", " \"\"\"\n", " A PolTensor workflow is initialized by the job of the\n", " ground-state of the system and an electric field amplitude.\n", "\n", " Parameters\n", " ----------\n", " ground_state: Job\n", " Reference job, specifying the base input parameters (where no\n", " electric field is applied), the geometry of the system and the\n", " base run directory of the calculations.\n", " ef_amplitude: float\n", " Electric field amplitude to be applied along each direction,\n", " in atomic units (default to 1.e-4).\n", " \"\"\"\n", " self._ground_state = ground_state\n", " # Set the queue of jobs thanks to ground_state and ef_amplitude \n", " pt_queue = [self.ground_state]\n", " for i, coord in enumerate([\"x\", \"y\", \"z\"]):\n", " # Add the electric field to the base input parameters\n", " elecfield = [0.0] * 3\n", " elecfield[i] = ef_amplitude\n", " new_inp = deepcopy(inp)\n", " new_inp[\"dft\"][\"elecfield\"] = elecfield\n", " # Define the correct run directory\n", " new_run_dir = os.path.join(ground_state.run_dir,\n", " \"EF_along_{}+/\".format(coord))\n", " # Append the new job to the end of the queue\n", " pt_queue.append(Job(inputparams=new_inp, posinp=pos,\n", " name=\"N2\", run_dir=new_run_dir))\n", " # This calls the __init__ method of the AbstractWorkflow class to\n", " # set the queue, the run method (that calls the post_proc method\n", " # automatically) and the default value of the post-processing\n", " # attributes to None.\n", " super().__init__(queue=pt_queue)\n", " \n", " @property\n", " def ground_state(self):\n", " \"\"\"\n", " Set ground_state as a semi-private attribute: it can be accessed,\n", " but it cannot be set.\n", " \n", " Returns\n", " -------\n", " Job\n", " Reference job, where no electric field is applied.\n", " \"\"\"\n", " return self._ground_state\n", " \n", " @property\n", " def pol_tensor(self):\n", " \"\"\"\n", " Set pol_tensor as a semi-private attribute: it can be accessed,\n", " but it cannot be set.\n", " \n", " Returns\n", " -------\n", " numpy array of shape (3, 3)\n", " Polarizability tensor of the system, using :attr:`ground_state`\n", " as a reference job.\n", " \"\"\"\n", " return self._pol_tensor\n", " \n", " @property\n", " def mean_polarizability(self):\n", " \"\"\"\n", " Set mean_polarizability as a semi-private attribute: it can be \n", " accessed, but it cannot be set.\n", " \n", " Returns\n", " -------\n", " numpy array of shape (3, 3)\n", " Polarizability tensor of the system, using :attr:`ground_state`\n", " as a reference job.\n", " \"\"\"\n", " return self._mean_polarizability\n", " \n", " def post_proc(self):\n", " \"\"\"\n", " Compute and set the polarizability tensor and the mean polarizability\n", " of the system.\n", " \"\"\"\n", " # Polarizability tensor\n", " pol_tensor = np.zeros((3, 3))\n", " delta_ef = max(self.queue[1].inputparams[\"dft\"][\"elecfield\"])\n", " d0 = np.array(self.queue[0].logfile.dipole)\n", " for i, job in enumerate(self.queue[1:4]):\n", " d1 = np.array(job.logfile.dipole)\n", " pol_tensor[:, i] = (d1 - d0) / delta_ef\n", " self._pol_tensor = pol_tensor\n", " # Mean polarizability\n", " self._mean_polarizability = np.trace(self.pol_tensor) / 3." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Correct error catched!\n", "Correct error catched!\n" ] } ], "source": [ "# Let us test that class with for the N2 molecule\n", "from mybigdft import Job, Posinp, Atom, InputParams\n", "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0935])]\n", "pos = Posinp(atoms, units=\"angstroem\", boundary_conditions=\"free\")\n", "inp = InputParams({\"dft\": {\"rmult\": [7, 9], \"hgrids\": 0.35}})\n", "\n", "# 1- Using non-default electric field amplitude:\n", "gs_1 = Job(posinp=pos, inputparams=inp, name=\"N2\",\n", " run_dir=\"N2/pol_tensor/non_default/\")\n", "pt_wf_1 = PolTensor(gs_1, ef_amplitude=5.e-4)\n", "# - correct electric field\n", "expected_efs_1 = [None, [5.e-4, 0, 0], [0, 5.e-4, 0], [0, 0, 5.e-4]]\n", "for job, ef in zip(pt_wf_1.queue, expected_efs_1):\n", " assert job.inputparams[\"dft\"].get(\"elecfield\") == ef\n", "# - correct run directories\n", "expected_run_dirs_1 = [\"non_default/\", \"non_default/EF_along_x+\",\n", " \"non_default/EF_along_y+\", \"non_default/EF_along_z+\"]\n", "for job, run_dir in zip(pt_wf_1.queue, expected_run_dirs_1):\n", " assert job.run_dir.endswith(run_dir)\n", "# - correct pol_tensor attribute\n", "assert pt_wf_1.pol_tensor is None\n", "try:\n", " pt_wf_1.pol_tensor = np.zeros((3, 3))\n", "except AttributeError as e:\n", " assert \"\"\"AttributeError(\"can't set attribute\"\"\" in repr(e)\n", " print(\"Correct error catched!\")\n", " \n", "# 2- Using default electric field amplitude:\n", "gs_2 = Job(posinp=pos, inputparams=inp, name=\"N2\",\n", " run_dir=\"N2/pol_tensor/default/\")\n", "pt_wf_2 = PolTensor(gs_2)\n", "# - correct electric field\n", "expected_efs_2 = [None, [1.e-4, 0, 0], [0, 1.e-4, 0], [0, 0, 1.e-4]]\n", "for job, ef in zip(pt_wf_2.queue, expected_efs_2):\n", " assert job.inputparams[\"dft\"].get(\"elecfield\") == ef\n", "# - correct run directories\n", "expected_run_dirs_2 = [\"default/\", \"default/EF_along_x+\",\n", " \"default/EF_along_y+\", \"default/EF_along_z+\"]\n", "for job, run_dir in zip(pt_wf_2.queue, expected_run_dirs_2):\n", " assert job.run_dir.endswith(run_dir)\n", "# - correct pol_tensor attribute\n", "assert pt_wf_2.pol_tensor is None\n", "try:\n", " pt_wf_2.pol_tensor = np.zeros((3, 3))\n", "except AttributeError as e:\n", " assert \"\"\"AttributeError(\"can't set attribute\"\"\" in repr(e)\n", " print(\"Correct error catched!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Use that class!\n", "\n", "The `PolTensor` class now being implemented, you can use it in order to (effortlessly) study some problems. Given the implementation of that class, you can solve some problems in the most concise manner, especially compared to what was possible if you were to use the schemes presented in the previous notebooks." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* ### The polarizability tensor must not depend on the rotation of the initial molecule\n", "\n", "The atoms of the molecule here lie along the $z$ axis and the results should not change if the atoms lie along the $x$ or $y$ axis, for instance. Use the default electric field amplitude." ] }, { "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/atoms_along_x\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "[[ 1.116730e+01 4.103100e-01 4.087800e-01]\n", " [ 5.400000e-04 1.075753e+01 -9.900000e-04]\n", " [-4.115100e-01 -4.115100e-01 1.484653e+01]]\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "[[ 1.116730e+01 4.103100e-01 4.087800e-01]\n", " [ 5.400000e-04 1.075753e+01 -9.900000e-04]\n", " [-4.115100e-01 -4.115100e-01 1.484653e+01]]\n" ] } ], "source": [ "# Compute the polarizability tensor when the molecule is along the x axis\n", "atoms_x = [Atom('N', [0, 0, 0]), Atom('N', [1.0935, 0, 0])]\n", "pos_x = Posinp(atoms_x, units=\"angstroem\", boundary_conditions=\"free\")\n", "run_dir_x = \"N2/pol_tensor/atoms_along_x/\"\n", "gs_x = Job(inputparams=inp, posinp=pos_x, name=\"N2\", run_dir=run_dir_x)\n", "pt_wf_x = PolTensor(gs_x)\n", "pt_wf_x.run(nmpi=6, nomp=3)\n", "print(pt_wf_x.pol_tensor)\n", "\n", "# Compute the polarizability tensor when the molecule is along the x axis\n", "atoms_y = [Atom('N', [0, 0, 0]), Atom('N', [0, 1.0935, 0])]\n", "pos_y = Posinp(atoms_y, units=\"angstroem\", boundary_conditions=\"free\")\n", "run_dir_y = \"N2/pol_tensor/atoms_along_y/\"\n", "gs_y = Job(inputparams=inp, posinp=pos_y, name=\"N2\", run_dir=run_dir_y)\n", "pt_wf_y = PolTensor(gs_x)\n", "pt_wf_y.run(nmpi=6, nomp=3)\n", "print(pt_wf_y.pol_tensor)" ] }, { "cell_type": "code", "execution_count": 4, "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", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/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", "[[ 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": [ "# Let us ensure that you obtained the expected results:\n", "# - both tensors are the same\n", "assert np.allclose(pt_wf_x.pol_tensor, pt_wf_y.pol_tensor)\n", "# - they have the same trace as the polarizability tensor where the molecule\n", "# is along the z axis\n", "gs_z = Job(posinp=pos, inputparams=inp, name=\"N2\", run_dir=\"N2/pol_tensor/\")\n", "pt_wf_z = PolTensor(gs_z)\n", "pt_wf_z.run(nmpi=6, nomp=3)\n", "print(pt_wf_z.pol_tensor)\n", "assert np.isclose(np.trace(pt_wf_x.pol_tensor), np.trace(pt_wf_z.pol_tensor))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the polarizability tensors are note the same if the molecule lies along the $z$ axis or not. The mean polarizability is the same, though." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* ### What is the influence of the electric field amplitude on the results?\n", "\n", "Compute the polarizability tensor of the N$_2$ molecule using various electric field amplitudes and comment on how this affects the value of polarizability tensor. Keep the same structure and input parameters as above." ] }, { "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_amplitude_1e-11\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-11/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-11/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-11/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-10\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-10/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-10/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-10/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-09\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-09/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-09/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-09/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-08\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-08/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-08/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-08/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-07\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-07/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-07/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-07/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-06\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-06/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-06/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-06/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-05\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-05/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-05/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_1e-05/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.0001\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.0001/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.0001/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.0001/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.001\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.001/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.001/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.001/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.01\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.01/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.01/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.01/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.1\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.1/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.1/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/ef_amplitude_0.1/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "Electric field amplitude: 1e-11 a.u.\n", "[[0. 0. 0.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "\n", "Electric field amplitude: 1e-10 a.u.\n", "[[10. 0. 0.]\n", " [ 0. 10. 0.]\n", " [ 0. 0. 10.]]\n", "\n", "Electric field amplitude: 1e-09 a.u.\n", "[[11. 0. 0.]\n", " [ 0. 11. 0.]\n", " [ 0. 0. 15.]]\n", "\n", "Electric field amplitude: 1e-08 a.u.\n", "[[10.8 0. 0. ]\n", " [ 0. 10.8 0. ]\n", " [ 0. 0. 15.2]]\n", "\n", "Electric field amplitude: 1e-07 a.u.\n", "[[10.76 0. 0. ]\n", " [ 0. 10.76 0. ]\n", " [ 0. 0. 15.25]]\n", "\n", "Electric field amplitude: 1e-06 a.u.\n", "[[ 1.0757e+01 0.0000e+00 -1.0000e-03]\n", " [ 0.0000e+00 1.0757e+01 -1.0000e-03]\n", " [-2.0000e-03 -2.0000e-03 1.5255e+01]]\n", "\n", "Electric field amplitude: 1e-05 a.u.\n", "[[ 1.07562e+01 0.00000e+00 -1.00000e-03]\n", " [ 0.00000e+00 1.07562e+01 -1.00000e-03]\n", " [-4.20000e-03 -4.20000e-03 1.52555e+01]]\n", "\n", "Electric field amplitude: 0.0001 a.u.\n", "[[ 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", "\n", "Electric field amplitude: 0.001 a.u.\n", "[[ 1.0756953e+01 4.1900000e-04 -7.9200000e-04]\n", " [ 4.1900000e-04 1.0756953e+01 -7.9200000e-04]\n", " [ 5.5100000e-04 5.5100000e-04 1.5260930e+01]]\n", "\n", "Electric field amplitude: 0.01 a.u.\n", "[[1.07220953e+01 2.30695000e-03 1.29460000e-03]\n", " [2.30695000e-03 1.07220953e+01 1.29460000e-03]\n", " [8.62080000e-03 8.62080000e-03 1.52781930e+01]]\n", "\n", "Electric field amplitude: 0.1 a.u.\n", "[[26.87530953 -0.10758047 -0.06919847]\n", " [-0.10758047 26.87530953 -0.06919847]\n", " [-0.0962627 -0.0962627 30.9737193 ]]\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is growing (delta= 3.34E-04) switch to SD\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Found an energy value lower than the FINAL energy (delta= 3.18E-04)\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: No convergence within the allowed number of minimization steps\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Wavefunctions not converged after cycle 1\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Self-consistent cycle did not meet convergence criteria\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is growing (delta= 1.62E-05) switch to SD\n", " warnings.warn(warning, UserWarning)\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is now decreasing again, coming back to DIIS\n", " warnings.warn(warning, UserWarning)\n" ] } ], "source": [ "ef_amplitudes = [1.e-11, 1.e-10, 1.e-9, 1.e-8, 1.e-7, 1.e-6, 1.e-5, 1.e-4, 1.e-3, 1.e-2, 1.e-1]\n", "# Loop over the electric field amplitudes to compute various\n", "# polarizability tensors\n", "pt_wfs = []\n", "for ef in ef_amplitudes:\n", " gs = Job(inputparams=inp, posinp=pos, name=\"N2\",\n", " run_dir=\"N2/pol_tensor/ef_amplitude_{}\".format(ef))\n", " pt_wf = PolTensor(gs, ef_amplitude=ef)\n", " pt_wf.run(nmpi=6, nomp=3)\n", " pt_wfs.append(pt_wf)\n", "\n", "for ef, pt_wf in zip(ef_amplitudes, pt_wfs):\n", " print(\"Electric field amplitude: {} a.u.\".format(ef))\n", " print(pt_wf.pol_tensor)\n", " print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The value of the amplitude of the electric field is subject to a certain balance: \n", "\n", "- if it is too high, the system undergoes a too large perturbation, leading to incorrect results (see how the two largest electric field amplitudes studied here give different results than the other ones; using a too large electric field even raises some warnings).\n", "- if it is too low, then the perturbation is too small and gives dipole differences that cannot be accurately measured (see how increasing the electric field amplitude gives more significant digits).\n", "\n", "That value must be chosen with great care, as it depends on the size of the system and on the grid extension used. The larger the system in a given direction (*e.g.*, a long, planar molecule such as an alkane), the larger the potential difference between both borders of the simulation box when an electric field is applied in that direction, meaning the larger the perturbation felt by the system. For that reason, it might even be better to apply an electric field with a smaller amplitude in that particular direction." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* ### How does the polarizability tensor depend on the system geometry?\n", "\n", "Modify the structure of the system by varying the distance between both atoms and comment on how this affects the value of polarizability tensor. Use the default electric field amplitude and the same base input parameters." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0735\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0735/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0735/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0735/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0835\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0835/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0835/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0835/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0935\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0935/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0935/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.0935/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1035\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1035/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1035/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1035/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1135\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1135/EF_along_x+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1135/EF_along_y+\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/distance_1.1135/EF_along_z+\n", "Logfile log-N2.yaml already exists!\n", "\n", "Distance: 1.0735 angstroem.\n", "[[ 1.076804e+01 1.105000e-02 9.520000e-03]\n", " [ 1.105000e-02 1.076804e+01 9.520000e-03]\n", " [-2.281000e-02 -2.281000e-02 1.523523e+01]]\n", "\n", "Distance: 1.0835 angstroem.\n", "[[ 1.076065e+01 3.660000e-03 2.130000e-03]\n", " [ 3.660000e-03 1.076065e+01 2.130000e-03]\n", " [-9.060000e-03 -9.060000e-03 1.524898e+01]]\n", "\n", "Distance: 1.0935 angstroem.\n", "[[ 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", "\n", "Distance: 1.1035 angstroem.\n", "[[ 1.075044e+01 -6.550000e-03 -8.080000e-03]\n", " [-6.550000e-03 1.075044e+01 -8.080000e-03]\n", " [ 6.920000e-03 6.920000e-03 1.526496e+01]]\n", "\n", "Distance: 1.1135 angstroem.\n", "[[ 1.074534e+01 -1.165000e-02 -1.318000e-02]\n", " [-1.165000e-02 1.074534e+01 -1.318000e-02]\n", " [ 1.428000e-02 1.428000e-02 1.527232e+01]]\n", "\n" ] } ], "source": [ "distances = [1.0735, 1.0835, 1.0935, 1.1035, 1.1135]\n", "# Loop over the interatomic distances to compute various\n", "# polarizability tensors\n", "pt_wfs = []\n", "for distance in distances:\n", " new_atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, distance])]\n", " new_pos = Posinp(new_atoms, units=\"angstroem\", boundary_conditions=\"free\")\n", " gs = Job(inputparams=inp, posinp=new_pos, name=\"N2\",\n", " run_dir=\"N2/pol_tensor/distance_{}\".format(distance))\n", " pt_wf = PolTensor(gs)\n", " pt_wf.run(nmpi=6, nomp=3)\n", " pt_wfs.append(pt_wf)\n", " \n", "for distance, pt_wf in zip(distances, pt_wfs):\n", " print(\"Distance: {} angstroem.\".format(distance))\n", " print(pt_wf.pol_tensor)\n", " print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The polarizability tensor slightly depends on the distance between the nitrogen atoms:\n", "\n", "- the $\\alpha_{xx}$ and $\\alpha_{yy}$ terms decrease while the $\\alpha_{zz}$ term increases when the distance is increased,\n", "- the non-diagonal elements are larger when the molecule is not in its optimal structure (that minimizes the forces, that is when the distance is 1.0935 $\\unicode[serif]{xC5}$),\n", "- the mean polarizability is only slightly altered (see the cell below)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Distance: 1.0735 angstroem\n", "12.257\n", "\n", "Distance: 1.0835 angstroem\n", "12.257\n", "\n", "Distance: 1.0935 angstroem\n", "12.257\n", "\n", "Distance: 1.1035 angstroem\n", "12.255\n", "\n", "Distance: 1.1135 angstroem\n", "12.254\n", "\n" ] } ], "source": [ "for distance, pt_wf in zip(distances, pt_wfs):\n", " print(\"Distance: {} angstroem\".format(distance))\n", " print(\"{:.3f}\".format(pt_wf.mean_polarizability))\n", " print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Improve that class!\n", "\n", "Another interest of defining a class is that you can easily improve it. You then only have to implement minor changes to your already performed scripts or notebooks to see those changes.\n", "\n", "For instance, you often find the mean polarizability in the literature instead of the tensor itself. It is defined as the mean value of the polarizability tensor diagonal elements. Given that you already compute that tensor in the post-processing procedure, the addition of another post-processing attribute named `mean_polarizability` seems like a good idea." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* ### Add the mean polarizability post-processing attribute\n", "\n", "The goal here is to compute the mean polarizability while post-processing the calculations and make it available via an attribute. You therefore have to perform three changes to your class:\n", "\n", "- add `\"mean_polarizability\"` to `POST_PROCESSING_ATTRIBUTES`\n", "- define the `mean_polarizability` attribute as a property (hint: use `pol_tensor` as a template).\n", "- compute the mean polarizabilty in the `post_proc` method (hint: the sum of the diagonal elements of a tensor can be easily computed with numpy: simply use `np.trace()`) and set the `mean_polarizability` attribute via its private counterpart (hint: again, use what is done for `pol_tensor` as a template).\n", "\n", "Once you are happy with your changes, run the whole notebook (click on \"Restart Kernel and run all cells\" under the \"Kernel\" tab). You should see that you did not need to modify any cell above to run the tests below. This is a reason why it is interesting to work with classes." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "assert np.isclose(pt_wf_z.mean_polarizability, 12.25712)\n", "assert pt_wf_x.mean_polarizability == pt_wf_y.mean_polarizability == pt_wf_z.mean_polarizability" ] } ], "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 }