手动计算高对称路径的方法两则

0：获得高对称点

详见文章“QE计算声子谱”第二节
经验证，以下两种方式的高对称路径结果完全相同

1：直接写脚本

import numpy as np
#每两个高对称点之间插入多少个点
number = 30
#scf.out里面的倒格失
vec_x=np.array([0.000000,  0.322290,  0.322290])
vec_y=np.array([0.322290,  0.000000,  0.322290])
vec_z=np.array([0.322290,  0.322290,  0.000000])
#高对称点的个数
number_of_path = 6
#高对称点的坐标
x0_vector=np.array([0, 0.5,  0,0.25,0,  0])
y0_vector=np.array([0,-0.5,  0,0.25,0,  0])
z0_vector=np.array([0, 0.5,0.5,0.25,0,0.5])

x_vector=[]
y_vector=[]
z_vector=[]
for n in range(number_of_path):
        x = x0_vector[n]*vec_x[0]+y0_vector[n]*vec_x[1]+z0_vector[n]*vec_x[2]
        y = x0_vector[n]*vec_y[0]+y0_vector[n]*vec_y[1]+z0_vector[n]*vec_y[2]
        z = x0_vector[n]*vec_z[0]+y0_vector[n]*vec_z[1]+z0_vector[n]*vec_z[2]
        x_vector.append(x)
        y_vector.append(y)
        z_vector.append(z)

print(x_vector)
print(y_vector)
print(z_vector)

f = open("path.dat", "w")
f.writelines(["{} {} {}\n".format(x_vector[0],y_vector[0],z_vector[0])])
for n in range(number_of_path-1):
        for i in range(number):
                x1 = (i+1)*(x_vector[n+1]-x_vector[n])/number+x_vector[n]
                x2 = (i+1)*(y_vector[n+1]-y_vector[n])/number+y_vector[n]
                x3 = (i+1)*(z_vector[n+1]-z_vector[n])/number+z_vector[n]
                f.writelines(["{} {} {}\n".format(x1,x2,x3)])
            
f.close()

2：从band.out里面读取

2.1进行一步自洽计算

&CONTROL
 calculation = 'scf'
  pseudo_dir = './'
     tstress = .true. ,
     tprnfor = .true. ,
/
&SYSTEM
       ibrav = 0
   celldm(1) = 1.88972688
         nat = 4
        ntyp = 2
     ecutwfc = 60.0
     ecutrho = 600
 occupations = 'smearing',
    smearing = 'mp',
     degauss = 0.01

/
&ELECTRONS
 mixing_beta = 0.3
 conv_thr = 3.000000e-10
/
ATOMIC_SPECIES
H  1.00784  h_pbe_v1.4.uspp.F.UPF
S  32.0590  s_pbe_v1.4.uspp.F.UPF

CELL_PARAMETERS
  -1.551398696   1.551398696   1.551398696
   1.551398696  -1.551398696   1.551398696
   1.551398696   1.551398696  -1.551398696

ATOMIC_POSITIONS (crystal)
H             0.0000000000        0.5000000000        0.5000000000
H             0.5000000000       -0.0000000000        0.5000000000
H             0.5000000000        0.5000000000        0.0000000000
S             0.0000000000        0.0000000000       -0.0000000000

K_POINTS automatic
32 32 32 0 0 0
                   

2.2进行一步band计算

&control
    calculation = 'bands',
    pseudo_dir = '.',
    outdir = '.',
    wf_collect = .true.
    verbosity = 'high'
    tstress = .t.,
    tprnfor = .t.,
 /
 &system
    ibrav = 0, celldm(1) =1.88972688,
    nat = 4, ntyp = 2,
    ecutwfc = 60, ecutrho = 600
    occupations='smearing', smearing='mp', degauss=0.02,
    nbnd = 90
 /
 &electrons
    mixing_beta = 0.7
    conv_thr = 1.0d-8
 /
ATOMIC_SPECIES
H  1.00784  h_pbe_v1.4.uspp.F.UPF
S  32.0590  s_pbe_v1.4.uspp.F.UPF

CELL_PARAMETERS
  -1.551398696   1.551398696   1.551398696
   1.551398696  -1.551398696   1.551398696
   1.551398696   1.551398696  -1.551398696

ATOMIC_POSITIONS (crystal)
H             0.0000000000        0.5000000000        0.5000000000
H             0.5000000000       -0.0000000000        0.5000000000
H             0.5000000000        0.5000000000        0.0000000000
S             0.0000000000        0.0000000000       -0.0000000000
K_POINTS crystal_b
6
0      0     0     30
0.5   -0.5   0.5   30
0      0     0.5   30
0.25   0.25  0.25  30
0      0     0     30
0      0     0.5   30

2.3用脚本从band.out里读取

#!/bin/bash

inputfile=band.out
numk=$(grep 'number of k points=' $inputfile | awk -F '=' '{print $2}' | awk '{print $1}')
# cartesian
echo $numk cartesian > path_cart.dat
grep -A $numk 'cart. coord. in units 2pi/alat' $inputfile | grep = | awk -F '(' '{print $3}' | awk -F ')' '{print $1}' | awk '{printf "%12.7f%12.7f%12.7f%6.1f\n",$1,$2,$3,1.0}' >> path_cart.dat
# crystal
echo $numk crystal > path_cryst.dat
grep -A $numk '                cryst. coord.' $inputfile | grep = | awk -F '(' '{print $3}' | awk -F ')' '{print $1}' | awk '{printf "%12.7f%12.7f%12.7f%12.7f\n",$1,$2,$3,1.0/'$numk'}' >> path_cryst.dat