Back to the main Index
This tutorial shows how to calculate the following physical properties related to strain:
The discussion is based on the tutorial on elastic properties available on the Abinit web site. More specifically, we will discuss how to
DdbRobot
to analyze the convergence.You might find additional material related to the present section in the following references (already mentioned in the official tutorial):
The first paper provides a detailed discussion of the theory underlying the incorporation of atom-relaxation corrections. We strongly recommend to read this article as this notebook will mainly focus on the usage of the python API assuming you are already familiar with the theoretical aspects. The second paper discusses in more details the DFPT treatment of strain perturbations in Abinit.
If you are already familiar with python and AbiPy-Abinit are already installed and configured, you may want to use directly the command line interface. There is a README.md file in the directory of this lesson explaining how to analyze the data from the shell using ipython and matplotlib.
Note: The code in this notebook requires abinit >= 8.9 and abipy >= 0.6
Before starting, we need to import the python modules and the functions needed in the notebook:
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import print_function, division, unicode_literals
import numpy as np
import warnings
warnings.filterwarnings("ignore") # to get rid of deprecation warnings
from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook
import abipy.flowtk as flowtk
# This line configures matplotlib to show figures embedded in the notebook.
# Replace `inline` with `notebook` in classic notebook
%matplotlib inline
# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib
#%matplotlib widget
and an useful function from the lesson_elastic
module required to generate our DFPT flows:
from lesson_elastic import make_scf_input
abilab.print_source(make_scf_input)
def make_scf_input(ngkpt=(4, 4, 4)):
"""
This function constructs the input file for the GS calculation of
AlAs in hypothetical wurzite (hexagonal) structure.
In principle, the stucture should be relaxed before starting the calculation,
here we use the *unrelaxed* geometry of the official tutorial.
Args:
ngkpt: K-mesh used both in the GS and in the DFPT part.
"""
# Initialize structure. Use enough significant digits
# so that Abinit will recognize the correct spacegroup
# (Hexagonal and rhombohedral lattices are a bit problematic).
structure = abilab.Structure.from_abivars(
acell=[7.5389648144E+00, 7.5389648144E+00, 1.2277795374E+01],
natom=4,
ntypat=2,
rprim=[ np.sqrt(0.75), 0.5, 0.0 ,
-np.sqrt(0.75), 0.5, 0.0,
0.0, 0.0, 1.0],
typat=[1, 1, 2, 2],
xred=[1/3, 2/3, 0,
2/3, 1/3, 1/2,
1/3, 2/3, 3.7608588373E-01,
2/3, 1/3, 8.7608588373E-01],
znucl=[13, 33],
)
pseudos = abidata.pseudos("13al.pspnc", "33as.pspnc")
gs_inp = abilab.AbinitInput(structure, pseudos=pseudos)
# Set other important variables (consistent with tutorial)
# All the other DFPT runs will inherit these parameters.
gs_inp.set_vars(
nband=8,
ecut=6.0,
ecutsm=0.5, # Important when performing structural optimization
# with variable cell. All DFPT calculations should use
# the same value to be consistent.
ngkpt=ngkpt,
nshiftk=1,
shiftk=[0.0, 0.0, 0.5], # This choice preserves the hexagonal symmetry of the grid.
diemac=9.0,
nstep=40,
paral_kgb=0,
tolvrs=1.0e-18,
)
return gs_inp
The function makes some assumptions for important parameters such as
the crystalline structure and the pseudos.
This is done on purpose to keep the code as simple as possible.
It should not be so difficult to generalize the implementation to take into account other cases.
Note how the function accepts an optional argument ngkpt
defining the k-mesh so that we can easily
change the sampling e.g. for convergence studies.
Let's start to play with our new function:
scf_input = make_scf_input()
scf_input
print(scf_input.structure)
Full Formula (Al2 As2) Reduced Formula: AlAs abc : 3.989448 3.989448 6.497130 angles: 90.000000 90.000000 120.000000 Sites (4) # SP a b c --- ---- -------- -------- -------- 0 Al 0.333333 0.666667 0 1 Al 0.666667 0.333333 0.5 2 As 0.333333 0.666667 0.376086 3 As 0.666667 0.333333 0.876086
scf_input.structure.plot();
We are using the same norm-conserving pseudopotentials of the official tutorial but this does not mean you should use them for production calculations (there must be a reason why the directory is called Psps_for_tests).
There are 16 valence electrons per unit cell hence nband
has been set to 8
(yes, we are dealing with a non-magnetic semiconductor):
for pseudo in scf_input.pseudos:
print(pseudo, "\n")
<NcAbinitPseudo: 13al.pspnc> summary: Troullier-Martins psp for element Al Thu Oct 27 17:31:05 EDT 1994 number of valence electrons: 3.0 maximum angular momentum: d angular momentum for local part: d XC correlation: LDA_XC_TETER93 supports spin-orbit: False radius for non-linear core correction: 2.09673076353074 hint for low accuracy: ecut: 0.0, pawecutdg: 0.0 hint for normal accuracy: ecut: 0.0, pawecutdg: 0.0 hint for high accuracy: ecut: 0.0, pawecutdg: 0.0 <NcAbinitPseudo: 33as.pspnc> summary: Troullier-Martins psp for element As Thu Oct 27 17:37:14 EDT 1994 number of valence electrons: 5.0 maximum angular momentum: p angular momentum for local part: p XC correlation: LDA_XC_TETER93 supports spin-orbit: False radius for non-linear core correction: 2.0573171556401 hint for low accuracy: ecut: 0.0, pawecutdg: 0.0 hint for normal accuracy: ecut: 0.0, pawecutdg: 0.0 hint for high accuracy: ecut: 0.0, pawecutdg: 0.0
The scf_input
represents the building block for our DFPT calculation.
As usual, AbiPy provides a Work
to compute elastic and piezoelectric properties
starting from an input representing a ground-state calculation
thus it is just a matter of calling make_scf_input
and pass the result to
ElasticWork.from_scf_input
to construct our flow.
Let's have a look at the actual implementation:
from lesson_elastic import build_flow
abilab.print_source(build_flow)
def build_flow(options=None):
"""
Create a `Flow` for phonon calculations. The flow has one work with:
- 1 GS Task
- 3 DDK Task
- 4 Phonon Tasks (Gamma point)
- 6 Elastic tasks (3 uniaxial + 3 shear strain)
The Phonon tasks and the elastic task will read the 3 DDK files produced at the beginning
"""
workdir = options.workdir if (options and options.workdir) else "flow_elastic"
flow = flowtk.Flow(workdir=workdir)
# Build input for GS calculation and register the first work.
scf_input = make_scf_input()
# Build work for elastic properties (clamped-ions)
# activate internal strain and piezoelectric part.
elast_work = flowtk.ElasticWork.from_scf_input(scf_input, with_relaxed_ion=True, with_piezo=True)
flow.register_work(elast_work)
return flow
Now we can call the function to build our flow:
flow = build_flow()
flow.show_info()
<Flow, node_id=351340, workdir=flow_elastic> Number of works: 1, total number of tasks: 14 Number of tasks with a given class: Task Class Number ------------ -------- ScfTask 1 DdkTask 3 PhononTask 4 ElasticTask 6
and use the get_graphviz
method to visualize the connection among the Tasks
:
flow.get_graphviz()
# matplotlib version based on networkx
#flow.plot_networkx(with_edge_labels=True);
In a nutshell:
WFK
file in the ScfTask
(red circle)DdkTasks
to compute $\dfrac{\partial u}{\partial{\bf k}}$ for the three different directions.ElasticTasks
needs the WFK
file to compute the six strain perturbations (3 for uniaxial and 3 for shear strain) while the DDK
files are required to compute the mixed 2nd-order derivatives with respect to strain and electric field needed for the piezoelectric tensor.Note that, contrarily to the approach used in the standard tutorial,
the AbiPy Work
does not use datasets.
The perturbations of interest (strain, atomic-perturbation, ddk, electric field)
are obtained with different Tasks
that can be executed in parallel.
To understand better this point, we can print a table with the most important Abinit variables defining
the DFPT calculation.
If the variable is not present in the input, the entry is set to None
.
flow.get_vars_dataframe("qpt", "rfphon", "rfatpol", "rfdir", "rfelfd", "rfstrs", "kptopt")
qpt | rfphon | rfatpol | rfdir | rfelfd | rfstrs | kptopt | class | |
---|---|---|---|---|---|---|---|---|
w0_t0 | None | None | None | None | None | None | None | ScfTask |
w0_t1 | (0, 0, 0) | None | None | (1, 0, 0) | 2 | None | 2 | DdkTask |
w0_t2 | (0, 0, 0) | None | None | (0, 1, 0) | 2 | None | 2 | DdkTask |
w0_t3 | (0, 0, 0) | None | None | (0, 0, 1) | 2 | None | 2 | DdkTask |
w0_t4 | (0, 0, 0) | 1 | [1, 1] | [1, 0, 0] | None | None | 2 | PhononTask |
w0_t5 | (0, 0, 0) | 1 | [1, 1] | [0, 0, 1] | None | None | 2 | PhononTask |
w0_t6 | (0, 0, 0) | 1 | [3, 3] | [1, 0, 0] | None | None | 2 | PhononTask |
w0_t7 | (0, 0, 0) | 1 | [3, 3] | [0, 0, 1] | None | None | 2 | PhononTask |
w0_t8 | (0, 0, 0) | None | None | [1, 0, 0] | None | 1 | 2 | ElasticTask |
w0_t9 | (0, 0, 0) | None | None | [0, 1, 0] | None | 1 | 2 | ElasticTask |
w0_t10 | (0, 0, 0) | None | None | [0, 0, 1] | None | 1 | 2 | ElasticTask |
w0_t11 | (0, 0, 0) | None | None | [1, 0, 0] | None | 2 | 2 | ElasticTask |
w0_t12 | (0, 0, 0) | None | None | [0, 1, 0] | None | 2 | 2 | ElasticTask |
w0_t13 | (0, 0, 0) | None | None | [0, 0, 1] | None | 2 | 2 | ElasticTask |
If the meaning of these variables is not clear, you can consult the Abinit documentation or access the documentation directly from python with e.g.:
abilab.docvar("rfstrs")
0
Used to run strain response-function calculations (e.g. needed to get elastic constants). Define, with rfdir, the set of perturbations.
See the possible restrictions on the use of strain perturbations, in the help:respfn.
Now we can generate the flow_elastic
directory with the input files by executing
the lesson_elastic.py
script.
Then use the abirun.py
script to launch the entire calculation with:
abirun.py flow_elastic scheduler
You will see that all PhononTasks
and ElasticTasks
are executed in parallel on your machine
once the three DdkTasks
are completed.
If you prefer to skip this part, you may want to jump to next section about the post-processing of the results. Note that the output files are already available in the repository so it is also possible to try the AbiPy post-processing tools without having to run the flow.
Our flow is completed and we have the final DDB file
in the outdata
directory of the Work
(AbiPy has automatically merged all the partial DDB files
at the end of the calculation by invoking anaddb
for you).
Let's open this DDB file with:
ddb = abilab.abiopen("flow_elastic/w0/outdata/out_DDB")
print(ddb)
================================= File Info ================================= Name: out_DDB Directory: /Users/gmatteo/git_repos/abitutorials/abitutorials/elastic/flow_elastic/w0/outdata Size: 27.03 kb Access Time: Sat Aug 18 16:03:31 2018 Modification Time: Mon Aug 13 17:47:44 2018 Change Time: Mon Aug 13 17:47:44 2018 ================================= Structure ================================= Full Formula (Al2 As2) Reduced Formula: AlAs abc : 3.989448 3.989448 6.497130 angles: 90.000000 90.000000 120.000000 Sites (4) # SP a b c cartesian_forces --- ---- -------- -------- -------- ------------------------------------------------- 0 Al 0.333333 0.666667 0 [-0.00000000e+00 -0.00000000e+00 4.06423231e-06] 1 Al 0.666667 0.333333 0.5 [-0.00000000e+00 -0.00000000e+00 4.06423231e-06] 2 As 0.333333 0.666667 0.376086 [-0.00000000e+00 -0.00000000e+00 -4.06423231e-06] 3 As 0.666667 0.333333 0.876086 [-0.00000000e+00 -0.00000000e+00 -4.06423231e-06] Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 12, has_timerev: True, symmorphic: True ================================== DDB Info ================================== Number of q-points in DDB: 1 guessed_ngqpt: [1 1 1] (guess for the q-mesh divisions made by AbiPy) Has total energy: True, Has forces: True Total energy: -551.6004759981904 eV [eV] Cartesian stress tensor in GPa with pressure -1.518e-07 (GPa): [[-1.11784777e-05 0.00000000e+00 0.00000000e+00] [ 0.00000000e+00 -1.11786387e-05 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 2.28125046e-05]] Has (at least one) atomic pertubation: True Has (at least one diagonal) electric-field perturbation: False Has (at least one) Born effective charge: True Has (all) strain terms: True Has (all) internal strain terms: True Has (all) piezoelectric terms: True
The DdbFile
object provides an easy-to-use interface that invokes anaddb
to post-process
the data stored in the DDB file.
All the methods that invoke anaddb
use the ana
prefix so it is not strange to
see that we can obtain the elastic and piezoelectric tensors by just calling:
edata = ddb.anaget_elastic(verbose=1)
ANADDB INPUT: asr 2 chneut 1 dieflag 0 elaflag 3 piezoflag 3 instrflag 1 workdir: /var/folders/89/47k8wfdj11x035svqf8qnl4m0000gn/T/tmpf3ep5ya9
Note that we are calling the method without arguments. This means that AbiPy will try to detect automatically how to set the anaddb input variables. The docstring of the method explains the logic used to set the variables.
abilab.print_doc(ddb.anaget_elastic)
def anaget_elastic(self, relaxed_ion="automatic", piezo="automatic",
dde=False, stress_correction=False, asr=2, chneut=1,
mpi_procs=1, workdir=None, manager=None, verbose=0, retpath=False):
"""
Call anaddb to compute elastic and piezoelectric tensors. Require DDB with strain terms.
By default, this method sets the anaddb input variables automatically
by looking at the 2nd-order derivatives available in the DDB file.
This behaviour can be changed by setting explicitly the value of:
`relaxed_ion` and `piezo`.
Args:
relaxed_ion: Activate computation of relaxed-ion tensors.
Allowed values are [True, False, "automatic"]. Defaults to "automatic".
In "automatic" mode, relaxed-ion tensors are automatically computed if
internal strain terms and phonons at Gamma are present in the DDB.
piezo: Activate computation of piezoelectric tensors.
Allowed values are [True, False, "automatic"]. Defaults to "automatic".
In "automatic" mode, piezoelectric tensors are automatically computed if
piezoelectric terms are present in the DDB.
NB: relaxed-ion piezoelectric requires the activation of `relaxed_ion`.
dde: if True, dielectric tensors will be calculated.
stress_correction: Calculate the relaxed ion elastic tensors, considering
the stress left inside cell. The DDB must contain the stress tensor.
asr: Anaddb input variable. See official documentation.
chneut: Anaddb input variable. See official documentation.
mpi_procs: Number of MPI processes to use.
workdir: Working directory. If None, a temporary directory is created.
manager: |TaskManager| object. If None, the object is initialized from the configuration file
verbose: verbosity level. Set it to a value > 0 to get more information
retpath: True to return path to anaddb.nc file.
Return:
|ElasticData| object if `retpath` is None else absolute path to anaddb.nc file.
"""
Let's print the object to get a summary of the most important results:
print(edata)
================================= Structure ================================= Full Formula (Al2 As2) Reduced Formula: AlAs abc : 3.989448 3.989448 6.497130 angles: 90.000000 90.000000 120.000000 Sites (4) # SP a b c --- ---- -------- -------- -------- 0 Al 0.333333 0.666667 0 1 Al 0.666667 0.333333 0.5 2 As 0.333333 0.666667 0.376086 3 As 0.666667 0.333333 0.876086 ============================== Anaddb Variables ============================== { "asr": 2, "chneut": 1, "dieflag": 0, "elaflag": 3, "instrflag": 1, "piezoflag": 3 } ========================= elastic tensors available ========================= [ELASTIC_RELAXED] relaxed-ion elastic tensor in Voigt notation (shape: (6, 6)) Units: GPa, set to zero below: 0.001, fit_to_structure: True xx yy zz yz xz xy xx 135.262182 54.450376 38.052927 0.00000 0.000000 0.000000 yy 54.450376 135.262181 38.052927 0.00000 0.000000 0.000000 zz 38.052927 38.052926 148.211029 0.00000 0.000000 0.000000 yz 0.000000 0.000000 0.000000 30.55071 0.000000 0.000000 xz 0.000000 0.000000 0.000000 0.00000 30.550709 0.000000 xy 0.000000 0.000000 0.000000 0.00000 0.000000 40.405903 [ELASTIC_CLAMPED] clamped-ion elastic tensor in Voigt notation (shape: (6, 6)) Units: GPa, set to zero below: 0.001, fit_to_structure: True xx yy zz yz xz xy xx 165.988592 40.464803 21.090298 0.000000 0.000000 0.000000 yy 40.464802 165.988592 21.090299 0.000000 0.000000 0.000000 zz 21.090298 21.090298 182.585743 0.000000 0.000000 0.000000 yz 0.000000 0.000000 0.000000 40.818194 0.000000 0.000000 xz 0.000000 0.000000 0.000000 0.000000 40.818195 0.000000 xy 0.000000 0.000000 0.000000 0.000000 0.000000 62.761895 ====================== piezoelectric tensors available ====================== [PIEZO_RELAXED] relaxed-ion piezoelectric tensor in Voigt notation (shape: (3, 6)) Units: c/m^2, set to zero below: 1e-05, fit_to_structure: True xx yy zz yz xz xy Px 0.000000 0.000000 0.000000 0.000000 -0.048288 0.0 Py 0.000000 0.000000 0.000000 -0.048288 0.000000 0.0 Pz -0.011872 -0.011872 0.064628 0.000000 0.000000 0.0 [PIEZO_CLAMPED] clamped-ion piezoelectric tensor in Voigt notation (shape: (3, 6)) Units: c/m^2, set to zero below: 1e-05, fit_to_structure: True xx yy zz yz xz xy Px 0.000000 0.000000 0.00000 0.000000 0.435488 0.0 Py 0.000000 0.000000 0.00000 0.435488 0.000000 0.0 Pz 0.384901 0.384901 -0.73943 0.000000 0.000000 0.0
Since the DDB file contains internal strain terms
and piezoeletric terms,
AbiPy set elaflag
to 3, instrflag
to 1 and piezoflag
to 3 so that anaddb
will compute both relaxed
and clamped-ion
tensors.
abilab.docvar("elaflag", executable="anaddb")
0
Flag for calculation of elastic and compliance tensors
The tensors are available as attributes of the edata
object:
edata.elastic_relaxed
xx | yy | zz | yz | xz | xy | |
---|---|---|---|---|---|---|
Voigt index | ||||||
xx | 135.262182 | 54.450376 | 38.052927 | 0.00000 | 0.000000 | 0.000000 |
yy | 54.450376 | 135.262181 | 38.052927 | 0.00000 | 0.000000 | 0.000000 |
zz | 38.052927 | 38.052926 | 148.211029 | 0.00000 | 0.000000 | 0.000000 |
yz | 0.000000 | 0.000000 | 0.000000 | 30.55071 | 0.000000 | 0.000000 |
xz | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 30.550709 | 0.000000 |
xy | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 40.405903 |
These are pymatgen tensors, more specifically ElasticTensor objects so we have access to several useful methods. To get the Voigt bulk modulus, use:
edata.elastic_relaxed.k_voigt
75.53864999272173
while the compliance tensor is easily obtained with:
edata.elastic_relaxed.compliance_tensor
ComplianceTensor([[[[ 9.12540966e-03 -6.63185608e-11 -1.35430963e-11] [-6.63185608e-11 -3.24901993e-03 5.18498777e-13] [-1.35430963e-11 5.18498777e-13 -1.50875299e-03]] [[-4.85950994e-11 6.18721480e-03 2.11517757e-13] [ 6.18721480e-03 7.09227015e-11 -9.82820729e-13] [ 2.11517757e-13 -9.82820729e-13 -1.46075330e-12]] [[-1.49197686e-11 2.11519838e-13 8.18311620e-03] [ 2.11519838e-13 8.81357826e-12 2.09478265e-12] [ 8.18311620e-03 2.09478265e-12 4.38686613e-12]]] [[[-4.85950994e-11 6.18721480e-03 2.11517757e-13] [ 6.18721480e-03 7.09227015e-11 -9.82820729e-13] [ 2.11517757e-13 -9.82820729e-13 -1.46075330e-12]] [[-3.24901988e-03 9.91846920e-11 7.06872590e-12] [ 9.91846920e-11 9.12540970e-03 -7.17330054e-13] [ 7.06872590e-12 -7.17330054e-13 -1.50875296e-03]] [[-1.03689808e-13 -1.01063498e-12 2.09477456e-12] [-1.01063498e-12 -7.20723372e-13 8.18311576e-03] [ 2.09477456e-12 8.18311576e-03 -2.67441612e-13]]] [[[-1.49197686e-11 2.11519838e-13 8.18311620e-03] [ 2.11519838e-13 8.81357826e-12 2.09478265e-12] [ 8.18311620e-03 2.09478265e-12 4.38686613e-12]] [[-1.03689808e-13 -1.01063498e-12 2.09477456e-12] [-1.01063498e-12 -7.20723372e-13 8.18311576e-03] [ 2.09477456e-12 8.18311576e-03 -2.67441612e-13]] [[-1.50875300e-03 -3.19970325e-11 4.99090803e-12] [-3.19970325e-11 -1.50875290e-03 1.92944671e-13] [ 4.99090803e-12 1.92944671e-13 7.52187567e-03]]]])
One can also use the elate online tool to analyse the elastic tensor.
To build a pandas DataFrame with properties derived from the elastic tensor and the associated structure, use:
edata.get_elastic_properties_dataframe(properties_as_index=True)
property | 0 | 1 | |
---|---|---|---|
0 | trans_v | 3.194459e+03 | 3.838052e+03 |
1 | long_v | 5.796035e+03 | 6.295200e+03 |
2 | snyder_ac | 5.767044e+01 | 8.693459e+01 |
3 | snyder_opt | 3.158141e-01 | 3.621134e-01 |
4 | snyder_total | 5.798626e+01 | 8.729670e+01 |
5 | clarke_thermalcond | 7.734051e-01 | 9.006348e-01 |
6 | cahill_thermalcond | 8.539415e-01 | 9.791319e-01 |
7 | debye_temperature | 3.760623e+02 | 4.477874e+02 |
8 | k_voigt | 7.553865e+01 | 7.553930e+01 |
9 | k_reuss | 7.553074e+01 | 7.553579e+01 |
10 | k_vrh | 7.553469e+01 | 7.553754e+01 |
11 | g_voigt | 3.951341e+01 | 5.767416e+01 |
12 | g_reuss | 3.761300e+01 | 5.366047e+01 |
13 | g_vrh | 3.856321e+01 | 5.566731e+01 |
14 | universal_anisotropy | 2.527308e-01 | 3.740360e-01 |
15 | homogeneous_poisson | 2.818554e-01 | 2.041909e-01 |
16 | y_mod | 9.886491e+10 | 1.340681e+11 |
For the meaning of the different quantities please consult the pymatgen module.
To construct a dataframe with the Voigt indices and the tensor elements use:
edata.get_elastic_voigt_dataframe(tol=1e-5)
voigt_cinds | elastic_relaxed | elastic_clamped | |
---|---|---|---|
0 | (0, 0) | 135.262182 | 165.988592 |
1 | (0, 1) | 54.450376 | 40.464803 |
2 | (0, 2) | 38.052927 | 21.090298 |
3 | (0, 3) | 0.000000 | 0.000000 |
4 | (0, 4) | 0.000000 | 0.000000 |
5 | (0, 5) | 0.000000 | 0.000000 |
6 | (1, 0) | 54.450376 | 40.464802 |
7 | (1, 1) | 135.262181 | 165.988592 |
8 | (1, 2) | 38.052927 | 21.090299 |
9 | (1, 3) | 0.000000 | 0.000000 |
10 | (1, 4) | 0.000000 | 0.000000 |
11 | (1, 5) | 0.000000 | 0.000000 |
12 | (2, 0) | 38.052927 | 21.090298 |
13 | (2, 1) | 38.052926 | 21.090298 |
14 | (2, 2) | 148.211029 | 182.585743 |
15 | (2, 3) | 0.000000 | 0.000000 |
16 | (2, 4) | 0.000000 | 0.000000 |
17 | (2, 5) | 0.000000 | 0.000000 |
18 | (3, 0) | 0.000000 | 0.000000 |
19 | (3, 1) | 0.000000 | 0.000000 |
20 | (3, 2) | 0.000000 | 0.000000 |
21 | (3, 3) | 30.550710 | 40.818194 |
22 | (3, 4) | 0.000000 | 0.000000 |
23 | (3, 5) | 0.000000 | 0.000000 |
24 | (4, 0) | 0.000000 | 0.000000 |
25 | (4, 1) | 0.000000 | 0.000000 |
26 | (4, 2) | 0.000000 | 0.000000 |
27 | (4, 3) | 0.000000 | 0.000000 |
28 | (4, 4) | 30.550709 | 40.818195 |
29 | (4, 5) | 0.000000 | 0.000000 |
30 | (5, 0) | 0.000000 | 0.000000 |
31 | (5, 1) | 0.000000 | 0.000000 |
32 | (5, 2) | 0.000000 | 0.000000 |
33 | (5, 3) | 0.000000 | 0.000000 |
34 | (5, 4) | 0.000000 | 0.000000 |
35 | (5, 5) | 40.405903 | 62.761895 |
Note that the Voigt indices are given following the Python (C) notation in which we start to count from zero.
This might be a bit confusing, especially when comparing with results reported in the literature. The reason why we opted with the 0-based notation is that it facilitates the integration between the DataFrame and other python methods. A similar approach is used in AbiPy when e.g. one has to specify the band or the phonon mode index.
At this point, it should be clear how to analyze the relaxed-ion piezoelectric tensor:
edata.piezo_relaxed
xx | yy | zz | yz | xz | xy | |
---|---|---|---|---|---|---|
Voigt index | ||||||
Px | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -0.048288 | 0.0 |
Py | 0.000000 | 0.000000 | 0.000000 | -0.048288 | 0.000000 | 0.0 |
Pz | -0.011872 | -0.011872 | 0.064628 | 0.000000 | 0.000000 | 0.0 |
edata.get_piezo_voigt_dataframe(tol=1e-6)
voigt_cinds | piezo_relaxed | piezo_clamped | |
---|---|---|---|
0 | (0, 0) | 0.000000 | 0.000000 |
1 | (0, 1) | 0.000000 | 0.000000 |
2 | (0, 2) | 0.000000 | 0.000000 |
3 | (0, 3) | 0.000000 | 0.000000 |
4 | (0, 4) | -0.048288 | 0.435488 |
5 | (0, 5) | 0.000000 | 0.000000 |
6 | (1, 0) | 0.000000 | 0.000000 |
7 | (1, 1) | 0.000000 | 0.000000 |
8 | (1, 2) | 0.000000 | 0.000000 |
9 | (1, 3) | -0.048288 | 0.435488 |
10 | (1, 4) | 0.000000 | 0.000000 |
11 | (1, 5) | 0.000000 | 0.000000 |
12 | (2, 0) | -0.011872 | 0.384901 |
13 | (2, 1) | -0.011872 | 0.384901 |
14 | (2, 2) | 0.064628 | -0.739430 |
15 | (2, 3) | 0.000000 | 0.000000 |
16 | (2, 4) | 0.000000 | 0.000000 |
17 | (2, 5) | 0.000000 | 0.000000 |
In this part of the tutorial, we discuss how to compute elastic and piezoelectric
tensors with different k-point meshes and how to use the DdbRobot
to analyze the results.
In principle, one should relax the structure with different k-point samplings and then use the relaxed structure to compute elastic and piezoelectric properties. This is easy with AbiPy but, for the time being, we ignore this point and use the same structure so that we can learn how to use Python to analyze multiple calculations. At the end of this tutorial you will find an example of flow in which the structural relaxation is performed with different k-meshes.
Since we do not have to change the structure, performing a convergence study with respect to k-points
is just a matter of creating multiple Works
inside a loop over k-meshes (our make_scf_input
is already accepting ngkpt
in input):
from lesson_elastic import build_ngkpt_convflow
abilab.print_source(build_ngkpt_convflow)
def build_ngkpt_convflow(options=None, ngkpt_list=([2, 2, 2], [4, 4, 4], [8, 8, 8])):
"""
Build and return a flow computing elastic and piezoelectric properties with
different k-point samplings given in `ngkpt_list`.
In principle, one should perform different structural relaxations for each `ngkpt`
and use the relaxed structures to compute elastic properties.
"""
workdir = options.workdir if (options and options.workdir) else "flow_elastic_ngkpt_conv"
flow = flowtk.Flow(workdir=workdir)
for ngkpt in ngkpt_list:
scf_input = make_scf_input(ngkpt=ngkpt)
elast_work = flowtk.ElasticWork.from_scf_input(scf_input, with_relaxed_ion=True, with_piezo=True)
flow.register_work(elast_work)
return flow
ngkpt_flow = build_ngkpt_convflow(options=None, ngkpt_list=([2, 2, 2], [4, 4, 4], [8, 8, 8]))
ngkpt_flow.get_graphviz()
#ngkpt_flow.get_vars_dataframe("ngkpt")
To generate the flow with the lesson_elastic.py
script, open the file,
comment the call to build_flow
and uncomment build_ngkpt_convflow
.
Then run the script and launch the calculation with:
abirun.py flow_elastic_ngkpt_conv scheduler
as usual.
There are several output files located inside the outdata
directories:
!find flow_elastic_ngkpt_conv/ -name "*_DDB"
flow_elastic_ngkpt_conv//w0/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t0/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t10/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t11/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t12/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t13/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t4/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t5/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t6/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t7/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t8/outdata/out_DDB flow_elastic_ngkpt_conv//w0/t9/outdata/out_DDB flow_elastic_ngkpt_conv//w1/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t0/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t10/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t11/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t12/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t13/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t4/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t5/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t6/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t7/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t8/outdata/out_DDB flow_elastic_ngkpt_conv//w1/t9/outdata/out_DDB flow_elastic_ngkpt_conv//w2/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t0/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t10/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t11/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t12/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t13/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t4/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t5/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t6/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t7/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t8/outdata/out_DDB flow_elastic_ngkpt_conv//w2/t9/outdata/out_DDB
Remember that our goal is to analyze the convergence of the elastic and piezoeletric properties
as function of nkpt
.
So we are mainly interested in the final DDB files located in the outdata
directories
of the works (w0/outdata
, w1/outdata
, w2/outdata
).
These are indeed the DDB files with all the information needed in anaddb to compute the tensors.
The code below tells our robot that we would like to analyze all the DDB files located in the output directories of the works:
robot = abilab.DdbRobot.from_dir_glob("./flow_elastic_ngkpt_conv/w*/outdata/")
robot
The DDB file are available in robot.abifile
.
Each DdbFile
object has a header (dictionary) containing metadata extracted from the DDB file.
To get the keywords in the header, use:
robot.abifiles[0].header.keys()
dict_keys(['version', 'lines', 'usepaw', 'natom', 'nkpt', 'nsppol', 'nsym', 'ntypat', 'occopt', 'nband', 'acell', 'amu', 'dilatmx', 'ecut', 'ecutsm', 'intxc', 'iscf', 'ixc', 'kpt', 'kptnrm', 'ngfft', 'nspden', 'nspinor', 'occ', 'rprim', 'dfpt_sciss', 'spinat', 'symafm', 'symrel', 'tnons', 'tolwfr', 'tphysel', 'tsmear', 'typat', 'wtk', 'xred', 'znucl', 'zion'])
We will be using these metavariables to construct our pandas Dataframe so that we can analyze
the convergence of our physical quantities with e.g. nkpt
.
Let's call anacompare_elastic
to construct a DataFrame (data
) with the elastic properties obtained
with the three DDB files and add the value of ddb.header["nkpt"]
.
elastdata_list
is a list of ElastData
object we can use afterwards to access the individual tensors:
data, elastdata_list = robot.anacompare_elastic(ddb_header_keys="nkpt")
data.keys()
Index(['trans_v', 'long_v', 'snyder_ac', 'snyder_opt', 'snyder_total', 'clarke_thermalcond', 'cahill_thermalcond', 'debye_temperature', 'k_voigt', 'k_reuss', 'k_vrh', 'g_voigt', 'g_reuss', 'g_vrh', 'universal_anisotropy', 'homogeneous_poisson', 'y_mod', 'tensor_name', 'formula', 'nkpt', 'natom', 'alpha', 'beta', 'gamma', 'a', 'b', 'c', 'volume', 'abispg_num', 'spglib_symb', 'spglib_num', 'spglib_lattice_type'], dtype='object')
data
trans_v | long_v | snyder_ac | snyder_opt | snyder_total | clarke_thermalcond | cahill_thermalcond | debye_temperature | k_voigt | k_reuss | ... | beta | gamma | a | b | c | volume | abispg_num | spglib_symb | spglib_num | spglib_lattice_type | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 3151.515931 | 5705.015870 | 55.194905 | 0.311229 | 55.506134 | 0.762578 | 0.841544 | 370.940942 | 73.572019 | 72.331727 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
1 | 3753.376687 | 6179.896889 | 81.728454 | 0.354736 | 82.083190 | 0.882070 | 0.959183 | 438.077881 | 73.667176 | 73.014097 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
2 | 3194.458694 | 5796.034976 | 57.670442 | 0.315814 | 57.986256 | 0.773405 | 0.853941 | 376.062257 | 75.538650 | 75.530735 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
3 | 3838.051894 | 6295.200071 | 86.934586 | 0.362113 | 87.296699 | 0.900635 | 0.979132 | 447.787424 | 75.539303 | 75.535785 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
4 | 3187.357488 | 5761.696011 | 56.983956 | 0.314556 | 57.298512 | 0.770980 | 0.850540 | 375.118016 | 74.267964 | 74.258479 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
5 | 3849.753466 | 6277.945635 | 87.049415 | 0.362273 | 87.411688 | 0.901308 | 0.979563 | 448.886026 | 74.268035 | 74.260840 | ... | 90.0 | 120.0 | 3.989448 | 3.989448 | 6.49713 | 89.552529 | 0 | P6_3mc | 186 | hexagonal |
6 rows × 32 columns
data[["k_voigt", "nkpt", "tensor_name"]]
k_voigt | nkpt | tensor_name | |
---|---|---|---|
0 | 73.572019 | 2 | elastic_relaxed |
1 | 73.667176 | 2 | elastic_clamped |
2 | 75.538650 | 8 | elastic_relaxed |
3 | 75.539303 | 8 | elastic_clamped |
4 | 74.267964 | 40 | elastic_relaxed |
5 | 74.268035 | 40 | elastic_clamped |
To plot the convergence of selected properties versus the number of k-points, use:
robot.plot_xy_with_hue(data, x="nkpt", y=["k_voigt", "g_voigt", "y_mod"], hue="tensor_name");
For more examples on the use of DDB and robots, see the DDB notebook
Now let's try something a bit more complicated.
Assume we want to plot the convergence of the individual elements of the tensor as as function of nkpt
.
In this case, we have to work a bit more to create a Dataframe with the elements in Voigt notation,
add the value of nkpt
associated to this tensor and finally concatenate the results in a single DataFrame:
df_list = []
for edata, ddb in zip(elastdata_list, robot.abifiles):
# Get dataframe with tensor elements in Voigt notation
df = edata.get_elastic_voigt_dataframe()
# Add metadata and store dataframe in df_list
df["nkpt"] = ddb.header["nkpt"]
df_list.append(df)
# Concatenate dataframes
import pandas as pd
data = pd.concat(df_list)
data.head()
voigt_cinds | elastic_relaxed | elastic_clamped | nkpt | |
---|---|---|---|---|
0 | (0, 0) | 1.526887e+02 | 1.778955e+02 | 2 |
1 | (0, 1) | 5.327625e+01 | 3.828480e+01 | 2 |
2 | (0, 2) | 2.952733e+01 | 1.722049e+01 | 2 |
3 | (0, 3) | -1.250329e-09 | 5.096060e-10 | 2 |
4 | (0, 4) | 7.739779e-08 | -2.507582e-09 | 2 |
To select only the (0, 0) elements, use the syntax:
c00 = data[data["voigt_cinds"] == (0, 0)]
c00
voigt_cinds | elastic_relaxed | elastic_clamped | nkpt | |
---|---|---|---|---|
0 | (0, 0) | 152.688699 | 177.895515 | 2 |
0 | (0, 0) | 135.262182 | 165.988592 | 8 |
0 | (0, 0) | 130.988919 | 162.816606 | 40 |
Now we can finally plot the (0, 0) tensor elements as function of nkpt
with:
c00.plot(x="nkpt", y=[k for k in C00 if k.startswith("elastic_")], subplots=True, style="-o");
Back to the main Index