特性阻抗的元参数：Abinit-11:介电张量声子频率和波恩有效电荷的计算

特性阻抗的元参数：Abinit-11:介电张量声子频率和波恩有效电荷的计算5.DFPT均匀电场的计算：（此步为重点，使用了多数据模式，以上步皆为预热和加深理解） Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000 Phonon energies in Hartree : 2.558878E-06 2.558879E-06 2.558880E-06 1.568559E-03 1.568559E-03 1.568559E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616089E-01 5.616092E-01 3.442590E 02 3.442590E 02 - 3.442590E 02好结果是有三个声子频率小于1，坏结果是其他三个声子频率与实验和理论不符，这是由于忽略了原子位移和

接着上一节，我们讨论DFPT的计算

1. 自洽计算 目的是获得电荷密度和波函数，tolvrs精度非常高，是为了获取较好的波函数，后续使用

#Specific to ground state calculation tolvrs 1.0d-18 # SCF stopping criterion ####################################################################### #Common input variables #Definition of the unit cell acell 3*10.61 # This is equivalent to 10.61 10.61 10.61 rprim 0.0 0.5 0.5 # In tutorials 1 and 2 these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is the default. #Definition of the atom types ntypat 2 # There are two types of atom znucl 13 33 # The keyword "znucl" refers to the atomic number #Definition of the atoms natom 2 # There are two atoms typat 1 2 # The first is of type 1 (Al) the second is of type 2 (As). xred # This keyword indicate that the location of the atoms 0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1. 0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2. #Gives the number of band explicitely (do not take the default) nband 4 # For an insulator #Definition of the planewave basis set ecut 3.0 # Maximal kinetic energy cut-off in Hartree #Definition of the k-point grid kptrlatt -4 4 4 # In cartesian coordinates this grid is simple cubic and 4 -4 4 # actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid. 4 4 -4 # It might as well be obtained through the use of ngkpt nshiftk and shiftk . #Definition of the SCF procedure nstep 15 # Maximal number of SCF cycles diemac 9.0 # it is worth to precondition the SCF cycle. # The dielectric constant of AlAs is smaller that the one of Si (=12).2

2.把输出的波函数文件中o_WFK修改为i_WFK，以便读取，沿着Al原子的x轴位移0.001作为一个扰动，输入文件如下：

# Crystalline AlAs : computation of the total energy # and forces in a distorted geometry irdwfk 1 # Read input waveFunctions #Specific to ground state calculation tolvrs 1.0d-12 # SCF stopping criterion ####################################################################### #Common input variables #Definition of the unit cell acell 3*10.61 # This is equivalent to 10.61 10.61 10.61 rprim 0.0 0.5 0.5 # In tutorials 1 and 2 these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is the default. #Definition of the atom types ntypat 2 # There are two types of atom znucl 13 33 # The keyword "znucl" refers to the atomic number of the #Definition of the atoms natom 2 # There are two atoms typat 1 2 # The first is of type 1 (Al) the second is of type 2 (As). xred # This keyword indicate that the location of the atoms will follow 0.001 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1. 0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2. #Gives the number of band explicitely (do not take the default) nband 4 # For an insulator there is no need to include conduction bands in response-function calculations #Definition of the planewave basis set ecut 3.0 # Maximal kinetic energy cut-off in Hartree #Definition of the k-point grid kptrlatt -4 4 4 # In cartesian coordinates this grid is simple cubic and 4 -4 4 # actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid. 4 4 -4 # It might as well be obtained through the use of ngkpt nshiftk and shiftk . #Definition of the SCF procedure nstep 15 # Maximal number of SCF cycles diemac 9.0 # Although this is not mandatory it is worth to

由于对称性降低，可以看到输出文件中nsym变小，K点增多，从输出文件中还可以看到总能的变化以及总能与改变坐标的梯度值dE/dt：

rms dE/dt= 3.5517E-03; max dE/dt= 5.0079E-03; dE/dt below (all hartree) 1 0.005007937776 0.002526310510 0.002526310510 2 -0.005007879256 -0.002526283046 -0.002526283046 ... total_energy : -9.76268124105730E 00

3.总能的二阶导数计算，记得把第2步输出的波函数文件o_WFK修改为i_WFK文件名，输入文件中加入了一些响应函数计算变量：

