Solution¶

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

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!

[1]:

import os
from copy import deepcopy
import numpy as np
from mybigdft import Job
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.
    """

    POST_PROCESSING_ATTRIBUTES = ["pol_tensor"]

    def __init__(self, ground_state, ef_amplitude=1.e-4):
        """
        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 = [self.ground_state]
        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(ground_state.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))
        # 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

    @property
    def mean_polarizability(self):
        """
        Set mean_polarizability 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._mean_polarizability

    def post_proc(self):
        """
        Compute and set the polarizability tensor and the mean polarizability
        of the system.
        """
        # Polarizability tensor
        pol_tensor = np.zeros((3, 3))
        delta_ef = max(self.queue[1].inputparams["dft"]["elecfield"])
        d0 = np.array(self.queue[0].logfile.dipole)
        for i, job in enumerate(self.queue[1:4]):
            d1 = np.array(job.logfile.dipole)
            pol_tensor[:, i] = (d1 - d0) / delta_ef
        self._pol_tensor = pol_tensor
        # Mean polarizability
        self._mean_polarizability = np.trace(self.pol_tensor) / 3.

[2]:

# 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 """AttributeError("can't set attribute""" in repr(e)
    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 """AttributeError("can't set attribute""" in repr(e)
    print("Correct error catched!")

Correct error catched!
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.

[3]:

# Compute the polarizability tensor when the molecule is along the x axis
atoms_x = [Atom('N', [0, 0, 0]), Atom('N', [1.0935, 0, 0])]
pos_x = Posinp(atoms_x, units="angstroem", boundary_conditions="free")
run_dir_x = "N2/pol_tensor/atoms_along_x/"
gs_x = Job(inputparams=inp, posinp=pos_x, name="N2", run_dir=run_dir_x)
pt_wf_x = PolTensor(gs_x)
pt_wf_x.run(nmpi=6, nomp=3)
print(pt_wf_x.pol_tensor)

# Compute the polarizability tensor when the molecule is along the x axis
atoms_y = [Atom('N', [0, 0, 0]), Atom('N', [0, 1.0935, 0])]
pos_y = Posinp(atoms_y, units="angstroem", boundary_conditions="free")
run_dir_y = "N2/pol_tensor/atoms_along_y/"
gs_y = Job(inputparams=inp, posinp=pos_y, name="N2", run_dir=run_dir_y)
pt_wf_y = PolTensor(gs_x)
pt_wf_y.run(nmpi=6, nomp=3)
print(pt_wf_y.pol_tensor)

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

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

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

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

[[ 1.116730e+01  4.103100e-01  4.087800e-01]
 [ 5.400000e-04  1.075753e+01 -9.900000e-04]
 [-4.115100e-01 -4.115100e-01  1.484653e+01]]
/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/pol_tensor/atoms_along_x
Logfile log-N2.yaml already exists!

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

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

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

[[ 1.116730e+01  4.103100e-01  4.087800e-01]
 [ 5.400000e-04  1.075753e+01 -9.900000e-04]
 [-4.115100e-01 -4.115100e-01  1.484653e+01]]

[4]:

# 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(nmpi=6, nomp=3)
print(pt_wf_z.pol_tensor)
assert np.isclose(np.trace(pt_wf_x.pol_tensor), np.trace(pt_wf_z.pol_tensor))

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

[[ 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]]

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.

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. Keep the same structure and input parameters as above.

[5]:

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
pt_wfs = []
for ef in ef_amplitudes:
    gs = Job(inputparams=inp, posinp=pos, name="N2",
             run_dir="N2/pol_tensor/ef_amplitude_{}".format(ef))
    pt_wf = PolTensor(gs, ef_amplitude=ef)
    pt_wf.run(nmpi=6, nomp=3)
    pt_wfs.append(pt_wf)

for ef, pt_wf in zip(ef_amplitudes, pt_wfs):
    print("Electric field amplitude: {} a.u.".format(ef))
    print(pt_wf.pol_tensor)
    print()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Electric field amplitude: 1e-11 a.u.
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]

Electric field amplitude: 1e-10 a.u.
[[10.  0.  0.]
 [ 0. 10.  0.]
 [ 0.  0. 10.]]

Electric field amplitude: 1e-09 a.u.
[[11.  0.  0.]
 [ 0. 11.  0.]
 [ 0.  0. 15.]]

Electric field amplitude: 1e-08 a.u.
[[10.8  0.   0. ]
 [ 0.  10.8  0. ]
 [ 0.   0.  15.2]]

Electric field amplitude: 1e-07 a.u.
[[10.76  0.    0.  ]
 [ 0.   10.76  0.  ]
 [ 0.    0.   15.25]]

