Compute the polarizability tensor

The goal of this notebook is to make you use the basic API of the MyBigDFT module. You will have to define the input parameters, the initial geometry and run some jobs in order to be able to extract the polarizability tensor of the N\(_2\) molecule from the various logfiles.

To that end, you will have to interact with the notebook by coding within a few cells. Each coding cell is followed by a testing cell, in order to assert whether you followed the instructions you were given. You should not have to modify those testing cells: you only have to modify the coding cell until the corresponding testing cell runs successfully. Once this is achieved, you can move on to the next coding cell. Please do not to modify the provided variable names: they are used in the testing cells!

To help you along this notebook, you might want to check:

If you feel lost, you can also go and check the solution of this exercise.

The polarizability tensor

You will focus on computing the polarizability tensor for the N\(_2\) molecule. It represents the response of the charges of a system (its dipole \(D\)) to the application of an external electric field \(E\).

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 (\(x\), \(y\) and \(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:

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

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

One way of computing the polarizability tensor is to perform three calculations where the system is subject to a different external electric field: one along the \(x\), \(y\) or \(z\) axis. You can extract the dipole of these calculations from the logfiles of the BigDFT calculations. To compute the variation of the dipole and of the electric field, you must also perform a reference calculation. The most convenient one is the ground state calculation, where all there is no external electric field. This means you have to perform four calculations in total before being able of computing the polarizability tensor.

Run the ground state calculation

You will first have to define the initial geometry and the input parameters given the following indications.

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.

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.

[ ]:

from mybigdft import Posinp, Atom, InputParams

# Let us start with the initial geometry
atoms = []  # Define the atoms of the system as a list of Atoms instance
gs_pos =  # Define a Posinp instance using the above list of atoms

# Let us continue with the input parameters
gs_inp =  # Define an InputParams instance

[ ]:

# Let us ensure that you followed the instructions:
# - The system is made of two nitrogen atoms
assert len(gs_pos) == 2 and all([atom.type == "N" for atom in gs_pos])
# - correct distance between the atoms
assert gs_pos.distance(0, 1) == 1.0935 and gs_pos.units == "angstroem"
# - correct boundary conditions
assert gs_pos.boundary_conditions == "free"
# - correct rmult
assert gs_inp["dft"]["hgrids"] == 0.35
# - correct hgrids
assert gs_inp["dft"]["rmult"] == [7, 9]
# - no external electric field is applied
assert "elecfield" not in gs_inp["dft"]

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

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.

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

[ ]:

from mybigdft import Job

# Define the job and run it, using the above indications
with Job() as gs_job:
    gs_job.run()

[ ]:

# Let us ensure that you followed the instructions:
# - correct posinp
assert gs_job.posinp == gs_pos == gs_job.logfile.posinp
# - correct input parameters
assert gs_job.inputparams == gs_inp == gs_job.logfile.inputparams
# - correct job name
assert gs_job.name == "N2"
# - correct run directory
assert gs_job.run_dir.endswith("N2/pol_tensor")

Run a job with an electric field in the \(x\) direction

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+/".

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.

[ ]:

# Define the new input parameters
EF_x_inp =

# Define the job and run it, using the above indications
with Job() as EF_x_job:
    EF_x_job.run()

[ ]:

# Let us ensure that you followed the instructions:
# - same basic input parameters
assert EF_x_inp["dft"]["rmult"] == gs_inp["dft"]["rmult"]
assert EF_x_inp["dft"]["hgrids"] == gs_inp["dft"]["hgrids"]
# - correct electric field
assert EF_x_job.logfile.inputparams["dft"]["elecfield"] == [1.e-4, 0.0, 0.0]
# - correct job name
assert EF_x_job.name == "N2"
# - correct run directory
assert EF_x_job.run_dir.endswith("N2/pol_tensor/EF_along_x+")

Start building the polarizability tensor

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 module to build the polarizability tensor.

[ ]:

import numpy as np

pol_tensor = np.zeros((3, 3))
# Update the polarizability tensor

print(pol_tensor)

[ ]:

expected_pol_tensor = [[ 1.075753e+01, 0.0, 0.0],
                       [ 5.400000e-04, 0.0, 0.0],
                       [-1.740000e-03, 0.0, 0.0]]
assert np.allclose(pol_tensor, expected_pol_tensor)

Apply an electric field in the \(y\) and \(z\) direction

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.

[ ]:

# Define and run a job with an electric field along the the y direction


# Define and run a job with an electric field along the the z direction


[ ]:

# Let us ensure that you followed the instructions:
# - correct electric fields applied
assert EF_y_job.logfile.inputparams["dft"]["elecfield"] == [0.0, 1.e-4, 0.0]
assert EF_z_job.logfile.inputparams["dft"]["elecfield"] == [0.0, 0.0, 1.e-4]
# - correct run directories
assert EF_y_job.run_dir.endswith("N2/pol_tensor/EF_along_y+")
assert EF_z_job.run_dir.endswith("N2/pol_tensor/EF_along_z+")

Complete the computation of the polarizability tensor

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.

[ ]:

# Define the rest of the polarizability tensor from the logfiles

print(pol_tensor)

[ ]:

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

If all the assertions ran successfully: congratulations! You made it to the end of this first exercise! The values can be compared to the one reported in table II and III of ` <https://aip.scitation.org/doi/10.1063/1.454327>`__. The value of \(\alpha_{zz}\) is reported to be between 14.8 and 15.0 a.u. (for atomic units), while the value of \(\alpha_{xx} = \alpha_{yy}\) is reported to be between 9.0 and 10.1 a.u. Our values respectively are 15.3 and 10.8 a.u.

You can now go on to the second exercise, where you will start to automate the workflow you successfully coded here.