静态非线性性质：

1.首先进行基态计算，收敛测试并优化结构，

2.PEAD方法：

# Linear and nonlinear response calculation for AlP# Perturbations: electric fields & atomic displacementsndtset 5#DATASET1 : scf calculation: GS WF in the IBZ#********************************************   prtden1    1        # save density on disk, will be used in other datasets    prtwf1    1        # save WF on disk, will be used in other datasets   kptopt1    1        # use Irreducible Brillouin Zone (all symmetry taken into account)   toldfe1    1.0d-12  # SCF convergence criteria (could be tolwfr or tolvrs)#DATASET2 : non scf calculation: GS WF in half BZ#*****************************************************   getden2    1        # use density from dataset 1   kptopt2    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   getwfk2    1        # use GS WF from dataset 1 (as input)     iscf2   -2        # non-self-consistent calculation   tolwfr2    1.0d-22  # convergence criteria on WF, need high precision for response    prtwf2    1        # save WF on disk, will be used in other datasets#DATASET3 : derivative of WF with respect to k points (d/dk)#**********************************************************   getwfk3    2        # use GS WF from dataset 2   kptopt3    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   rfelfd3    2        # compute 1st-order WF derivatives (d/dk)    rfdir3    1 1 1    # compute all directions   tolwfr3    1.0d-22  # convergence criteria on 1st-order WF    prtwf3    1        # save 1st-order WF on disk, will be used in other datasets#DATASET4 : response functions (2nd derivatives of E)#           and corresponding 1st order WF derivatives#           phonons, electric fields, and strains are all done#**************************************************************   getwfk4    2        # use GS WF from dataset 2   kptopt4    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   getddk4    3        # use ddk WF from dataset 3 (needed for electric field)   rfphon4    1        # compute 1st-order WF derivatives with respect to atomic displacements...  rfatpol4    1 2      # ...for all atoms (so here only 1 and 2)    rfdir4    1 1 1    # compute all directions   rfelfd4    3        # compute 1st-order WF derivatives with respect to electric field   rfstrs4    3        # compute 1st-order WF derivatives with respect to strains   tolvrs4    1.0d-12  # SCF convergence criteria (could be tolwfr)  prepanl4    1        # make sure that response functions are correctly prepared for a non-linear computation    prtwf4    1        # save 1st-order WF on disk, will be used in other datasets   prtden4    1        # save 1st-order density on disk, will be used in other datasets#DATASET5 : 3rd derivatives of E#*********************************   getwfk5           2     # use GS WF from dataset 2  get1den5           4     # use 1st-order densities from dataset 4   get1wf5           4     # use 1st-order WFs from dataset 4   kptopt5           2     # use Half Brillouin Zone (only time-reversal symmetry taken into account)optdriver5           5     # compute 3rd order derivatives of the energy  d3e_pert1_elfd5    1     # activate electric field for 1st perturbation...  d3e_pert1_phon5    1     # ...and also atomic displacements... d3e_pert1_atpol5    1 2   # ...for all atoms (so here only 1 and 2)...   d3e_pert1_dir5    1 1 1 # ...for all directions (for both atomic displacements and electric fields)  d3e_pert2_elfd5    1     # activate electric field for 2nd perturbation...   d3e_pert2_dir5    1 1 1 # ...for all directions  d3e_pert3_elfd5    1     # activate electric field for 3rd perturbation...   d3e_pert3_dir5    1 1 1 # ...for all directions#Definition of the unit cell# these cell parameters were derived from a relaxation run done with the# current ecut and kpt values. The current ecut value used is very low but# is done to speed the calculations.# acell      7.2728565836E+00  7.2728565836E+00  7.2728565836E+00 Bohr # ecut 5# acell      7.2511099467E+00  7.2511099467E+00  7.2511099467E+00 Bohr # ecut 30acell      7.1391127387E+00  7.1391127387E+00  7.1391127387E+00 Bohr # ecut 2.8rprim      0.0000000000E+00  7.0710678119E-01  7.0710678119E-01           7.0710678119E-01  0.0000000000E+00  7.0710678119E-01           7.0710678119E-01  7.0710678119E-01  0.0000000000E+00#Definition of the atom types and pseudopotentials ntypat 2 # two types of atoms znucl 15 13 # the atom types are Phosphorous and Aluminum#Definition of the atoms natom 2 # two atoms in the cell typat 1 2 # type 1 is Phosphorous, type 2 is Aluminum (order defined by znucl above and pseudos list) nband 4 # nband is restricted here to the number of filled bands only, no empty bands. nbdbuf 0xred1/4 1/4 1/40 0 0#Numerical parameters of the calculation : planewave basis set and k point gridecut  2.8 # this value is very low but is used here to achieve very low calculation times.ecutsm 0.5dilatmx 1.05ngkpt 6 6 6 # must be checked carefully for convergence, Raman calculations converge slowly with kptnshiftk 4 # this Monkhorst-Pack shift pattern is used so that the symmetry of the shifted grid          # is correct. A gamma-centered grid would also have the correct symmetry but would be          # less efficient.shiftk 0.5 0.5 0.5       0.5 0.0 0.0       0.0 0.5 0.0       0.0 0.0 0.5# scf parametersnstep 8  # limit number of steps for the tutorial, in general this should be set to its# suppress printing of density, wavefunctions, etc except what is# explicitly requested above in the ndtset sectionprtwf 0prtden 0prteig 0