Electric field amplitude: 1e-06 a.u.
[[ 1.0757e+01  0.0000e+00 -1.0000e-03]
 [ 0.0000e+00  1.0757e+01 -1.0000e-03]
 [-2.0000e-03 -2.0000e-03  1.5255e+01]]

Electric field amplitude: 1e-05 a.u.
[[ 1.07562e+01  0.00000e+00 -1.00000e-03]
 [ 0.00000e+00  1.07562e+01 -1.00000e-03]
 [-4.20000e-03 -4.20000e-03  1.52555e+01]]

Electric field amplitude: 0.0001 a.u.
[[ 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]]

Electric field amplitude: 0.001 a.u.
[[ 1.0756953e+01  4.1900000e-04 -7.9200000e-04]
 [ 4.1900000e-04  1.0756953e+01 -7.9200000e-04]
 [ 5.5100000e-04  5.5100000e-04  1.5260930e+01]]

Electric field amplitude: 0.01 a.u.
[[1.07220953e+01 2.30695000e-03 1.29460000e-03]
 [2.30695000e-03 1.07220953e+01 1.29460000e-03]
 [8.62080000e-03 8.62080000e-03 1.52781930e+01]]

Electric field amplitude: 0.1 a.u.
[[26.87530953 -0.10758047 -0.06919847]
 [-0.10758047 26.87530953 -0.06919847]
 [-0.0962627  -0.0962627  30.9737193 ]]

/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is growing (delta= 3.34E-04) switch to SD
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Found an energy value lower than the FINAL energy (delta= 3.18E-04)
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: No convergence within the allowed number of minimization steps
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Wavefunctions not converged after cycle  1
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: Self-consistent cycle did not meet convergence criteria
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is growing (delta= 1.62E-05) switch to SD
  warnings.warn(warning, UserWarning)
/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/iofiles/logfiles.py:270: UserWarning: The energy value is now decreasing again, coming back to DIIS
  warnings.warn(warning, UserWarning)

The value of the amplitude of the electric field is subject to a certain balance:

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

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.

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.

[6]:

distances = [1.0735, 1.0835, 1.0935, 1.1035, 1.1135]
# Loop over the interatomic distances to compute various
# polarizability tensors
pt_wfs = []
for distance in distances:
    new_atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, distance])]
    new_pos = Posinp(new_atoms, units="angstroem", boundary_conditions="free")
    gs = Job(inputparams=inp, posinp=new_pos, name="N2",
             run_dir="N2/pol_tensor/distance_{}".format(distance))
    pt_wf = PolTensor(gs)
    pt_wf.run(nmpi=6, nomp=3)
    pt_wfs.append(pt_wf)

for distance, pt_wf in zip(distances, pt_wfs):
    print("Distance: {} angstroem.".format(distance))
    print(pt_wf.pol_tensor)
    print()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Distance: 1.0735 angstroem.
[[ 1.076804e+01  1.105000e-02  9.520000e-03]
 [ 1.105000e-02  1.076804e+01  9.520000e-03]
 [-2.281000e-02 -2.281000e-02  1.523523e+01]]

Distance: 1.0835 angstroem.
[[ 1.076065e+01  3.660000e-03  2.130000e-03]
 [ 3.660000e-03  1.076065e+01  2.130000e-03]
 [-9.060000e-03 -9.060000e-03  1.524898e+01]]

Distance: 1.0935 angstroem.
[[ 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]]

Distance: 1.1035 angstroem.
[[ 1.075044e+01 -6.550000e-03 -8.080000e-03]
 [-6.550000e-03  1.075044e+01 -8.080000e-03]
 [ 6.920000e-03  6.920000e-03  1.526496e+01]]

Distance: 1.1135 angstroem.
[[ 1.074534e+01 -1.165000e-02 -1.318000e-02]
 [-1.165000e-02  1.074534e+01 -1.318000e-02]
 [ 1.428000e-02  1.428000e-02  1.527232e+01]]

The polarizability tensor slightly depends on the distance between the nitrogen atoms:

the \(\alpha_{xx}\) and \(\alpha_{yy}\) terms decrease while the \(\alpha_{zz}\) term increases when the distance is increased,
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}\)),
the mean polarizability is only slightly altered (see the cell below).

[7]:

for distance, pt_wf in zip(distances, pt_wfs):
    print("Distance: {} angstroem".format(distance))
    print("{:.3f}".format(pt_wf.mean_polarizability))
    print()

Distance: 1.0735 angstroem
12.257

Distance: 1.0835 angstroem
12.257

Distance: 1.0935 angstroem
12.257

Distance: 1.1035 angstroem
12.255

Distance: 1.1135 angstroem
12.254

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.

[8]:

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