STATIC DIELECTRIC TENSOR (including local field effects in DFT)``` 电子贡献或者
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
```MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))``` ```MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION``` 离子贡献 ```BORN EFFECTIVE CHARGES (in e, cummulative output)``` ```ELASTIC MODULI IONIC CONTR (kBar)``` ```PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)``` ### 频率依赖介电函数：读取WAVECAR IPA： ```javascript {.line-numbers} ALGO = Exact NBANDS = 64 LOPTICS = .TRUE. CSHIFT = 0.100 NEDOS = 2000 #and you might try with the following #LPEAD = .TRUE.
frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects) and frequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)
上面计算完介电常数后，一般正交晶系。 grep "PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)" OUTCAR -A10 & grep "PIEZOELECTRIC TENSOR for field in x, y, z (C/m^2)" OUTCAR -A10
$E = 6062.434883706029P^6 + 2437.756598493882P^4 + -340.5806845162749*P^2 + 10.339600058315137$ 其中参数 A, B, C 分别为
$A = -0.68116 eV *(m^2/C)^2 B = 9.75103 eV *(m^2/C)^4 C = 36.37461 eV *(m^2/C)^6$
绘制能量曲线的脚本
{.line-numbers}
1 2 3 4 5 6 7 8 9 10 11
withopen('energy.dat','r') asf: content = f.readlines() res = [float(i.strip('\n')) for i in content] minres = min(res) res = [1e3*(i-minres) for i in res] import matplotlib.pyplotas plt plt.figure() plt.plot(res,'b.-') plt.xlabel('displacement') plt.ylabel('Energy ($meV$)') plt.show()
### 考虑极化量子 c = 4.0335000000000001 tmp = [] N = 20 for j inrange(N): start = -N/2+j tmp1 = [i+start*c for i in res] # print(tmp1[0]) tmp.append(tmp1) plt.figure() for i intmp: plt.plot(i,'.')
The vasp wiki tells that one can set manual uccupancy with FERWE and FERDO:
One has to converge a WAVECAR before proceeding further.
Let say your system contain 2042 electrons. Assume the ground state is a closed shell singlet. It looks enough to do the calculation without any spin polarization with ISPIN=1 since both the ground and the excited state must have the same spin singlet wavefunction. But in order to move an electron from the HOMO to the LUMO one has to preconverge the ground state with ISPIN=2.
To do an excitation, one has to promote an electron from spatial orbital #1021 to #1022. To do so, one has to
ISTART = 1
ISTART=2 also might work, but ISTART =1 is not a simple vasp continuation. ISTART=1 will transform the old WAVECAR into the new job’s bases, or truncate the bands, or create new ones if the continuation job have more NBANDS setting than the original. Keep in mind NBANDS sighly varies with the number of MPI threads on which VASP is being executed, since NBANDS/THREADS must be an integer. VASP will automatically increase this value accordingly regardless your NBANDS is fixed.
ISMEAR = -2
ISPIN = 2 # 1020 1021 1022 1023 just a commend line
FERWE = 10201.0 10.0 11.0 10000.0 #The syntax is BANDS*OCC
FERDO = 10201.0 11.0 10.0 10000.0
The 1000*0.0 means zero electron for the next bands. Note this will do not override the NBANDS = 1150 setting is this will dont calculate 2022 bands, the bands after 1150 defined with FERWE and FERDO will be simply ignored.
The following SCF cycle alterations also required for HSE06 deltaSCF calculations:
ALGO = All TIME = 0.4 LDIAG = .FALSE. LSUBROT = .FALSE.
shf3 the shellpack have the option to copy the WAVECAR along the with the inputs. OTHER=”diamSiV.WAVECAR.gz” Just uncomment the OTHER in your *.guide file!
#!/bin/bash path=`pwd` #记录当前文件夹位置并定义到变量 path中，称为初始目录 pwd for cij in`ls -F | grep /$` #ls -F| grep /$ 命令为返回当前文件夹中所有子目录 #将所有子目录添加到变量cij中，并准备进行for循环 do cd ${path}/$cij #进入cij变量所代表的子目录中 pwd for s in strain_* #将cij代表的子目录中所有strain文件夹添加到变量s中 #此部添加变量s的方式可更改为添加变量cij的方式，但不建议这么做 do cd ${path}/$cij/$s #进入初始目录下变量cij和变量s组成的目录，由于for循环的设置，将遍历所有可计算的strain文件夹 pwd mpirun -np 16 vasp_std #执行vasp计算 cd ${path} #返回初始目录 pwd done done
>dfc2 DFC2--a generator of renormalized harmonic FCs from SCPH outputs. XML file containing original Fc2 : STo222.xml Output xml filename with anharmonic correction : ST0222_SCPH2-2_300K.xml FC2 correction file fromSCPH calculation :STo_scph2-2.scph_dfc2 Target temperature : 300
New XML file STo222_SCPH2-2_300K.xml was created successfully.
&cell ## primitive cell is best 6.2507550151332 -0.5 0.5 0.5 0.5 -0.5 0.5 0.5 0.5 -0.5 /
&kpoint 1 # KPMODE = 1: get phonon band 2: get others properties G 0.0 0.0 0.0 H 0.5 -0.5 0.5 201 H 0.5 -0.5 0.5 P 0.25 0.25 0.25 201 P 0.25 0.25 0.25 G 0.0 0.0 0.0 201 G 0.0 0.0 0.0 N 0.0 0.0 0.5 201 /
[scf1355@ln151%nc-l ah]$ ~/soft/alamode/tools/dfc2 DFC2 -- a generator of renormalized harmonic FCs from SCPH outputs. XML file containing original FC2 : crte2-442.xml Output xml filename with anharmonic correction : crte2-new.xml FC2 correction file from SCPH calculation : crte2-442.scph_dfc2 Target temperature : 300
The accuracy of the anharmonic force constants can be checked by increasing the size of the training dataset. If the physical properties, such as thermal conductivity and free energy, change significantly with the size of the training dataset, it indicates that the training is insufficient. If the change is minor, your force constants may be OK.
报错：Warning: The following force constant doesn’t exist in the original file: You can neglect this warning as long as the correction term (the last column) is small. This warning is raised when the absolute value of the correction is larger than 1.0e-10 and the corresponding force constant does not exist in the original FCSXML file. In your case, the correction term is only slightly above the tolerance. So, it should not cause any problems hopefully. When the correction term is much larger than 1.0e-10, this warning indicates a possible error in the inputs. Likely errors are the input structure is not fully symmetrized, or the numerical digits of lattice vectors and fractional coordinates are insufficient. the given KMESH_INTERPOLATE is not commensurate with the supercell size of the original FCSXML.
In the current implementation of alamode, only the former choice, renormalized 2FCs + conventional (old) 3FCs, is supported. Another possible choice may be renormalized 2FCs + renormalized 3FCs, but it is not supported within alamode yet. When MODE = RTA, the anphon code only reads the 2FCs and 3FCs from given XML files. By contrast, when MODE = SCPH, the code additionally reads the 4FCs because they are necessary to compute the renormalized 2FCs based on the self-consistent phonon theory. Therefore, the second combination of renormalized 2FCs + high order FCs using SCPH is not possible.
I performed the SCPH calculation. It is not converged at low temperatures (0K, 100K, 200K), while it achieved convergence at high temperatures (300~1200K).A potential solution is to use a smaller MIXALPHA and a smaller DT.
I would like to use second-order, third-order, and fourth-order force constants to obtain anharmonic potentials and calculate the potential energy as a function of atomic displacement. My goal is to generate a plot similar to the one in the paper “Orbitally driven giant phonon anharmonicity in SnSe” (Nature Phys. 2015).Please use https://github.com/ttadano/alamode/blob/develop/tools/calcpes.py.
now that I have obtained the renormalized force constant file xxx.xml at a finite temperature, I now want to convert the renormalized force constant file directly to FORCE_CONSTANS
How to convert a second-order force constant matrix obtained by almode to FORCE_CONSTANTS in phonopy format？
python convert.py alamode.xml > FORCE_CONSTANTS
Could you give me some suggestion? How shoud I covert this with phonopy?
FULL_FORCE_CONSTANTS = .TRUE. WRITE_FORCE_CONSTANTS = .TRUE. Note: the old FORCE_CONSTANTS will be overwrited
How small the fitting error should be?
It depends on the Taylor expansion potential and displacement magnitude you choose. In the standard harmonic calculation where –mag=0.01 is used in displace.py and the all harmonic interactions are considered (cutoff = None), the fitting error is usually less than 5% (~1–2% in most cases). In the calculation of third-order force constants with –mag=0.04, the fitting error should be small as well. Indeed, in many cases, we obtain much smaller fitting errors (< 1%) than the harmonic case. In the temperature-dependent effective potential method, we try to fit the harmonic potential to the displacement-force datasets sampled by ab initio molecular dynamics at finite temperature. Therefore, the fitting error tends to be much larger (> 10%).
h_ref=`grep " $homo_ref " $j | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'` l_ref=`grep " $lumo_ref " $j | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'` echo "The band gap (in eV) without zero-point vibrations is:" echo "$h_ref $l_ref" |awk '{print ($2-$1)}' echo "The band gap (in eV) including zero-point vibrations is:" echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}' echo "The zero-point renormalization of the band gap (in eV) is:" echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $h_ref $l_ref" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-($11-$10)}'
输出如下：
1 2 3 4 5 6
The band gap (in eV) without zero-point vibrations is: 4.4049 The band gap (in eV) including zero-point vibrations is: 4.05102 The zero-point renormalization of the band gap (in eV) is: -0.353883
bandshift=0 gwrun=-1 dgbd=-1 val=-1 con=-1 gwldadiff=-1 test=-1 while [[ $# -gt 0 ]] do key="$1" case $key in esac shift done
if [ -f "helpscript.perl" ]; then rm helpscript.perl fi
cat > helpscript.perl <<EOF #!/bin/perl
use strict; use warnings;
my \$zahler=0; my @entropy=0; my \$ezp=0; my \$fhelmholtz=0; my \$uenergy; my \$kboltzmann=8.6173303*10**(-5.0); my \$ntemp=8; my \$tmax=700; my \$tmin=0; my \$tstep=100; for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { \$entropy[\$itemp]=0; } while (<>) { chomp; \$_=~s/^/ /; my @help=split(/[\t,\s]+/); \$zahler=\$zahler+1; if (\$zahler == 1) {\$uenergy=\$help[1];} else { my \$homega=\$help[2]/1000; \$ezp=\$ezp+\$homega*0.5; for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { my \$temp=(\$itemp-1)*\$tstep+\$tmin; my \$kbt=\$kboltzmann*\$temp; if (\$temp < 0.0000001) { } else { \$entropy[\$itemp]=\$entropy[\$itemp]-\$kboltzmann*log(1-exp(-\$homega/\$kbt)); \$entropy[\$itemp]=\$entropy[\$itemp]+\$kboltzmann*\$homega/\$kbt*(1/(exp(\$homega/\$kbt)-1)); } } } last if eof; }
\$ezp=\$ezp; #/ \$zahler;
for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { my \$temp=(\$itemp-1)*\$tstep+\$tmin; \$entropy[\$itemp]=\$entropy[\$itemp]; # / \$zahler; \$fhelmholtz=\$uenergy+\$ezp-\$temp*\$entropy[\$itemp]; printf ("%15.8e %15.8e\n",\$temp,\$fhelmholtz); } printf ("#NOTEMP %15.8e\n",\$uenergy); EOF
if [ ! -d "$DIR" ];then echo Please specify the source directory. exit fi
forNAMEin 3RD.POSCAR.* do NUM="${NAME#3RD.POSCAR.}" # the number of the 3nd POSCAR file CODE=$(head -1 $NAME) # the code of the POSCAR file # find the source POSCAR file in the source directory INFO=$(grep $CODE $DIR/3RD.POSCAR.* | head -1) [ x"$INFO" == x ] && echo "NOT FOUND for job-$NUM" && continue SOURCE_NUM="${INFO%:*}" # the number of the source POSCAR file SOURCE_NUM="${SOURCE_NUM##*.}"
#module load spglib/master #thirdorder_vasp.py sow 331 -10 #生成2*2*2的超胞 0.6表示截断半径为6埃。 #sow a b c m abc表示生成的超胞大小，m为正数的时候表示截断半径（单位：nm）（例如0.6，表示截断半径为6.0埃）；为负数时表示第几近邻截断（此时应取负整数，例如-3表示第三近邻）
for i in 3RD.POSCAR.* do s=$(echo $i| cut -d "." -f 3) d=job-$s mkdir $d cp $i $d/POSCAR cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d pushd $d popd done
批量提交任务
{.line-numbers}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#!/bin/sh root_path='pwd'
#在每个job文件夹里进行静态计算。 for i in {001..1050} #按顺序从job-001到1050做静态vasp计算,根据生成的文件夹设定 do cd job-$i echo "job-$i", time sbatch sub.sh wait cd .. done
Besides the routine output files from previous ShengBTE program, FourPhonon generates these output files:
BTE.Numprocess_4ph: number of allowed four-phonon scattering processes, for each irreducible q point and phonon band
BTE.P4: phase space available for four-phonon processes, for each irreducible q point and phonon band
BTE.P4_total: total volume in phase space available for four-phonon processes
BTE.P4_plusplus*, BTE.P4_plusminus*, BTE.P4_minusminus*: similar to BTE.P4 but only includes contributions from ++/+-/- - processes
*Note for four-phonon scatterings, there are three different channels: recombination (++), redistribution (+-) and splitting (- -) processes.
Under temperature-dependent directories: (the unit of scattering rates is ps−1 if not specified)
BTE.WP4: weighted phase space available for four-phonon processes ([rad/ps]$^{-5}$, 2nd column) vs angular frequency (rad/ps, 1st column) for each irreducible q point and phonon band
BTE.WP4_plusplus*, BTE.WP4_plusminus*, BTE.WP4_minusminus*: similar to BTE.WP4 but only includes contributions from ++/+-/- - processes
BTE.w_3ph: three-phonon scattering rates for each irreducible q point and phonon band, this file replaces the original output BTE.w_anharmonic. Absorption and emission processes are written out into BTE.w_3ph_plus and BTE.w_3ph_minus
BTE.w_4ph: four-phonon scattering rates for each irreducible q point and phonon band. Similarly, we provide the contributions from different channels: BTE.w_4ph_plusplus, BTE.w_4ph_plusminus and BTE.w_4ph_minusminus
BTE.w_4ph_normal: four-phonon scattering rates of normal processes, for each irreducible q point and phonon band. This file has 5 columns and each one represents: angular frequency in rad/ps, recombination channel, redistribution channel, splitting channel, overall scattering rates from normal processes
BTE.w_4ph_Umklapp: four-phonon scattering rates of Umklapp processes, for each irreducible q point and phonon band. The format is the same as BTE.w_4ph_normal
BTE.kappa*: thermal conductivity results are written out as usual
1987年，Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交，不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法（或者称为 density perturbation functional theory，DFPT），直接计算出原子的移动而导致的势场变化，再进一步构造出动力学矩阵。进而计算出声子谱。
# -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np from ase.ioimport read, write import os
# 获取虚频对应的OUTCAR部分 vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行 #print(vib_lines) vib_dis = [] for line invib_lines: #model_positions = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 vib_dis.append(vib_infor) vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
#!/usr/bin/env python # -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np import os
# 获取虚频对应的OUTCAR部分 num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()]) vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
model_positions = [] vib_dis = [] for line invib_lines: position = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 model_positions.append(position) vib_dis.append(vib_infor)
###保存结构 f_out = open('POSCAR_new','w') f_out.writelines(lines_pos[:8]) f_out.write('Cartesian\n') for i innew_positions: f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n') f_out.close()
随后进行数据提取，以vaspkit1.12版本为例，打开vaspkit，输入73，打开VASP2other Interface 子功能，选择739 Sort Phonon Band Structure for Phononpy，运行完成即可获得group velocity相关数据以及各原子的声子谱数据；在运行过程中我们可以看到软件说明find group velocity并写入输出文件中。 1.5版本为789功能。
generate 1000 supercells with displacements at 300K.
New in v2.25 Symfc (symfc/symfc) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use symfc, its python module has to be installed via conda-forge or pip.
New in v2.3 ALM (ttadano/ALM) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use ALM, its python module has to be installed via conda-forge or building it. The installation instruction is found here.
When ALM is used, please cite the paper: T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).
Dielectric constants, piezoelectric constants and polar-phonon frequency：
DFPT方法，Note, DFPT cannot be used with hybrid exchange-correlation functionals. In these cases the LCALCEPS flag should be used in combination with IBRION = 6
amset run –directory path/to/files 运行即可 其他设置可以通过命令行加入： amset run –interpolation-factor 20 画图： amset plot rates mesh_99x99x99.h5
Output filesThe transport results will be written to transport_99x99x99.json. The 99x99x99 is thefinal interpolated mesh upon which scattering rates and the transport properties arecomputed.The write_mesh option is set to True, meaning scattering rates for all doping levelsand temperatures will be written to the mesh_99x99x99.h5 file.
zero_weighted_kpoints： 设置方式： keep: Keep zero-weighted k-points in the band structure. drop: Drop zero-weighted k-points, keeping only the weighted k-points. prefer: Drop weighted-kpoints if zero-weighted k-points are present in the calculation (useful for cheap hybrid calculations).
Use amset plot band --stats vasprun.xml: This calculates the effective mass in the same way as sumo but using a band structure interpolated from a uniform calculation.
orbitalWise 这个参数一般加在cohpBetween或者cohpGenerator之后，用以把原子之间的相互作用具体分到s p d 这些轨道上，比如计算s-s, p-p, d-d之间的pCOOP, pCOHP； 比如：cohpBetween atom 1 atom 2 orbitalWise
cohpGenerator from 1.4 to 1.5 type Ga type Sr orbitalWise
在 /vaspkit.0.72/examples/LDOS_PDOS/Partial_DOS_of_CO_on_Ni_111_surface 目录下启动 vaspkit ，输入命令 115 选择子功能 The Sum of Projected DOS for Selected Atoms and orbitals 。我们需要提取的是O的s，p轨道，C的s，p轨道以及表面的dxy和d3z2-r2(dz2)轨道。首先提示你选择元素（累加），第一次输入元素O，回车后提示你输入提取的轨道名（累加），第一次输入轨道s。回车后重复以上两个操作。如果想结束输入，在元素选择行直接按回车键结束输入。
元素行接受自由格式输入，你可以输入以下内容元素行接受自由格式输入，你可以输入以下内容1-3 4 Ni表示选择元素1,2,3,4和Ni元素的PDOS进行累加。轨道行只支持标准输入，如果使用了LORBIT=10不投影轨道，你只能选择 s p d f中的一个或多个，如果使用了LORBIT=11投影了轨道，你可以从s p d f和s py pz px dxy dyz dz2 dxz dx2 f-3 f-2 f-1 f0 f1 f2 f3中选择一个或多个输入。轨道行还支持输入all计算所有轨道的DOS之和。
比如元素行输入C O ,对应的轨道选择s px，那么提取的就是所有C和O元素的s轨道和px轨道的PDOS之和。 PDOS_USER.dat中的轨道名 为#Energy C&O_s&px，显而易懂。
以 CO 分子吸附在 Ni(100) 表面为例： 步骤一：优化 CO/Ni(100) 的结构； 步骤二：分别计算 Ni(100) 和 CO 的单点能 CO 和 Ni(100) 片段的坐标从 CO/Ni(100) 的 CONTCAR 里直接截取，不要再结构优化！！计算时也要保证三次自洽计算所采用的 FFT mesh 一致（NGXF,NGYF,NGZF） 步骤三：用 vaspkit 314 功能做电荷差分。依次输入三个片段的 CHGCAR 路径:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
314 ======== File Options ======== Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./co/CHGCAR ./slab/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./co/CHGCAR File... -->> (04) Reading Charge Density From ./co/CHGCAR File... -->> (05) Reading Structural Parameters from ./slab/CHGCAR File... -->> (06) Reading Charge Density From ./slab/CHGCAR File... -->> (07) Written CHGDIFF.vasp File!
生成 CHGDIFF.vasp 包含电荷密度差的信息，可以直接导入到 VESTA 里作图。 图片
Vaspkit计算变形电荷密度差
变形电荷密度差：自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)。 以 CO 分子为例： 步骤一：先自洽计算优化 CO 分子。 步骤二：新键文件夹，在优化好的结构基础上用 ICHARG = 12 做非自洽计算。 步骤三：用 vaspkit 314 功能做电荷差分。依次输入两个 CHGCAR 路径。
1 2 3 4 5 6 7 8 9 10 11 12 13 14
314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./deform/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./deform/CHGCAR File... -->> (04) Reading Charge Density From ./deform/CHGCAR File... -->> (05) Written CHGDIFF.vasp File!
314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./noele/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./noele/CHGCAR File... -->> (04) Reading Charge Density From ./noele/CHGCAR File... -->> (05) Written CHGDIFF.vasp File! 步骤四：得到 CHGDIFF.vasp 文件，导入到 VESTA 里做图。
ALGO = Damped HFSCREEN = 0.2 LHFCALC = .True. LMAXFOCK = 4 PRECFOCK = Fast TIME = 0.4
+SOC计算电子结构：
结构优化（PBE）
静态自洽scf计算（PBE+SOC）KPOINTS和POTCAR文件不变，写波函数和电荷密度
非自洽计算（PBE+SOC） Spin-Orbit Coupling Calculation
{.line-numbers}
1 2 3 4 5 6 7 8 9
ISPIN = 2 NELMIN = 6 (Min electronic SCF steps) LSORBIT = .TRUE. (ActivateSOC) GGA_COMPAT = .FALSE. (Apply spherical cutoff on gradient field) #VOSKOWN = 1 (Enhances the magnetic moments and the magnetic energies) LMAXMIX = 2 (For d elements increase LMAXMIX to 4, f: LMAXMIX = 6) SAXIS = 001 (Directionof the magnetic field) #MAGMOM = 003 (Setthis parameters manually, Local magnetic moment parallel to SAXIS, 3*NIONS*1.0for non-collinear magnetic systems) NBANDS = 32 (Setthis parameters manually, 2 * number of bands of collinear-run)
调节SOC强度
调节自旋轨道耦合作用的大小
要调优VASP中的自旋轨道耦合强度，需要修改源文件vasp_source_code_path/src/relativistic.F 并重新编译。如果要将强度减少一半，可以在relativistic.F的129行L●S项乘以0.5d0，如下： DO I=0,1DO J=0,1DO M =1,2LL+1DO MP=1,2LL+1 DLLMM(LMP+MP-1,LM+M-1,J+2I+1)=DLLMM(LMP+MP-1,LM+M-1,J+2I+1)+ & 0.5d0SUMLS(M,MP,I+2*J+1,LL) !!!line 129 relativistic.F fileEND DOEND DOEND DOEND DO
# Frequency dependent dielectric tensor including # local field effects within the RPA (default) or # including changes in the DFT xc-potential (LRPA=.FALSE.). # N.B.: beware one first has to have done a # calculation with ALGO=Exact, LOPTICS=.TRUE. # and a reasonable number of virtual states (see above) ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50 # be sure to take the same number of bands as for # the LOPTICS=.TRUE. calculation, otherwise the # WAVEDER file is not read correctly NBANDS = 64 # Add this to update the quasiparticle energies # in the Green's function (GW0) #NELM = 4 ISMEAR = 0 SIGMA = 0.05 EDIFF = 1E-8
k-point 1 : 0.0000 0.0000 0.0000 band No. KS-energies QP-energies sigma(KS) V_xc(KS) V^pw_x(r,r') Z occupation .........
LWANNIER90计算能带结构
重新开始第三步
INCAR：
1 2 3 4 5 6 7 8
ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50 #LRPA = .FALSE. ## be sure to take the same number of bands as for ## the LOPTICS=.TRUE. calculation, otherwise the ## WAVEDER file is not read correctly NBANDS = 64 ##VASP2WANNIER90 LWANNIER90=.TRUE.
This run generates the wannier90 standard output (wannier90.wout) and the file wannier90.chk needed for the wannier interpolation (next step)
Obtain bandstructure (Wannier interpolation) Uncomment the bandstructure plot flags in wannier90.win and rerun (restart) wannier90:
wannier90.x wannier90
This run generates the following bandstructure files which can be visualized using xmgrace or gnuplot:
wannier90_band.agr
wannier90_band.dat
wannier90_band.gnu
to plot the band structure using gnuplot: gnuplot -persist wannier90_band.gnu
Beyond the random-phase-approximation
To include local field effects beyond the random-phase-approximation in the description of the frequency dependent dielectric response function (local field effects in DFT) add the following line to your INCAR file:
LRPA = .FALSE. and again restart from the WAVECAR and WAVEDER files from step 2.
Beyond G0W0: GW0
The most usual step beyond single-shot GW (G0W0) is to iterate the quasi-particle energies in the Greens functions. This is the socalled GW0 approximation. To have VASP do, for instance, 4 iterations of the QP-energies in G, add the following line to the INCAR file:
NELM = 4 and again restart from the WAVECAR and WAVEDER files from step 2.
To quickly find the QP-energy of the highest lying occupied state after 4 iterations of the QP energies in G, type:
，属于四方晶系。