包含了电场、原子位移和力的响应，前4个数据是典型的线性响应计算，第5个数据是能量的3阶导数，用于计算非线性响应，例如非线性光学系数以及拉曼张量。强制性的定义：nbdbuf 0 nband 4 (= number of valence bands)，

2.full DFPT计算线性和非线性响应，作为对比前4个数据和上面相似的，第5和6计算了需要的二阶波函数，和数据3相似，

# Linear and nonlinear response calculation for AlP# Perturbations: electric fields & atomic displacements# Use of 'full DFPT' method for third derivatives,# so we need to solve Second Order Sternheimer equationndtset 7#DATASET1 : scf calculation: GS WF in the IBZ#********************************************   prtden1    1        # save density on disk, will be used in other datasets    prtwf1    1        # save WF on disk, will be used in other datasets   kptopt1    1        # use Irreducible Brillouin Zone (all symmetry taken into account)   toldfe1    1.0d-12  # SCF convergence criteria (could be tolwfr or tolvrs)#DATASET2 : non scf calculation: GS WF in half BZ#*****************************************************   getden2    1        # use density from dataset 1   kptopt2    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   getwfk2    1        # use GS WF from dataset 1 (as input)     iscf2   -2        # non-self-consistent calculation   tolwfr2    1.0d-22  # convergence criteria on WF, need high precision for response    prtwf2    1        # save WF on disk, will be used in other datasets#DATASET3 : derivative of WF with respect to k points (d/dk)#**********************************************************   getwfk3    2        # use GS WF from dataset 2   kptopt3    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   rfelfd3    2        # compute 1st-order WF derivatives (d/dk)    rfdir3    1 1 1    # compute all directions   tolwfr3    1.0d-22  # convergence criteria on 1st-order WF    prtwf3    1        # save 1st-order WF on disk, will be used in other datasets#DATASET4 : response functions (2nd derivatives of E)#           and corresponding 1st order WF derivatives#           phonons, electric fields, and strains are all done#**************************************************************   getwfk4    2        # use GS WF from dataset 2   kptopt4    2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   getddk4    3        # use ddk WF from dataset 3 (needed for electric field)   rfphon4    1        # compute 1st-order WF derivatives with respect to atomic displacements...  rfatpol4    1 2      # ...for all atoms (so here only 1 and 2)    rfdir4    1 1 1    # compute all directions   rfelfd4    3        # compute 1st-order WF derivatives with respect to electric field   rfstrs4    3        # compute 1st-order WF derivatives with respect to strains   tolvrs4    1.0d-12  # SCF convergence criteria (could be tolwfr)  prepanl4    1        # make sure that response functions are correctly prepared for a non-linear computation    prtwf4    1        # save 1st-order WF on disk, will be used in other datasets   prtden4    1        # save 1st-order density on disk, will be used in other datasets# DATASET5 : 2nd order WF derivatives with respect to k-points (d/dkdk)#********************************************************************** rf2_dkdk5   1        # compute 2nd-order WF derivatives (d/dkdk)   getwfk5   2        # use GS WF from dataset 2   kptopt5   2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)   getddk5   3        # use 1st-order WF from dataset 3  prepanl5   1        # compute only needed directions for the non-linear computation   tolwfr5   1.0d-22  # convergence criteria on 2nd-order WF    prtwf5   1        # save 2nd-order WF on disk, will be used in other datasets# DATASET6 : 2nd order WF derivatives with respect to electric field and k-points (d/dkde)#***************************************************************************************** rf2_dkde6   1        # compute 2nd-order WF derivatives (d/dkde)   getwfk6   2        # use GS WF from dataset 2   getddk6   3        # use 1st-order WF (d/dk) from dataset 3  get1den6   4        # use 1st-order densities from dataset 4 getdelfd6   4        # use 1st-order WF (d/de) from dataset 4  getdkdk6   5        # use 2nd-order WF (d/dkdk) from dataset 5  prepanl6   1        # compute only directions needed for nonlinear computation   tolwfr6   1.0d-22  # convergence criteria on 2nd-order WF   kptopt6   2        # use Half Brillouin Zone (only time-reversal symmetry taken into account)    prtwf6   1        # save 2nd-order WF on disk, will be used in other datasets#DATASET7 : 3rd derivatives of E#*********************************   getwfk7           2     # use GS WF from dataset 2   getddk7           3     # use 1st-order WF (d/dk) from dataset 3  get1den7           4     # use 1st-order densities from dataset 4   get1wf7           4     # use 1st-order WFs from dataset 4  getdkde7           6     # use 2nd-order WF (d/dkde) from dataset 6   kptopt7           2     # use Half Brillouin Zone (only time-reversal symmetry taken into account)optdriver7           5     # compute 3rd order derivatives of the energy  usepead7           0     # no use of pead method, so use full DFPT instead  d3e_pert1_elfd7    1     # activate electric field for 1st perturbation...  d3e_pert1_phon7    1     # ...and also atomic displacements... d3e_pert1_atpol7    1 2   # ...for all atoms (so here only 1 and 2)...   d3e_pert1_dir7    1 1 1 # ...for all directions (for both atomic displacements and electric fields)  d3e_pert2_elfd7    1     # activate electric field for 2nd perturbation...   d3e_pert2_dir7    1 1 1 # ...for all directions  d3e_pert3_elfd7    1     # activate electric field for 3rd perturbation...   d3e_pert3_dir7    1 1 1 # ...for all directions#Definition of the unit cell# these cell parameters were derived from a relaxation run done with the# current ecut and kpt values. The current ecut value used is very low but# is done to speed the calculations.# acell      7.2728565836E+00  7.2728565836E+00  7.2728565836E+00 Bohr # ecut 5# acell      7.2511099467E+00  7.2511099467E+00  7.2511099467E+00 Bohr # ecut 30acell      7.1391127387E+00  7.1391127387E+00  7.1391127387E+00 Bohr # ecut 2.8rprim      0.0000000000E+00  7.0710678119E-01  7.0710678119E-01           7.0710678119E-01  0.0000000000E+00  7.0710678119E-01           7.0710678119E-01  7.0710678119E-01  0.0000000000E+00#Definition of the atom types and pseudopotentials ntypat 2 # two types of atoms znucl 15 13 # the atom types are Phosphorous and Aluminum#Definition of the atoms natom 2 # two atoms in the cell typat 1 2 # type 1 is Phosphorous, type 2 is Aluminum (order defined by znucl above and pseudos list) nband 4 # nband is restricted here to the number of filled bands only, no empty bands. nbdbuf 0# atomic positions in units of cell vectorsxred1/4 1/4 1/40 0 0#Numerical parameters of the calculation : planewave basis set and k point gridecut  2.8 # this value is very low but is used here to achieve very low calculation times.        # in a production environment this should be checked carefully for convergence and        # a more reasonable value is probably around 30ecutsm 0.5dilatmx 1.05ngkpt 6 6 6 # must be checked carefully for convergence, Raman calculations converge slowly with kptnshiftk 4 # this Monkhorst-Pack shift pattern is used so that the symmetry of the shifted grid          # is correct. A gamma-centered grid would also have the correct symmetry but would be          # less efficient.shiftk 0.5 0.5 0.5       0.5 0.0 0.0       0.0 0.5 0.0       0.0 0.0 0.5# scf parametersnstep 8  # limit number of steps for the tutorial, in general this should be set to its         # default (30) or higher# suppress printing of density, wavefunctions, etc except what is# explicitly requested above in the ndtset sectionprtwf 0prtden 0prteig 0

