Solution

Here is the solution to the first exercise regarding 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 = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.0935])]
gs_pos = Posinp(atoms, units="angstroem", boundary_conditions="free")

# Let us continue with the input parameters
gs_inp = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35}})

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

import os
from mybigdft import Job

# Define the job and run it, using the above indications
base_run_dir = os.path.relpath("N2/pol_tensor")
with Job(posinp=gs_pos, inputparams=gs_inp, name="N2",
         run_dir=base_run_dir) as gs_job:
    gs_job.run(nmpi=6, nomp=3)

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

# 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 = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35,
                        "elecfield": [1.e-4, 0.0, 0.0]}})

# Define the job and run it, using the above indications
with Job(posinp=gs_pos, inputparams=EF_x_inp, name="N2",
         run_dir=os.path.join(base_run_dir, "EF_along_x+")) as EF_x_job:
    EF_x_job.run(nmpi=6, nomp=3)

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

# 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
d_0 = np.array(gs_job.logfile.dipole)
d_x = np.array(EF_x_job.logfile.dipole)
delta_ef = 1.e-4
pol_x_column = (d_x - d_0) / delta_ef
pol_tensor[:, 0] = pol_x_column
print(pol_tensor)

[[ 1.075753e+01  0.000000e+00  0.000000e+00]
 [ 5.400000e-04  0.000000e+00  0.000000e+00]
 [-1.740000e-03  0.000000e+00  0.000000e+00]]

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
EF_y_inp = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35,
                        "elecfield": [0.0, 1.e-4, 0.0]}})
with Job(posinp=gs_pos, inputparams=EF_y_inp, name="N2",
         run_dir=os.path.join(base_run_dir, "EF_along_y+")) as EF_y_job:
    EF_y_job.run(nmpi=6, nomp=3)

# Define and run a job with an electric field along the the z direction
EF_z_inp = InputParams({"dft": {"rmult": [7, 9], "hgrids": 0.35,
                        "elecfield": [0.0, 0.0, 1.e-4]}})
with Job(posinp=gs_pos, inputparams=EF_z_inp, name="N2",
         run_dir=os.path.join(base_run_dir, "EF_along_z+")) as EF_z_job:
    EF_z_job.run(nmpi=6, nomp=3)

/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!

# 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
for i, job in enumerate([EF_y_job, EF_z_job]):
    d_i = np.array(job.logfile.dipole)
    pol_tensor[:, i+1] = (d_i - d_0) / delta_ef
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)