Skip to content

PWscf Operation Procedure

by Ziyue Lin

I. Kinetic Energy Cutoff (Encut) Test

(1) Prepare Pseudopotential File

The required file format is *.pbe-n-rrkjus_psl.1.0.0.UPF. The specific filename will vary.

You can download pseudopotentials from the official Quantum Espresso website (URL: http://www.quantum-espresso.org/). Navigate to PSEUDOPOTENTIALS, click on the desired element in the periodic table, find the file matching the *.pbe-n-rrkjus_psl.1.0.0.UPF format, right-click, and select "Save Link As..." to download it.

  • NCPP: Norm-Conserving Pseudopotentials
  • USPP: Ultrasoft pseudopotentials
  • PAW: Projector Augmented-Wave method

(2) Prepare the script file run.pbs

(Note: In the original document, comments were marked in blue, and parameters to be modified were in red. This has been translated into English comments.)

Important: The script may fail if it is in DOS/Windows text format, showing an error like qsub: script is written in DOS/Windows text format.

To fix this, convert the file format to Unix: - Open the file with vi/vim: vi xxx.pbs - In command mode, type: :set ff=unix - Alternatively: :set fileformat=unix - Save and quit with :wq

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#!/bin/bash
#
#PBS -q CT6
#PBS -N encut # Job name
#PBS -l nodes=1:ppn=36 # Modify this line and -q according to the number of cores of your machine
#PBS -j oe
#PBS -V
cd $PBS_O_WORKDIR
# The following two lines may not be necessary for all machines
# ulimit -s 40000 
# ulimit -n 4096
# This needs to be set to the specific installation path of the Quantum Espresso executables.
# You can find this path by checking the submission scripts in a colleague's working directory.
export PW_ROOT=/home/ddf/program/espresso-5.4.0/bin 

mkdir -p ./tmp
rm -rf ./tmp/*

for ec in 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120
do
cat > F.scf.in_$ec << EOF
&control
calculation='scf',
restart_mode='from_scratch',
prefix='F' # Can be left unchanged; if modified, replace all instances of 'F' in the script.
pseudo_dir = '.',
outdir='./tmp'
tstress=.t.,
tprnfor=.t.
/
&system
ibrav = 0,
celldm(1)= 1.88972688,
nat= 8, # Number of atoms in the structure (count the lines under ATOMIC_POSITIONS)
ntyp= 1, # Number of atomic species
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.02
ecutwfc=$ec,
/
&electrons
mixing_beta = 0.7
conv_thr = 1.0d-8
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF # Element, Atomic Mass, Pseudopotential file
# If there are other elements, add them on new lines.
CELL_PARAMETERS # Structure with symmetry applied
2.605373874519238 0.000000000000000 0.000000000000000
0.000000000000000 2.605373874519239 0.000000000000000
0.000000000000000 0.000000000000000 2.605373874519239
ATOMIC_POSITIONS {crystal}
F 0.0000000000000000 0.5000000000000000 0.7500000000000000
F 0.0000000000000000 0.5000000000000000 -0.7500000000000000
F 0.7500000000000000 0.0000000000000000 0.5000000000000001
F -0.7500000000000000 0.0000000000000000 0.5000000000000000
F 0.5000000000000000 0.7500000000000000 0.0000000000000000
F 0.5000000000000000 -0.7500000000000001 0.0000000000000000
F 0.0000000000000000 0.0000000000000000 0.0000000000000000
F 0.5000000000000000 0.4999999999999999 0.5000000000000000
# The above two sections are obtained from the cell file's lattice parameters.
K_POINTS {automatic}
12 12 12 0 0 0
EOF
mpirun -n 36 $PW_ROOT/pw.x <F.scf.in_$ec >F.scf.out_$ec
done
rm -rf ./tmp*

(3) Prepare the getE.sh file

1
2
3
4
5
6
7
8
9
#!/bin/sh
echo "happy everyday....."
for a in 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120
do
pre=$(awk '/kinetic-energy cutoff/{ print $4 }' F.scf.out_$a |tail -n 1)
ent=$(awk '/! total energy/{print $5}' F.scf.out_$a |tail -n 1)
echo $pre $ent >> pe.txt
done
echo "end"

(4) Submit the job: qsub run.pbs

(5) If the output files show no errors, proceed to the next step.

(6) Make the script executable and run it

First, grant executable permissions: chmod a+x getE.sh. Then, run the script: ./getE.sh.

This will generate a pe.txt file with content like this: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
40.0000 -376.71234415
45.0000 -376.72445036
50.0000 -376.72957826
55.0000 -376.73705113
60.0000 -376.74526493
65.0000 -376.75168878
70.0000 -376.75550688
75.0000 -376.75729302
80.0000 -376.75786603
85.0000 -376.75799146
90.0000 -376.75812468
95.0000 -376.75837947
100.0000 -376.75875647
105.0000 -376.75914959
110.0000 -376.75949136
115.0000 -376.75973600
120.0000 -376.75987688
The cutoff energy can be chosen when the energy difference between consecutive steps is less than approximately 0.001 Ry. Here, we choose encut = 80 Ry.

II. Smearing Width (degauss) Test

(1) Prepare scf.in file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
&control
calculation='scf'
restart_mode='from_scratch',
prefix='F', # Same as before
pseudo_dir = '.',
outdir='./tmp'
/
&system
ibrav = 0,
celldm(1)=1.88972688,
nat= 8, # Same as before
ntyp= 1, # Same as before
ecutwfc=80, # Value obtained from the Encut test
ecutrho=800, # Typically 4x ecutwfc for NCPP, 10x for Ultrasoft PP
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.01
la2F = .true.,
/
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS # The structure must be symmetrized beforehand to avoid errors
2.60537300000000 0.00000000000000 0.00000000000000
0.00000000000000 2.60537300000000 0.00000000000000
0.00000000000000 0.00000000000000 2.60537300000000
ATOMIC_POSITIONS {crystal}
F 0.000000000117762 0.499999999338523 0.749986558395997
F 0.000000000339613 0.500000000256769 -0.749986558605385
F 0.749987711283413 0.000000000252767 0.499999999798142
F -0.749987712380392 0.000000000217478 0.500000000505681
F 0.500000000854287 0.749988023424662 0.000000002496860
F 0.500000000234474 -0.749988023462086 -0.000000002382903
F -0.000000000905482 -0.000000000090146 -0.000000000183659
F 0.500000000456325 0.500000000062032 0.499999999975266
K_POINTS {automatic} # A very dense mesh can be slow and hard to clean up; consider machine performance
10 10 10 0 0 0
EOF

(2) Prepare ph.in file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
Electron-phonon coefficients for F
&inputph
prefix='F',
fildvscf='Fdv',
amass(1)= 18.9984, # Mass of the F atom; if other atoms exist, add amass(2) = xxxx, etc.
outdir='./tmp',
fildyn='F.dyn',
electron_phonon='interpolated'
el_ph_sigma=0.005,
el_ph_nsigma=10,
trans=.true.,
ldisp=.false.
tr2_ph = 1.0d-12 /
0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00

(3) Prepare Pseudopotential File (same as for the encut test)

(4) Prepare getdegauss_1.sh file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#!/bin/bash
#by authors mayanbin 20150924
date=`date +%Y%m%d`
echo -e "\033[032m$date\033[0m"
awk '/Gaussian Broadening:/{print $3}' elph.* > degauss.dat
sed -e "s/( /a/g" elph.* | awk '/lambda/{print $2}' > lambda.dat
wait
# The first number is 3 times the number of atoms (e.g., 3 * 8 = 24)
awk '{print $1}' lambda.dat | head -n 24 > 0.005.dat 
# The ranges below should cover all groups of 24 lines
awk '{print $1}' lambda.dat | sed -n '25,48p' > 0.010.dat 
awk '{print $1}' lambda.dat | sed -n '49,72p' > 0.015.dat
awk '{print $1}' lambda.dat | sed -n '73,96p' > 0.020.dat
awk '{print $1}' lambda.dat | sed -n '97,120p' > 0.025.dat
awk '{print $1}' lambda.dat | sed -n '121,144p' > 0.030.dat
awk '{print $1}' lambda.dat | sed -n '145,168p' > 0.035.dat
awk '{print $1}' lambda.dat | sed -n '169,192p' > 0.040.dat
awk '{print $1}' lambda.dat | sed -n '193,216p' > 0.045.dat
awk '{print $1}' lambda.dat | tail -n 24 > 0.050.dat
wait
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.005.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.010.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.015.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.020.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.025.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.030.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.035.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.040.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.045.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.050.dat >> log
echo -e "\033[034mend\033[0m"
cp log ../log_$i
cp degauss.dat ../
rm log* *.dat job.pbs*
cd ../
done
paste degauss.dat log_0.005 log_0.01 log_0.015 log_0.02 log_0.025 log_0.03 log_0.035 log_0.04 log_0.045 log_0.05 >> pe.txt
rm log_* degauss.dat

(5) Prepare rmdat.sh file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
#!/bin/bash
#by mayanbin
rm pe.txt
for i in 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
do
cd $i
rm log*
rm *.dat
rm -r tmp
cd ../
done

(6) Prepare getdate_2.sh file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#!/bin/bash
#by mayanbin
for i in 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
do
cp getdegauss_1.sh getlambda_2.sh $i
cd $i
awk '/Gaussian Broadening:/{print $3}' elph.* > degauss.dat
sed -e "s/( /a/g" elph.* | awk '/lambda/{print $3}' > lambda.dat
wait
# The first number is 3 times the number of atoms (e.g., 3 * 8 = 24)
awk '{print $1}' lambda.dat | head -n 24 > 0.005.dat 
# The ranges below should cover all groups of 24 lines, same as before
awk '{print $1}' lambda.dat | sed -n '25,48p' > 0.010.dat
awk '{print $1}' lambda.dat | sed -n '49,72p' > 0.015.dat
awk '{print $1}' lambda.dat | sed -n '73,96p' > 0.020.dat
awk '{print $1}' lambda.dat | sed -n '97,120p' > 0.025.dat
awk '{print $1}' lambda.dat | sed -n '121,144p' > 0.030.dat
awk '{print $1}' lambda.dat | sed -n '145,168p' > 0.035.dat
awk '{print $1}' lambda.dat | sed -n '169,192p' > 0.040.dat
awk '{print $1}' lambda.dat | sed -n '193,216p' > 0.045.dat
awk '{print $1}' lambda.dat | tail -n 24 > 0.050.dat
wait
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.005.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.010.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.015.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.020.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.025.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.030.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.035.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.040.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.045.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.050.dat >> log
cp log ../log_$i
cp degauss.dat ../
rm log* *.dat job.pbs*
cd ../
done
paste degauss.dat log_0.005 log_0.01 log_0.015 log_0.02 log_0.025 log_0.03 log_0.035 log_0.04 log_0.045 log_0.05 >> pe.txt
rm log_* degauss.dat

(7) Prepare getlambda_2.sh file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#!/bin/bash
#by mayanbin 20150924
awk '{print $1}' lambda.dat | head -n 24 > 0.005.dat
awk '{print $1}' lambda.dat | sed -n '25,48p' > 0.010.dat
awk '{print $1}' lambda.dat | sed -n '49,72p' > 0.015.dat
awk '{print $1}' lambda.dat | sed -n '73,96p' > 0.020.dat
awk '{print $1}' lambda.dat | sed -n '97,120p' > 0.025.dat
awk '{print $1}' lambda.dat | sed -n '121,144p' > 0.030.dat
awk '{print $1}' lambda.dat | sed -n '145,168p' > 0.035.dat
awk '{print $1}' lambda.dat | sed -n '169,192p' > 0.040.dat
awk '{print $1}' lambda.dat | sed -n '193,216p' > 0.045.dat
awk '{print $1}' lambda.dat | tail -n 24 > 0.050.dat
wait
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.005.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.010.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.015.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.020.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.025.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.030.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.035.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.040.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.045.dat >> log
awk 'BEGIN{sum=0}{sum+=$1}END{print sum}' 0.050.dat >> log

(8) Prepare submission script job.pbs

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#!/bin/bash
#
#PBS -q CT6 # To be modified
#PBS -N Delph
#PBS -l nodes=1:ppn=36 # Modify nodes, cores per node, and path as needed
#PBS -j oe
#PBS -V
cd $PBS_O_WORKDIR
export PW_ROOT=/share/apps/compiler/QE/espresso-5.4.0/bin/
mkdir -p ./tmp
rm -rf ./tmp/*
# running program
# SCF at dense k-mesh, good enough for electronic DOS
mpirun -n 36 $PW_ROOT/pw.x <scf.in >scf.out
mpirun -n 36 $PW_ROOT/ph.x <ph.in >ph.out
# clean
rm -rf ./tmp/*

(9) Prepare cpelph.sh

(This step copies the lambda files so that getdate_2.sh can run.)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
cp 0.005/elph_dir/elph.inp_lambda.1 0.005/
cp 0.01/elph_dir/elph.inp_lambda.1 0.01/
cp 0.015/elph_dir/elph.inp_lambda.1 0.015/
cp 0.02/elph_dir/elph.inp_lambda.1 0.02/
cp 0.025/elph_dir/elph.inp_lambda.1 0.025/
cp 0.03/elph_dir/elph.inp_lambda.1 0.03/
cp 0.035/elph_dir/elph.inp_lambda.1 0.035/
cp 0.04/elph_dir/elph.inp_lambda.1 0.04/
cp 0.045/elph_dir/elph.inp_lambda.1 0.045/
cp 0.005/elph_dir/elph.inp_lambda.1 0.05/

(10) Prepare degauss_1.sh

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
#!/bin/bash
#by authors mayanbin
for d in 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05
do
mkdir $d
cp *.UPF ph.in job.pbs $d
cd $d
mkdir tmp
cd ../
cat ./scf.in | sed -e "s/degauss=0.02/degauss=$d/g" \
> $d/scf.in
cd $d
qsub job.pbs
cd ../
done

(11) Submit the script: qsub job.pbs (or bsub on some clusters)

This will produce *.dyn, Delph.o*, ph.out, scf.out, a tmp folder, and an elph_dir folder. Check the output files for errors before proceeding.

(12) Submission Sequence

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Grant permissions and submit to automatically generate a series of jobs
chmod +x degauss_1.sh
./degauss_1.sh

chmod +x cpelph.sh
./cpelph.sh

chmod +x getdate_2.sh
./getdate_2.sh
# This creates pe.txt. Import it into Origin to plot and determine the smearing width.

As shown in the figure, the smearing width is chosen as 0.03 (i.e., the position where several lines converge).

III. Structure Optimization

(1) Prepare Pseudopotential File (same as before)

(2) Apply Symmetry using phonopy

i) Manually convert the lattice parameters from CELL format to POSCAR format.

CELL format: 

1
2
3
4
5
6
7
CELL_PARAMETERS (alat= 1.88972688)
2.6086092816666668 0.0000000000000000 0.0000000000000000
0.0000000000000000 2.6086092816666668 0.0000000000000000
0.0000000000000000 0.0000000000000000 2.6086092816666668
ATOMIC_POSITIONS {crystal}
F 0.0000000000000000 0.5000000000000000 0.7500000000000000
...

POSCAR format: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
System Name
1.0
2.605373874519238 0.000000000000000 0.000000000000000
0.000000000000000 2.605373874519239 0.000000000000000
0.000000000000000 0.000000000000000 2.605373874519239
F
8
Direct
0.0000000000000000 0.5000000000000000 0.7500000000000000
...

ii) Save the structure as a file named POSCAR.

iii) Execute: phonopy --symmetry POSCAR

iv) This will generate PPOSCAR and BPOSCAR files.

v) Copy the parameters from the generated PPOSCAR file back into the script below, in the CELL file format.

(3) Prepare the script pwelph.pbs

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#!/bin/bash
#
#PBS -q CT6
#PBS -N pw-opt
#PBS -l nodes=1:ppn=36
#PBS -j oe
#PBS -V
cd $PBS_O_WORKDIR
export PW_ROOT=/share/apps/compiler/QE/espresso-5.4.0/bin/
# Add for specific servers if needed
#ulimit -s 40000 
#ulimit -n 4096
# running program
cat > F.scf.in << EOF
&control
calculation='vc-relax',
restart_mode='from_scratch',
dt=30
prefix='F.scf'
pseudo_dir = '.',
outdir='./tmp'
tstress=.t.,
tprnfor=.t.
nstep=10000
/
&system
ibrav = 0,
celldm(1)=1.88972688
nat= 8, # These four items are the same as in the smearing test scf.in
ntyp= 1,
ecutwfc=80,
ecutrho=800.
occupations='smearing',
smearing='mp',
degauss=0.03 # Result from the smearing width test
nosym = .f. # Set to .f. to apply symmetry operations
/
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
/
&IONS
ion_dynamics = 'damp'
upscale = 20
/
&CELL
cell_dynamics = 'damp-pr'
press = 30000 # Set pressure (unit: kbar, 1 kbar = 0.1 GPa)
/
ATOMIC_SPECIES # All below are the same as before
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS
2.60537300000000 0.00000000000000 0.00000000000000
0.00000000000000 2.60537300000000 0.00000000000000
0.00000000000000 0.00000000000000 2.60537300000000
ATOMIC_POSITIONS {crystal}
F 0.000000000117762 0.499999999338523 0.749986558395997
...
K_POINTS {automatic}
10 10 10 0 0 0
EOF
mpiexec -n 36 $PW_ROOT/pw.x <F.scf.in >F.scf.out

(4) Submit the script: qsub pwelph.pbs

(5) Format Conversion

The lattice parameters obtained in scf.out need to be converted back to POSCAR format. Use phonopy to apply symmetry again. (An SCF calculation might fail without applying symmetry, which can be solved by setting nosym = .f..) Then, input the final structure, in CELL format, into the subsequent SCF calculation.

IV. First Self-Consistent (SCF) Calculation

The purpose of this first SCF step is threefold: 1. Determine the number of symmetry operations to verify that the symmetry is correct. Find this under the Sym. Ops. field in scf.out. (Reference for point groups vs. number of symmetry operations: http://img.chem.ucl.ac.uk/sgp/large/sgp.htm) 2. Determine the number of q-points, which corresponds to the number of .dyn files. Find this under the cart. coord. field. 3. Check the total stress to see if it corresponds to the desired pressure. Find this under the total stress field.

Create a new directory, e.g., 3-scf

(1) Prepare Pseudopotential File

(2) Prepare scf.in (q-points must match ph.in)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
&control
calculation='scf'
restart_mode='from_scratch',
prefix='F',
pseudo_dir = './',
outdir='./tmp'
tstress=.t.,
verbosity='high'
/
&system
ibrav = 0,
celldm(1)=1.88972688,
nat= 8, # Number of atoms
ntyp= 1, # Number of atomic species
ecutwfc=80, # From cutoff energy test
ecutrho=800, # Typically 4x ecutwfc for NCPP, 10x for Ultrasoft PP
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.03 # From smearing width test
la2F = .true.,
/
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF # Atomic mass, pseudopotential, lattice parameters, etc.
CELL_PARAMETERS (alat= 1.88972688)
2.6086092816666668 0.0000000000000000 0.0000000000000000
...
ATOMIC_POSITIONS {crystal}
F 0.0000000000000000 0.5000000000000000 0.7500000000000000
...
K_POINTS {automatic}
4 4 4 0 0 0 # Set this to the q-point grid specified in ph.in (e.g., nq1=4, nq2=4, nq3=4)

(3) Prepare Submission Script

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
#!/bin/bash
#PBS -q CT1
#PBS -l nodes=1:ppn=12
#PBS -j oe
#PBS -V
#PBS -N scf-3
cd $PBS_O_WORKDIR
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/pw.x <scf.in> scf.out
#clean
rm -rf tmp

(4) Submit the script

V. Second Self-Consistent (SCF) Calculation

Create another new directory, e.g., 4-scf

(1) Prepare Pseudopotential File

(2) Generate Specific q-points for Band Structure

(An example is provided in the Genq_qpoints folder.)

First, create a new directory Genq_qpoints.

i) Prepare necessary files.

ii) First, import the optimized and symmetrized POSCAR into a visualization software (like VESTA or MS) to identify the high-symmetry points of the Brillouin zone.

Input the following into a file named syml:

1
2
3
4
5
6
7
5 # Number of high-symmetry points
20 20 20 20 # Insert 20 points in each of the 4 segments between the 5 high-symmetry points
X 0.500 0.000 0.000 # High-symmetry path
R 0.500 0.500 0.500
M 0.500 0.500 0.000
G 0.000 0.000 0.000
R 0.500 0.500 0.500

iii) Make the generator executable (chmod +x a.out) and run it (./a.out).

After execution, the specific k/q-points will be written to inp.kpt. Copy this list into the k-point section of the next scf.in file.

(3) Prepare scf.in (with explicit q/k-points)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
&control
calculation='scf'
restart_mode='from_scratch',
prefix='F',
pseudo_dir = './',
outdir='./tmp'
tstress=.t.,
verbosity='high'
/
&system
ibrav = 0,
celldm(1)=1.88972688,
nat= 8, 
ntyp= 1, 
ecutwfc=80, 
ecutrho=800, 
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.03 
la2F = .true.,
/
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (alat= 1.88972688)
2.6086092816666668 0.0000000000000000 0.0000000000000000
...
ATOMIC_POSITIONS {crystal}
F 0.0000000000000000 0.5000000000000000 0.7500000000000000
...
K_POINTS {crystal} # Note: using 'crystal' coordinates here
81 # The total number of k-points from inp.kpt
# Paste the values from the inp.kpt file here.
0.500000 0.000000 0.000000 1.00
0.500000 0.025000 0.025000 1.00
......
0.500000 0.050000 0.050000 1.00

(4) Prepare Submission Script

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
#!/bin/bash # Script is the same as the first SCF step
#PBS -q CT1
#PBS -l nodes=1:ppn=12
#PBS -j oe
#PBS -V
#PBS -N scf-4
cd $PBS_O_WORKDIR
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/pw.x <scf.in> scf.out
#clean
rm -rf tmp

(5) Submit the script

VI. Electron-Phonon Coupling and SCF Calculation

(1) Prepare Pseudopotential File (same as before)

(2) Prepare ph.in file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
Electron-phonon coefficients for F
&inputph
prefix = 'F'
fildvscf = 'dvscf'
fildyn = 'F.dyn'
amass(1)= 18.9984,
outdir='./tmp',
electron_phonon='interpolated',
el_ph_sigma=0.05,
el_ph_nsigma=10,
alpha_mix=0.5
recover=.true.
trans=.true.,
ldisp=.true.
nq1=4, # A coarse q-point grid of 4x4x4. A value between 3-5 is common.
nq2=4, # Larger values are more accurate but slower, and can help eliminate imaginary frequencies.
nq3=4,
tr2_ph = 1.0d-12
/

(3) Prepare scf.in file

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
&control
calculation='scf'
...
/
&system
...
ecutwfc=80,
ecutrho=800,
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.03
la2F = .true.,
/
&electrons
...
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS (alat= 1.88972688)
...
ATOMIC_POSITIONS {crystal}
...
K_POINTS {automatic}
8 8 8 0 0 0 # Set to 2x the q-point grid in ph.in (2 * 4 = 8)

(4) Prepare fit.scf.in file (same as scf.in, but k-point mesh is doubled again)

1
2
3
4
5
6
7
&control
...
/
&system
...
K_POINTS {automatic}
16 16 16 0 0 0 # Set to 4x the q-point grid in ph.in (4 * 4 = 16)

(5) Prepare Submission Script

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#!/bin/bash
#PBS -q CT1
#PBS -l nodes=1:ppn=12
#PBS -j oe
#PBS -V
#PBS -N QE-E.P.
cd $PBS_O_WORKDIR
# Modify paths and node settings as needed
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/pw.x <fit.scf.in> fit.scf.out
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/pw.x <scf.in> scf.out
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/ph.x <ph.in> ph.out
#clean
rm -rf tmp

Finally, submit the script.

VII. Phonon Spectrum and its Density of States

Now, for post-processing, prepare the following scripts:

(1) Create and run run_q2r_5.4.0

1
2
3
4
5
6
7
#!/bin/sh
cat > q2r.in << EOF
&input
fildyn='F.dyn', zasr='crystal', la2F = .true., flfrc='F.fc'
/
EOF
/share/apps/compiler/QE/espresso-5.4.0/bin/q2r.x < q2r.in > q2r.out

Then execute the script: 

1
2
chmod +x run_q2r_5.4.0
./run_q2r_5.4.0

(2) Create and submit run_dos-5.4.0

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
#!/bin/sh
#PBS -q CT1
#PBS -l nodes=1:ppn=12
#PBS -j oe
#PBS -V
#PBS -N phonon-dos
cd $PBS_O_WORKDIR
####################################################################
cat > matdyn.in.dos << EOF
&input
asr='simple',
amass(1)= 18.9984,
flfrc='F.fc',
flfrq='F.freq',
dos=.true.
fldos='phdos.dat',
nk1=30, nk2=30, nk3=30,
ndos=5000
/
EOF
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/PHonon/PH/matdyn.x < matdyn.in.dos > matdyn.out.dos

Then submit the script: qsub run_dos_5.4.0

(3) Create and run run_frequency

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
#!/bin/sh
#PBS -q CT1
#PBS -l nodes=1:ppn=12
#PBS -j oe
#PBS -V
#PBS -N phonon-frq
cd $PBS_O_WORKDIR
cat > matdyn.in.freq << EOF
&input
asr='simple', amass(1)= 18.9984,
flfrc='F.fc', flfrq='F.freq',la2F = .true., dos=.false.,
/
81
0.1916730 0.0000000 0.0000000 0.0246914
0.1916730 0.0095836 0.0095836 0.0246914
...
0.1916730 0.1916730 0.1916730 0.0246914
EOF
mpirun -np 12 /share/apps/compiler/QE/espresso-5.4.0/bin/matdyn.x < matdyn.in.freq > matdyn.out

Note on k-points: The list of k-points above should be taken from the output of the second SCF calculation (4-scf). Look for the section in scf.out that starts with number of k points= 81 and cart. coord. in units 2pi/alat. Be sure to use the Cartesian coordinates (cart. coord.), not the crystal coordinates. You may need to delete a column of data from what you copy.

Tip for column selection: - Linux (vim): Press Ctrl+v (visual block mode), use arrow keys to select the column, then press d to delete. - Windows: Paste the text into a text editor that supports block selection (like Notepad++ or VS Code) by holding Alt while selecting with the mouse.

Finally, execute the script: 

1
2
chmod +x run_frequency
./run_frequency

This will output the phonon spectrum and DOS data. - F.freq.gp is used for plotting the phonon spectrum. - phdos.dat is used for plotting the phonon density of states.

The resulting plot should look something like this:

VIII. Superconductivity Calculation

First, navigate into the elph_dir folder.

(1) Prepare the script run_lambda

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
#!/bin/sh
cat > lambda.in << EOF
80 0.52 1 # Line 1: [max_freq_THz] [smearing_width] [smearing_method=1]. The max frequency should be higher than the highest phonon frequency. A smearing of ~0.5 is typical.
11 # Line 2: Number of k-points. Get the Cartesian k-points and weights from the first SCF step (3-scf) output file, under "cart. coord.".
0.0000000 0.0000000 0.0000000 0.0160000
0.0000000 0.0000000 0.0766692 0.0960000
...
elph.inp_lambda.1
elph.inp_lambda.2
...
elph.inp_lambda.11
0.13 ! mu_star, the Coulomb pseudopotential. Typically between 0.1 - 0.13.
! Used in the modified Allen-Dynes formula for T_c (via omega_log)
EOF
/share/apps/compiler/QE/espresso-5.4.0/bin/lambda.x < lambda.in > lambda.out

(2) Make executable and run

1
2
chmod +x run_lambda
./run_lambda

(3) Result Analysis: The calculation produces lambda.out and lambda.dat files.

In the image above, the goal is to make the values outside and inside the parentheses (in the red box) as close as possible. Ideally, they should be identical up to the last few digits. This is mainly controlled by adjusting the parameters in lambda.in. The maximum frequency value has a large impact; setting it 1-2 THz higher than the actual maximum phonon frequency is a good starting point. The smearing width also affects convergence. To determine the superconducting transition temperature (Tc), it is recommended to plot Tc versus the Gaussian broadening (blue box) in Origin to check for convergence. The converged value is the calculated Tc for the structure.

IX. Calculating Only Phonon Spectrum and DOS

For this, you only need one SCF step. The scf.in and ph.in files are simplified by removing some parameters.

scf.in: (&system namelist)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
&system
ibrav = 0,
celldm(1)=1.88972688,
nat= 8,
ntyp= 1,
ecutwfc=80,
ecutrho=800,
occupations='smearing',
smearing='methfessel-paxton',
degauss=0.03
/

ph.in:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
Phonon calculation for F
&inputph
tr2_ph=1.0d-12,
prefix='F',
ldisp=.true.,
nq1=5, nq2=5, nq3=5
amass(1)= 18.9984,
outdir='./tmp',
fildyn='F.dyn',
alpha_mix=0.5
/

X. Superconductivity Calculation - Tetrahedron Method

All necessary files are assumed to be in the compressed archive Tc.tar.gz.

(1) Structure Optimization (Same as Section III)

Create a directory opt and perform structure optimization, which outputs an scf.out file.

(2) Self-Consistent Calculation

Prepare the files: *.scf.in, job.pbs, and pseudopotential files.

*.scf.in

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
&CONTROL
calculation = 'scf' ,
prefix='RbF',
pseudo_dir = './',
outdir='./tmp'
/
&SYSTEM
ibrav = 0,
celldm(1) = 1.88972688,
nat = 12,
ntyp = 2,
ecutwfc = 80 ,
ecutrho = 320 ,
occupations = 'tetrahedra_opt' ,
/
&ELECTRONS
/
ATOMIC_SPECIES
F 18.9984 F.pbe-n-kjpaw_psl.1.0.0.UPF
Rb 85.4678 Rb.pbe-spn-kjpaw_psl.1.0.0.UPF
CELL_PARAMETERS
...
ATOMIC_POSITIONS (crystal)
...
K_POINTS automatic
8 8 8 0 0 0
# Modify as in previous scf.in files

job.pbs

1
2
3
4
5
6
7
8
9
#!/bin/bash
#
#PBS -q CT3
#PBS -N RbF2_scf
#PBS -l nodes=1:ppn=28
#PBS -j oe
#PBS -V
cd $PBS_O_WORKDIR
mpirun -n 28 /share/apps/compiler/QE/espresso-6.4/bin/pw.x <RbF.scf.in >RbF.scf.out

Submit the script: qsub job.pbs. This will output a tmp folder.

(3) Phonon Calculation

Create a new directory phonon-2 for the second step. Copy the tmp folder from the SCF step into this directory.

Prepare the files: *.ph.in, job.pbs, and pseudopotential files.

*.ph.in

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
RbF phonon
&INPUTPH
prefix = 'RbF',
outdir = './tmp',
fildyn = 'RbF.dyn'
fildvscf = 'dv',
fildrho = 'drho',
ldisp = .true.,
lshift_q = .true.,
nq1 = 4,
nq2 = 4,
nq3 = 4,
/
# Set nq1, nq2, nq3 appropriately; 4 4 4 is common. The principle is similar to the previous method.

job.pbs

1
2
3
4
5
6
7
8
#!/bin/bash
#PBS -q CT3
#PBS -l nodes=1:ppn=28
#PBS -j oe
#PBS -V
#PBS -N RbF7_ph
cd $PBS_O_WORKDIR
mpirun -n 28 /share/apps/compiler/QE/espresso-6.4/bin/ph.x < RbF.ph.in >RbF.ph.out

Submit the script. This will generate a new tmp folder and a series of .dyn files.

(4) Electron-Phonon Calculation

Create a new directory elph-3. Copy the tmp folder and all .dyn files from phonon-2 into this directory.

Prepare the files: *.elph.in, job.pbs, pseudopotential files, and alpha2x.pbs.

*.elph.in

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
&INPUTPH
prefix = 'RbF',
outdir = './tmp',
fildyn = 'RbF.dyn'
fildvscf = 'dv',
fildrho = 'drho',
ldisp = .true.,
lshift_q = .true.,
nq1 = 4,
nq2 = 4,
nq3 = 4,
electron_phonon = "lambda_tetra"
nk1 = 4,
nk2 = 4,
nk3 = 4,
/
&INPUTa2F
nfreq = 800
/

job.pbs

1
2
3
4
5
6
7
8
#!/bin/bash
#PBS -q CT3
#PBS -l nodes=1:ppn=28
#PBS -j oe
#PBS -V
#PBS -N RbF6_elph
cd $PBS_O_WORKDIR
mpirun -np 28 /share/apps/compiler/QE/espresso-6.4/bin/ph.x < RbF.elph.in >RbF.elph.out

alpha2x.pbs

1
2
3
4
5
6
7
8
#!/bin/bash
#PBS -q CT3
#PBS -l nodes=1:ppn=28
#PBS -j oe
#PBS -V
#PBS -N RbF6_elph
cd $PBS_O_WORKDIR
mpirun -np 28 /share/apps/compiler/QE/espresso-6.4/bin/alpha2f.x < RbF.elph.in > alpha2f.out

First, submit job.pbs. After it completes and all *.dyn*.elph.* files have been generated, submit alpha2x.pbs.

Use gnuplot *.McMillan.gp to plot the results. Observe where the curve converges to find Tc. Other parameters are available in alpha2f.out.