Code your own polarizability tensor workflow, part 2

The main goal of this notebook is to show how the code implemented in the previous exercise can be easily converted to a fully functional class deriving from the AbstractWorkflow class.

We will stick to the polarizability tensor problem, in order to build on what was achieved in the previous exercises. The structure of this notebook is very similar to the previous ones. You can use the same resources to help you (especially the), while the solution to this exercise can be found here.

The AbstractWorkflow class

You already saw that using the Workflow class allows to get a clean separation of the procedure to compute a quantity of interest:

  • first, the definition of the jobs to be run,

  • then the jobs are run,

  • and finally, the desired output can be computed from the results of the calculations.

The goal of the AbstractWorkflow class is to ease the implementation of this type of workflow in a class that is meant to be general, i.e., that can be applicable to a large class of systems, using a large set of possible input parameters. Such a class is meant to be used simply by:

  • initializing an instance of the class, setting the queue of jobs automatically,

  • running all the jobs of the workflow instance and performing the post-processing, all these actions being done by applying the run method.

And that’s it!

Note that the post-processing procedure, which is defined in the class body, is automatically performed after all the jobs in the workflow queue. The quantities of interest (here, the polarizability tensor) are even used to set semi-private attributes, so that they can be easily accessed afterwards (don’t worry, the last sentence should get clearer by the end of the notebook).

Another interest of writing a class is that it enables you to add more functionalities to your class afterwards, such as adding another quantity to compute while running the post-processing procedure. Such changes would be readily available in your scripts or notebooks where you used that class: no need to modify those files one by one to make that functionality available there.

To write such a class you will need to:

  • define which are the quantities to be computed while post-processing the results by:

    • setting the value of the POST_PROCESSING_ATTRIBUTES list, containing the names of all the quantities of interest computed while post-processing,

    • defining a property getter for each of these attributes, so as to access them afterwards,

  • override the __init__ method, in order to make sure than the queue of jobs is correctly initialized,

  • define the post_proc method to compute and then set the values of the quantities of interest.

Create a class

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.

Part of the body of the class is already written, you must only code the following:

  • add the string "pol_tensor" to the POST_PROCESSING_ATTRIBUTES list,

  • initialize the queue of jobs in the __init__ method,

  • define the post_proc method to compute the value of the pol_tensor attribute.

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).

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!

[ ]:

import numpy as np
# Do not forget to import some relevant stuff...
from mybigdft.workflows.workflow import AbstractWorkflow

class PolTensor(AbstractWorkflow):
    """
    This workflow allows to compute the (electronic) polarizability
    tensor of a given system.

    The polarizability tensor represents the response of the charges of
    a system (its dipole) to the application of an external electric
    field.

    To compute this polarizability tensor, some BigDFT calculations are
    performed, where the system is subject to an external electric
    field along each direction of space (:math:`x`, :math:`y` and
    :math:`z`). The elements of the polarizability tensor are then
    defined by the ratio of the delta of the dipole in one direction and
    the delta of the electric field amplitudes:

    .. math::

        \alpha_{ij} = \frac{\Delta D_i}{\Delta E_j}

    where :math:`i, j \in \{x, y, z\}`, :math:`\Delta D_i` is the
    variation of the dipole along the :math:`i` direction and
    :math:`\Delta E_j` is the variation of the electric field amplitude
    along the :math:`j` direction.
    """

    # Add "pol_tensor" to the list
    POST_PROCESSING_ATTRIBUTES = []

    # Add ef_amplitude optional argument (default value: 1.e-4)
    def __init__(self, ground_state, ):
        """
        A PolTensor workflow is initialized by the job of the
        ground-state of the system and an electric field amplitude.

        Parameters
        ----------
        ground_state: Job
            Reference job, specifying the base input parameters (where no
            electric field is applied), the geometry of the system and the
            base run directory of the calculations.
        ef_amplitude: float
            Electric field amplitude to be applied along each direction,
            in atomic units (default to 1.e-4).
        """
        self._ground_state = ground_state
        # Set the queue of jobs thanks to ground_state and ef_amplitude
        pt_queue = []
        # This calls the __init__ method of the AbstractWorkflow class to
        # set the queue, the run method (that calls the post_proc method
        # automatically) and the default value of the post-processing
        # attributes to None.
        super().__init__(queue=pt_queue)

    @property
    def ground_state(self):
        """
        Set ground_state as a semi-private attribute: it can be accessed,
        but it cannot be set.

        Returns
        -------
        Job
            Reference job, where no electric field is applied.
        """
        return self._ground_state

    @property
    def pol_tensor(self):
        """
        Set pol_tensor as a semi-private attribute: it can be accessed,
        but it cannot be set.

        Returns
        -------
        numpy array of shape (3, 3)
            Polarizability tensor of the system, using :attr:`ground_state`
            as a reference job.
        """
        return self._pol_tensor

    def post_proc(self):
        """
        Compute and set the polarizability tensor of the system.
        """
        pol_tensor = np.zeros((3, 3))
        # Compute the polarizability tensor thanks to the jobs in the queue
        self._pol_tensor = pol_tensor

[ ]:

# Let us test that class with for the N2 molecule
from mybigdft import Job, Posinp, Atom, InputParams
atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0935])]
pos = Posinp(atoms, units="angstroem", boundary_conditions="free")
inp = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35}})

# 1- Using non-default electric field amplitude:
gs_1 = Job(posinp=pos, inputparams=inp, name="N2",
           run_dir="N2/pol_tensor/non_default/")