3.合并DDB：用上面任意2种产生的输出文件处理，例如PEAD方法，产生了4和5两个DDB文件，

tnlo_3.ddb.outlinear + nonlinear response calculation of AlP2tnlo_2o_DS4_DDBtnlo_2o_DS5_DDB

4.分析DDB：files:

tnlo_4.abitnlo_4.abotnlo_3.ddb.outtnlo_4_thm_dummytnlo_4_gkk_dummytnlo_4_ep_dummytnlo_4_ddk_dummy

输入文件：

# General information#**********************   rfmeth  1   enunit  0   eivec   1   chneut  2   selectz 0# Flags#*******   dieflag  1   ifcflag  0 ! This is the default value   thmflag  0 ! This is the default value    nlflag  1   elaflag  3 piezoflag  3 instrflag  1# Nonlinear response calculation#*******************************   ramansr  1    alphon  1    prtmbm  1# Wavevector list number 1**************************                                                           nph1l  1     qph1l  0.0 0.0 0.0 1.0 # Wave vector list no. 2#***********************     nph2l  1     qph2l  1.0 0.0 0.0 0.0

其中应变的微扰：elaflag, piezoflag, and instrflag；非线性响应需要设置：

 nlflag  1     #打开非线性响应ramansr 1    #原子位移下电子介电常数的一阶变化 alphon  1     prtmbm  1   #获得离子对电光系数的声子贡献