#Response-function calculation with q=0 rfphon 1 # 声子的响应函数，=1 rfatpol 1 1 # 原子极化的响应函数，定义将要位移的原子，1~natom rfdir 1 0 0 # 响应函数的方向，3个轴 nqpt 1 # One wavevector is to be considered，q点的数量 qpt 0 0 0 # This wavevector is q=0 (Gamma) kptopt 2 # Automatic generation of k points taking into account the time-reversal symmetry only tolvrs 1.0d-8 # SCF stopping criterion irdwfk 1 # Read the ground-state wavefunctions ####################################################################### #Common input variables #Definition of the unit cell acell 3*10.61 # This is equivalent to 10.61 10.61 10.61 rprim 0.0 0.5 0.5 # In tutorials 1 and 2 these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is the default. #Definition of the atom types ntypat 2 # There are two types of atom znucl 13 33 # The keyword "znucl" refers to the atomic number of the #Definition of the atoms natom 2 # There are two atoms typat 1 2 # The first is of type 1 (Al) the second is of type 2 (As). xred # This keyword indicate that the location of the atoms 0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1. 0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2. #Gives the number of band explicitely (do not take the default) nband 4 # For an insulator (if described correctly as an insulator #Definition of the planewave basis set ecut 3.0 # Maximal kinetic energy cut-off in Hartree #Definition of the k-point grid kptrlatt -4 4 4 # In cartesian coordinates this grid is simple cubic and 4 -4 4 # actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid. 4 4 -4 # It might as well be obtained through the use of #Definition of the SCF procedure nstep 15 # Maximal number of SCF cycles diemac 9.0 # Although this is not mandatory

从输出文件中可查询到2DEtotal的值。输出文件中还包含DDB文件，输出文件的内容诠释参加abinit官网：https://docs.abinit.org/guide/respfn/#ddb。使用Abipy的命令如下，可检查精度收敛情况

abiopen.py trf1_3.abo --expose -sns=talk

4.在Gamma点（q=0）计算动力学矩阵：读取第二步产生的波函数，input中修改rfatplo和rfdir变量：所有原子在所有方向上移动

# Crystalline AlAs : computation of the dynamical matrix at Gamma #Response-function calculation with q=0 rfphon 1 # Will consider phonon-type perturbation rfatpol 1 2 # All the atoms will be displaced rfdir 1 1 1 # Along all reduced coordinate axis nqpt 1 # One wavevector is to be considered qpt 0 0 0 # This wavevector is q=0 (Gamma) kptopt 2 # Automatic generation of k points taking into account the time-reversal symmetry only tolvrs 1.0d-8 # SCF stopping criterion irdwfk 1 # Read the ground-state wavefunctions ####################################################################### #Common input variables #Definition of the unit cell acell 3*10.61 # This is equivalent to 10.61 10.61 10.61 rprim 0.0 0.5 0.5 # In tutorials 1 and 2 these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is the default. #Definition of the atom types ntypat 2 # There are two types of atom znucl 13 33 # The keyword "znucl" refers to the atomic number of the #Definition of the atoms natom 2 # There are two atoms typat 1 2 # The first is of type 1 (Al) the second is of type 2 (As). xred # This keyword indicate that the location of the atoms 0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1. 0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2. #Gives the number of band explicitely (do not take the default) nband 4 # For an insulator #Definition of the planewave basis set ecut 3.0 # Maximal kinetic energy cut-off in Hartree #Definition of the k-point grid kptrlatt -4 4 4 # In cartesian coordinates this grid is simple cubic and 4 -4 4 # actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid. 4 4 -4 # It might as well be obtained through the use of ngkpt nshiftk and shiftk . #Definition of the SCF procedure nstep 15 # Maximal number of SCF cycles diemac 9.0

关于对称性的设置kptopt变量：基态计算=1；q=0响应函数对应kptopt=2; q不为0时的响应函数计算kptopt=3. 在输出文件中找到声子频率结果；

Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000 Phonon energies in Hartree : 2.558878E-06 2.558879E-06 2.558880E-06 1.568559E-03 1.568559E-03 1.568559E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616089E-01 5.616092E-01 3.442590E 02 3.442590E 02 - 3.442590E 02

好结果是有三个声子频率小于1，坏结果是其他三个声子频率与实验和理论不符，这是由于忽略了原子位移和电场的耦合，因此需要处理均匀电场的微扰：

5.DFPT均匀电场的计算：（此步为重点，使用了多数据模式，以上步皆为预热和加深理解）

rfelfd: 电场的响应变量

