Solution

Here is the solution to the second exercise regarding the polarizability tensor.

Define your worflow

As already mentioned, a workflow is mainly made of a queue of BigDFT jobs. This means you need to define a list of Job instances, to be run sequentially afterwards.

We want to compute the polarizability tensor of the N\(_2\) molecule using the Workflow class, just like in the previous exercise. The main input parameters being given, you only have to add Job instances to the pt_queue list. You can do so by re-using the same jobs as in the previous notebook (especially the same electric fields and run directories): first, a ground state job (without electric field) and then three calculations with an electric field, along the \(x\), \(y\) and \(z\) axis respectively. Try to use a for loop to add the last three jobs to the queue.

Only the base classes of the MyBigDFT package are required to do that, and might also require the help of extra packages such as numpy (you should not need it yet, though).

import os
from copy import deepcopy
from mybigdft import Workflow, InputParams, Posinp, Atom, Job
import numpy as np

# Main input variables:
# - use the optimized structure for the N2 molecule (found using the
#   LDA exchange-correlation functional)
atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0935])]
pos = Posinp(atoms, units="angstroem", boundary_conditions="free")
# - use some non-default input parameters to get high quality results
inp = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35}})
# Electric field amplitude to be applied in each direction
ef_amplitude = 1.e-4

# Initialize the queue of jobs in order to be able to compute the
# polarizability tensor:
# - Add the ground state job (without any electric field)
base_run_dir = os.path.relpath("N2/pol_tensor/")
gs = Job(inputparams=inp, posinp=pos, name="N2",
         run_dir=base_run_dir)
pt_queue = [gs]
# - Add the jobs with an electric field
for i, coord in enumerate(["x", "y", "z"]):
    # Add the electric field to the base input parameters
    elecfield = [0.0] * 3
    elecfield[i] = ef_amplitude
    new_inp = deepcopy(inp)
    new_inp["dft"]["elecfield"] = elecfield
    # Define the correct run directory
    new_run_dir = os.path.join(gs.run_dir,
                               "EF_along_{}+/".format(coord))
    # Append the new job to the end of the queue
    pt_queue.append(Job(inputparams=new_inp, posinp=pos,
                        name="N2", run_dir=new_run_dir))
pt_wf = Workflow(queue=pt_queue)

# Let us ensure that you followed the instructions:
# - there are four jobs in the queue
assert len(pt_wf.queue) == 4
# - correct electric fields are applied
assert "elecfield" not in pt_wf.queue[0].inputparams["dft"]
assert pt_wf.queue[1].inputparams["dft"]["elecfield"] == [1.e-4, 0.0, 0.0]
assert pt_wf.queue[2].inputparams["dft"]["elecfield"] == [0.0, 1.e-4, 0.0]
assert pt_wf.queue[3].inputparams["dft"]["elecfield"] == [0.0, 0.0, 1.e-4]
# - correct run directories
assert pt_wf.queue[0].run_dir.endswith("N2/pol_tensor")
assert pt_wf.queue[1].run_dir.endswith("N2/pol_tensor/EF_along_x+")
assert pt_wf.queue[2].run_dir.endswith("N2/pol_tensor/EF_along_y+")
assert pt_wf.queue[3].run_dir.endswith("N2/pol_tensor/EF_along_z+")

Run the jobs of your workflow

Nothing important to code here, as you might only want to change the value of the attributes to suit your needs.

Note that you do not have to use the context manager: no need for a with statement before using the run method of a Workflow instance as it is actually taken care of internally. This is one of the advantage of using a workflow.

pt_wf.run(nmpi=6, nomp=3)

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor
Logfile log-N2.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_x+
Logfile log-N2.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_y+
Logfile log-N2.yaml already exists!

/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/EF_along_z+
Logfile log-N2.yaml already exists!

Note that, if you succeeded in performing the previous notebook, the calculations were not run again.

# Let us ensure that you followed the instructions:
assert pt_wf.is_completed

Post-process your results

Now that you ran the workflow, you can use the results to extract meaningful quantities. Let us define a post_proc function to do just that. Its only argument must be a workflow instance. You should be able to do the post-processing without having to add another argument to the post_proc function. Try to loop over the three calculations with an electric field to get their dipoles.

def post_proc(workflow):
    """
    Use the logfile of the calculations run by the workflow to
    extract the polarizability tensor.

    Parameters
    ----------
    workflow: Workflow
        Workflow of run calculations allowing to compute the
        polarizability tensor of a given system.

    Returns
    -------
    numpy array
        Polarizability tensor of the system.
    """
    pol_tensor = np.zeros((3, 3))
    delta_ef = 1.e-4
    d0 = np.array(workflow.queue[0].logfile.dipole)
    for i, job in enumerate(workflow.queue[1:4]):
        d1 = np.array(job.logfile.dipole)
        pol_tensor[:, i] = (d1 - d0) / delta_ef
    return pol_tensor

# Run the post-processing of the workflow to compute
pol_tensor = post_proc(pt_wf)
print(pol_tensor)

[[ 1.075753e+01  5.400000e-04 -9.900000e-04]
 [ 5.400000e-04  1.075753e+01 -9.900000e-04]
 [-1.740000e-03 -1.740000e-03  1.525630e+01]]

# Let us ensure that you obtained the expected result:
expected_pol_tensor = [[ 1.075753e+01,  5.400000e-04, -9.900000e-04],
                       [ 5.400000e-04,  1.075753e+01, -9.900000e-04],
                       [-1.740000e-03, -1.740000e-03,  1.525630e+01]]
assert np.allclose(pol_tensor, expected_pol_tensor)