pt_wf_1 = PolTensor(gs_1, ef_amplitude=5.e-4)
#   - correct electric field
expected_efs_1 = [None, [5.e-4, 0, 0], [0, 5.e-4, 0], [0, 0, 5.e-4]]
for job, ef in zip(pt_wf_1.queue, expected_efs_1):
    assert job.inputparams["dft"].get("elecfield") == ef
#   - correct run directories
expected_run_dirs_1 = ["non_default/", "non_default/EF_along_x+",
                       "non_default/EF_along_y+", "non_default/EF_along_z+"]
for job, run_dir in zip(pt_wf_1.queue, expected_run_dirs_1):
    assert job.run_dir.endswith(run_dir)
#   - correct pol_tensor attribute
assert pt_wf_1.pol_tensor is None
try:
    pt_wf_1.pol_tensor = np.zeros((3, 3))
except AttributeError as e:
    assert repr(e) == """AttributeError("can't set attribute",)"""
    print("Correct error catched!")

# 2- Using default electric field amplitude:
gs_2 = Job(posinp=pos, inputparams=inp, name="N2",
           run_dir="N2/pol_tensor/default/")
pt_wf_2 = PolTensor(gs_2)
#   - correct electric field
expected_efs_2 = [None, [1.e-4, 0, 0], [0, 1.e-4, 0], [0, 0, 1.e-4]]
for job, ef in zip(pt_wf_2.queue, expected_efs_2):
    assert job.inputparams["dft"].get("elecfield") == ef
#   - correct run directories
expected_run_dirs_2 = ["default/", "default/EF_along_x+",
                       "default/EF_along_y+", "default/EF_along_z+"]
for job, run_dir in zip(pt_wf_2.queue, expected_run_dirs_2):
    assert job.run_dir.endswith(run_dir)
#   - correct pol_tensor attribute
assert pt_wf_2.pol_tensor is None
try:
    pt_wf_2.pol_tensor = np.zeros((3, 3))
except AttributeError as e:
    assert repr(e) == """AttributeError("can't set attribute",)"""
    print("Correct error catched!")

Use that class!

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.

  • The polarizability tensor must not depend on the rotation of the initial molecule

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.

[ ]:

# Compute the polarizability tensor when the molecule is along the x axis
atoms_x = []
pos_x = Posinp(atoms_x, )
run_dir_x = "N2/pol_tensor/atoms_along_x/"
pt_wf_x =
pt_wf_x.run()

# Compute the polarizability tensor when the molecule is along the x axis
atoms_y = []
pos_y = Posinp(atoms_y, )
run_dir_y = "N2/pol_tensor/atoms_along_y/"
pt_wf_y =
pt_wf_y.run()

[ ]:

# Let us ensure that you obtained the expected results:
# - both tensors are the same
assert np.allclose(pt_wf_x.pol_tensor, pt_wf_y.pol_tensor)
# - they have the same trace as the polarizability tensor where the molecule
#   is along the z axis
gs_z = Job(posinp=pos, inputparams=inp, name="N2", run_dir="N2/pol_tensor/")
pt_wf_z = PolTensor(gs_z)
pt_wf_z.run()
print(pt_wf_z.pol_tensor)
assert np.isclose(np.trace(pt_wf_x.pol_tensor), np.trace(pt_wf_z.pol_tensor))

  • What is the influence of the electric field amplitude on the results?

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.

[ ]:

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]
# Loop over the electric field amplitudes to compute various
# polarizability tensors

  • How does the polarizability tensor depend on the system geometry?

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.

[ ]:

distances = [1.0735, 1.0835, 1.0935, 1.1035, 1.1135]
# Loop over the interatomic distances to compute various
# polarizability tensors

Improve that class!

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.

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.

  • Add the mean polarizability post-processing attribute

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:

  • add "mean_polarizability" to POST_PROCESSING_ATTRIBUTES

  • define the mean_polarizability attribute as a property (hint: use pol_tensor as a template).

  • 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).

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.

[ ]:

assert np.isclose(pt_wf_z.mean_polarizability, 12.25712)
assert pt_wf_x.mean_polarizability == pt_wf_y.mean_polarizability == pt_wf_z.mean_polarizability

  • Further step (optional)

You could add the possibility to perform two BigDFT calculations per space coordinates in order to get more accurate results, while removing the reference job without any electric field. This means that six jobs must be run instead of four: two per space coordinate - one with a positive electric field amplitude, the other with a negative one (the run directory in that case must end by x- if the electric field is along the \(x\) direction). This is especially relevant when the system under consideration has no dipole: BigDFT would still find a residual dipole, that might be large enough to give messy results (depending on the required accuracy).

This exercise is optional as it might be a bit tedious for beginners. The main difficulty is to keep the ability of running only four jobs to get the polarizability tensor and the mean polarizability: you will have to add an argument order to the __init__ method and make sure it triggers to correct behavior:

  • if order=1, then four jobs should be performed,

  • if order=2, then six jobs should be performed.

This means you must make sure that the initialization of the queue depends on that parameter, ans so must be the computation of the polarizability tensor at the post-processing level.

For those interested in doing it, you may want to use the same scheme as above: first use the Workflow class to set up the workflow using six different calculations before actually implementing that into the PolTensor class. You might even want to translate the workflow with six jobs to a new class deriving from AbstractWorkflowbefore actually modifying the PolTensor class. This should allow you to see which part of the code you have to modify, and to do it in a minimal fashion.