# Crystalline AlAs : computation of the response to homogeneous # electric field and atomic displacements at q=0 ndtset 3 #Ground state calculation kptopt1 1 # Automatic generation of k points taking into account the symmetry tolvrs1 1.0d-18 # SCF stopping criterion #Response Function calculation : d/dk rfelfd2 2 # Activate the calculation of the d/dk perturbation rfdir2 1 0 0 # Need to consider the perturbation in the x-direction only # This is rather specific due to the high symmetry of the AlAs crystal # In general just use rfdir 1 1 1 nqpt2 1 qpt2 0.0 0.0 0.0 # This is a calculation at the Gamma point getwfk2 -1 # Uses as input the output wf of the previous dataset kptopt2 2 # Automatic generation of k points using only the time-reversal symmetry to decrease iscf2 -3 # The d/dk perturbation must be treated in a non-self-consistent way tolwfr2 1.0d-22 # Must use tolwfr for non-self-consistent calculations #Here the value of tolwfr is very low. #Response Function calculation : electric field perturbation and phonons rfphon3 1 # Activate the calculation of the atomic dispacement perturbations rfatpol3 1 2 # All the atoms will be displaced rfelfd3 3 # Activate the calculation of the electric field perturbation rfdir3 1 1 1 # All directions are selected. However symmetries will be used to decrease # the number of perturbations so only the x electric field is needed # (and this explains why in the second dataset rfdir was set to 1 0 0). nqpt3 1 qpt3 0.0 0.0 0.0 # This is a calculation at the Gamma point getwfk3 -2 # Uses as input wfs the output wfs of the dataset 1 getddk3 -1 # Uses as input ddk wfs the output of the dataset 2 kptopt3 2 # Automatic generation of k points # using only the time-reversal symmetry to decrease # the size of the k point set. tolvrs3 1.0d-8 ####################################################################### #Common input variables #Definition of the unit cell acell 3*10.61 # This is equivalent to 10.61 10.61 10.61 rprim 0.0 0.5 0.5 # In tutorials 1 and 2 these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is the default. #Definition of the atom types ntypat 2 # There are two types of atom znucl 13 33 # The keyword "znucl" refers to the atomic number of the. #Definition of the atoms natom 2 # There are two atoms typat 1 2 # The first is of type 1 (Al) the second is of type 2 (As). xred 0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1. 0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2. #Gives the number of band explicitely (do not take the default) nband 4 #Definition of the planewave basis set ecut 3.0 # Maximal kinetic energy cut-off in Hartree #Definition of the k-point grid kptrlatt -4 4 4 # In cartesian coordinates this grid is simple cubic and 4 -4 4 # actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid. 4 4 -4 # It might as well be obtained through the use of ngkpt nshiftk and shiftk . #Definition of the SCF procedure nstep 15 # Maximal number of SCF cycles diemac 9.0

包含了电场和原子位移扰动，在第二个数据输出文件中可以看到：显示了一个微扰

==> initialize data related to q vector <== The list of irreducible perturbations for this q vector is: 1) idir= 1 ipert= 3 ================================================================================ -------------------------------------------------------------------------------- Perturbation wavevector (in red.coord.) 0.000000 0.000000 0.000000 Perturbation : derivative vs k along direction 1

以及f-sum接近1：

dfpt_looppert : ek2= 1.6833336546E 01 f-sum rule ratio= 1.0028582975E 00

在第三个数据结果中，考虑了3个扰动：

==> initialize data related to q vector <== The list of irreducible perturbations for this q vector is: 1) idir= 1 ipert= 1 2) idir= 1 ipert= 2 3) idir= 1 ipert= 4

然后给出了介电张量：由于材料的各向同性，介电常数为9.7606052146

Dielectric tensor in cartesian coordinates j1 j2 matrix element dir pert dir pert real part imaginary part 1 4 1 4 9.7606052146 -0.0000000000 1 4 2 4 0.0000000000 -0.0000000000 1 4 3 4 0.0000000000 -0.0000000000 2 4 1 4 0.0000000000 -0.0000000000 2 4 2 4 9.7606052146 -0.0000000000 2 4 3 4 0.0000000000 -0.0000000000 3 4 1 4 0.0000000000 -0.0000000000 3 4 2 4 0.0000000000 -0.0000000000 3 4 3 4 9.7606052146 -0.0000000000

其次给出了波恩有效电荷，有电场扰动获得的以及声子扰动获得的，Al原子是2.104，As原子是-2.127，未完全满足电中性法则，可以通过增加ecut，改善k点收敛，使正负电荷更接近0。最后，给出了声子频率：

Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000 Phonon energies in Hartree : 2.558878E-06 2.558879E-06 2.558880E-06 1.568559E-03 1.568559E-03 1.568559E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616089E-01 5.616092E-01 3.442590E 02 3.442590E 02 - 3.442590E 02 Phonon at Gamma with non-analyticity in the direction (cartesian coordinates) 1.00000 0.00000 0.00000 Phonon energies in Hartree : 2.558878E-06 2.558879E-06 4.068910E-06 1.568559E-03 1.568559E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616089E-01 8.930225E-01 3.442590E 02 3.442590E 02 - 3.795979E 02 Phonon at Gamma with non-analyticity in the direction (cartesian coordinates) 0.00000 1.00000 0.00000 Phonon energies in Hartree : 2.558878E-06 2.558880E-06 4.068909E-06 1.568559E-03 1.568559E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616092E-01 8.930223E-01 3.442590E 02 3.442590E 02 - 3.795979E 02 Phonon at Gamma with non-analyticity in the direction (cartesian coordinates) 0.00000 0.00000 1.00000 Phonon energies in Hartree : 2.558879E-06 2.558880E-06 4.068909E-06 1.568559E-03 1.568559E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.616089E-01 5.616092E-01 8.930223E-01 3.442590E 02 3.442590E 02 - 3.795979E 02

其中包含了4个部分，第1个部分未考虑电场，剩余3个部分考虑了沿着不同轴向的电场，四部分数据相似，说明电场对当前材料影响不大。以上结果数据与实验的差异性可以通过改善ecut，kpoint和赝势来改善数据。

*前路漫漫，仅供参考。