Wavevector用于获得拉曼系数，运行后在输出文件中可以找到以下物理参数：

Al的波恩有效电荷：

Z*_Al = 2.207969 (2.233928)

Gamma点的光子声子频率（Optical phonon frequencies ）：

w_TO (cm^-1) = 463.2713 (417.4934)w_LO (cm^-1) = 534.6882 (493.4411)

线性光学介电常数：

Electronic dielectric tensor = 6.12216718 (6.10643103)

静态介电常数：

Relaxed ion dielectric tensor = 8.15521897 (8.53019279)

压电系数：

clamped ion (C/m^2) = -0.58864089 (-0.58710170)relaxed ion (C/m^2) =  0.26741085 ( 0.11459002)

非线性光学系数：

d_36 (pm/V)  = 8.934453 (10.174996) [21]

电光系数：

Electronic EO constant (pm/V): -0.953493194  (-1.091488899)Full Ionic EO constant (pm/V):  0.536131045  ( 0.662692165)     Total EO constant (pm/V): -0.417362150  (-0.428796933)

拉曼性质，横向和纵向：

alpha(TO) = 0.004315217 (0.004379774)alpha(LO) = 0.006863765 (0.007243040)

d chi_23/d tau_1 (Bohr^-1, Al)  = -0.043617186 (-0.043054192)

以上（）括号内数字是更高ecut下的计算结果.

nband的值是价带填满的值，在使用的赝势文件中可以找到，要注意乘以原子个数，PEAD算法不兼容PAW赝势；（usepead=1）教程的例子是半导体或绝缘体，对于金属不太适合。

*长路漫漫修远兮，一起学习与进步。学到新东西会继续更新。