Skip to content

Wien2k Pseudopotential Testing

I. Structure Optimization

First, create a new folder in the current directory, for example named vasp. In this folder, first perform multi-pressure point optimization. (Multi-pressure point optimization steps can be found in the VASP documentation. This document uses the vasp software for geometry optimization as preprocessing example)

II. Preparing Files

Three scripts: 0-Rename.sh, 1-sub.pbs, 2-job.sh

(All input files are in Wien2k_INPUT.rar)

chmod +x 0-Rename.sh, 2-job.sh (Grant permissions)

III. Script Contents

0-Rename.sh (depends on phonopy)

 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
#!/bin/sh
#Edit by Linziyue
#This series contains 0, 1, 2 three scripts (in order). Before starting, you need to calculate multi-pressure point optimization in the vasp folder first, then execute this script for copying and format conversion
#Note: after copying, pad the cif filename (pressure) to the same number of digits for easier script processing, for example, change 500 to 0500
#The loops in scripts 1, 2 follow the above rules
#The final output volume units have been converted to cubic angstroms

if [ -d CIF ];then
rm -r CIF
fi
mkdir CIF

for i in 0 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000
#Pressures to extract, unit is kbar
do
cp vasp/CONTCAR_$i CIF/
cd CIF/
phonopy -c CONTCAR_$i --symmetry --tolerance=0.01
rm BPOSCAR
mv PPOSCAR CONTCAR_$i
cabal poscar cif <CONTCAR_$i> CONTCAR_$i.cif
mv CONTCAR_$i.cif $i.cif
cd ..
echo "***********"
echo "Job Rename $i done"
echo "***********"
done

echo "*************"
echo "All Jobs Done"
echo "*************"

1-sub.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
#!/bin/bash
#Edit by Linziyue
#PBS -N wien2k_cal
#PBS -l nodes=1:ppn=12
#PBS -j n
#PBS -e ${PBS_JOBNAME}.err
#PBS -o ${PBS_JOBNAME}.out

cd $PBS_O_WORKDIR
NP=`cat $PBS_NODEFILE|wc -l`
cat $PBS_NODEFILE > nodes.info

for i in 0500 0750 1000 1250 1500 1750 2000 2250 2500 2750 3000
do
if [ -d "$i" ];then
rm -r $i
fi
mkdir $i
cp CIF/$i.cif $i;
cd $i ; cif2struct $i.cif

init_lapw1 << EOF
N
a
c
c
c
13
0.99
c
2000
0
c
n
EOF

sed -i "s/EF= 0.50000/EF= 2.0/" $i.in1
#Whether to add these two lines and how to modify the values is still unclear to me. Generally, you can execute this script normally (outputting broyd files in each folder normally) by adding or removing these two lines. Students who understand how to modify this parameter are welcome to discuss further with me.
sed -i "s/7.00/8.00/" $i.in1
run_lapw >1og 2>&1 &
cd ..
done
sleep 3m

2-job.sh

 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
#!/bin/sh
#Edit by Linziyue
echo "Let's Rock"

if [ -f job.txt ];then
rm job.txt
fi

if [ -d scf ];then
rm -r scf
fi
mkdir scf

for i in 0500 0750 1000 1250 1500 1750 2000 2250 2500 2750 3000
do
cp $i/$i.scf scf/
echo "***********"
echo "Job $i done"
echo "***********"
done

echo "*************"
echo "All Jobs Done"
echo "*************"

cd scf
grepline :ENE '*.scf' 1 > scf.analysis
grepline :VOL '*.scf' 1 >> scf.analysis
#Copy struct file from any pressure
cp ../3000/*.struct scf.struct

eplot_lapw << EOF
vol
0
B
EOF

IV. Submitting Scripts

./0-Rename.sh

qsub 1-sub.pbs

./2-job.sh

Output results are shown in the figure. After executing 0-Rename, a CIF folder is generated.

After submitting 1-sub.pbs, each pressure point will output a folder with many broyd files.

After executing 2-job.sh, an scf folder will be output containing the data needed for plotting.

The file name for plotting is: scf.eosfit. As shown in the right figure, the first column is volume and the third column is pressure.

Pressure unit is GPa, volume conversion to cubic angstroms requires multiplying by 0.148187.

V. Linear Fitting

B-M equation of state E-V curve fitting tutorial:

http://blog.sciencenet.cn/blog-3388193-1152915.html

Click Tools - Fitting Function Generator

For fitting the P-V equation of state, the expression is:

x=(3/2)B0((V0/y)(7.0/y)-(V0/y)(5.0/3))(1+(3.0/4)(B1-4)*((V0/y)^(2.0/3)-1))

Where x is pressure, y is volume; parameters are B0, B1, V0.

Note: The above link's corresponding blog post is currently blocked and cannot access the specific content of the B-M equation of state E-V curve fitting tutorial. If you need to view related tutorials, it is recommended to contact the author or search for other resources.

The above link is invalid, see the tutorial below. Just look at the fitting section. In the AI era, using Python for fitting is recommended

Birch-Murnaghan Equation of State Fitting WIEN2k E-V Curve Tutorial

1. Introduction

In density functional theory (DFT) based first-principles calculations, by calculating the total energy (E) for a series of unit cell volumes (V), an energy-volume (E-V) curve can be obtained. The Birch-Murnaghan equation of state is a widely used empirical formula for fitting this E-V curve to extract material equilibrium physical properties at absolute zero, mainly including:

  • Equilibrium volume \(V_0\)
  • Equilibrium energy \(E_0\)
  • Bulk modulus \(B_0\) (characterizes a material's ability to resist uniform compression)
  • First-order pressure derivative of bulk modulus \(B_0'\)

This tutorial will guide you through the entire process from WIEN2k calculation to final fitting.

2. Basic Principle: Birch-Murnaghan Equation of State

The most commonly used form is the third-order Birch-Murnaghan equation of state, expressed as:

\[ E(V) = E_0 + \frac{9V_0 B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^{\frac{2}{3}} - 1\right]^3 B_0' + \left[\left(\frac{V_0}{V}\right)^{\frac{2}{3}} - 1\right]^2 \left[6 - 4\left(\frac{V_0}{V}\right)^{\frac{2}{3}}\right] \right\} \]

Where:

  • \(E(V)\) is the total energy at volume \(V\).
  • \(E_0\) is the ground state equilibrium energy (fitting parameter).
  • \(V_0\) is the equilibrium volume (fitting parameter).
  • \(B_0\) is the equilibrium bulk modulus (fitting parameter, unit usually GPa).
  • \(B_0'\) is the first-order pressure derivative of \(B_0\) (fitting parameter, usually a dimensionless number between 3.5-5).

Simplified form (second-order): When \(B_0'\) is fixed at 4 (\(B_0''=0\)), the equation can be simplified to: $$ E(V) = E_0 + \frac{9V_0 B_0}{16} \left[\left(\frac{V_0}{V}\right)^{\frac{2}{3}} - 1\right]^2 \left{ 2 + \left[6 - 4\left(\frac{V_0}{V}\right)^{\frac{2}{3}}\right] \right} $$ The second-order form reduces one fitting parameter and is more stable when \(B_0'\) is close to 4 or when data points are few, but accuracy may be slightly lower. Generally, the third-order form is recommended.

3. WIEN2k Calculation Steps: Generating E-V Data Points

The core idea is: fix the unit cell shape and atomic positions, only uniformly scale the unit cell volume, perform self-consistent calculations at each volume, and extract total energy.

  1. Create initial structure: Use structgen or import from crystal database to create an initial structure file (.struct) close to equilibrium.
  2. Volume scaling:
  3. Determine a set of volume scaling factors (e.g., 0.94, 0.96, 0.98, 1.00, 1.02, 1.04, 1.06).

  4. For each scaling factor \(c\), use the following command to scale lattice constants:

    1
    init_lapw -numax 7 -rkmax 7 -v 1.00 # First initialize original structure

    Note: For each different volume, you need to perform calculations in a separate directory to avoid file conflicts.

    1
2
3
4
5
6
7
    # Assume your initial structure is in directory 'Fe_1.00'
mkdir ../Fe_$c
cp *.def Fe_$c/
cp *.struct Fe_$c/
cd ../Fe_$c
instgen_lapw // Use instgen_lapw to scale structure, input scaling factor c as prompted
# Or manually edit lattice constants (a, b, c) in .struct file
    3. Execute calculations: - In the directory corresponding to each volume, execute standard WIEN2k self-consistent calculation loop:

    1
2
    init_lapw -b -v $c # Set RKMAX and other parameters
run_lapw -ec 0.0001 -i 100 # Execute SCF calculation, set convergence criteria and maximum iterations
    - Ensure each calculation reaches energy convergence (check the :ENE change in the last line of case.scf). 4. Extract data: - After calculation completion, find the converged total energy (unit: Ry) from the last line of case.scf file in each directory. - Volume \(V\) can be calculated from lattice constants in .struct file (unit: \(a.u.^3\) or \(\AA^3\), ensure units are consistent). - Organize data into a two-column file (e.g., eV_data.dat):

    1
2
3
4
5
    # Volume(a.u.^3)    Energy(Ry)
150.0              -1000.12345
155.0              -1000.23456
160.0              -1000.34567
...

4. Fitting Methods

Method 1: Using Python (scipy.optimize.curve_fit)

This is a very flexible and powerful method.

 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
68
69
70
71
72
73
74
75
76
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 1. Define third-order B-M equation function
def birch_murnaghan_3rd_order(V, E0, V0, B0, B0p):
    """
    V: Volume
    E0: Equilibrium energy
    V0: Equilibrium volume
    B0: Bulk modulus (GPa)
    B0p: First-order derivative of B0
    Returns: Energy
    """
    # Convert B0 from GPa to Ry/(a.u.)^3 or other energy/volume units if necessary
    # Conversion factor depends on your unit system.
    # For example: 1 Ry = 13.6056980659 eV, 1 a.u.^3 = (0.529177)^3 A^3
    # If your energy is Ry and volume is (a.u.)^3, B0 unit is also Ry/(a.u.)^3, no conversion needed.
    # But typically B0 is reported in GPa in literature, so unit conversion is needed.
    # Unit conversion factor (1 Ry/Bohr^3 to GPa): 
    # 1 Ry/Bohr^3 = 14710.5076 GPa (approximate value)
    # Therefore, inside the function, we use B0_internal = B0 / 14710.5 to convert input GPa value to units needed for calculation
    B0_internal = B0 / 14710.5 # If input B0 is GPa, energy is Ry, volume is Bohr^3

    eta = (V0 / V) ** (2/3)
    energy = E0 + (9 * V0 * B0_internal / 16) * (
                (eta - 1)**3 * B0p + 
                (eta - 1)**2 * (6 - 4 * eta)
            )
    return energy

# 2. Load your data
data = np.loadtxt('eV_data.dat') # Ensure file has no comment lines, or use skiprows=1
V = data[:, 0]   # First column is volume
E = data[:, 1]   # Second column is energy

# 3. Provide initial guess values [E0, V0, B0(GPa), B0p]
# This is very important! Good initial values help fitting converge.
# E0_guess: Minimum energy in data
# V0_guess: Volume corresponding to minimum energy
# B0_guess: Estimate based on material (e.g., iron~150 GPa, diamond~440 GPa)
# B0p_guess: Usually set to 4.0
initial_guess = [np.min(E), V[np.argmin(E)], 150, 4.0]

# 4. Perform fitting
popt, pcov = curve_fit(birch_murnaghan_3rd_order, V, E, p0=initial_guess, maxfev=8000)
# popt: Optimal parameters [E0_opt, V0_opt, B0_opt, B0p_opt]
# pcov: Covariance matrix for calculating errors

# 5. Calculate standard errors of parameters
perr = np.sqrt(np.diag(pcov))

# 6. Output results
E0_opt, V0_opt, B0_opt_GPa, B0p_opt = popt
E0_err, V0_err, B0_err, B0p_err = perr

print("=== Birch-Murnaghan 3rd Order Fitting Results ===")
print(f"E0 = {E0_opt:.9f} ± {E0_err:.9f} (Ry)")
print(f"V0 = {V0_opt:.6f} ± {V0_err:.6f} (a.u.^3)")
print(f"B0 = {B0_opt_GPa:.6f} ± {B0_err:.6f} (GPa)")
print(f"B0' = {B0p_opt:.6f} ± {B0p_err:.6f}")

# 7. Plot fitted curve and data points
V_fit = np.linspace(np.min(V), np.max(V), 1000)
E_fit = birch_murnaghan_3rd_order(V_fit, *popt)

plt.figure(figsize=(10, 6))
plt.scatter(V, E, color='blue', label='WIEN2k Data', s=20)
plt.plot(V_fit, E_fit, 'r-', label=f'B-M Fit: $B_0$={B0_opt_GPa:.2f} GPa')
plt.xlabel('Volume (a.u.$^3$)')
plt.ylabel('Energy (Ry)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('E-V Curve Fitting with Birch-Murnaghan EOS')
plt.savefig('BM_fit.png', dpi=300)
plt.show()

Unit conversion note: The conversion factor 14710.5 in the above code applies when energy is in Ry and volume is in Bohr\(^3\) (a.u.\(^3\)). If your data units are different (such as eV and \(\AA^3\)), you need to use different conversion factors. This is one of the most error-prone areas in fitting.

Method 2: Using OriginLab or Gnuplot

These graphical software also support nonlinear fitting.

  1. Prepare data: Import V and E two-column data into software.
  2. Input equation: In nonlinear fitting tool, input third-order B-M equation formula.

  3. Origin: Create a new function in Fitting Function Organizer, input a formula like the following (variable is x, parameters are E0, V0, B0, B0p):

    1
    E0 + (9*V0*B0/16) * ( (((V0/x)^(2/3)-1)^3 * B0p ) + (((V0/x)^(2/3)-1)^2 * (6-4*(V0/x)^(2/3)) )
    - Gnuplot:

    1
    birch(x, E0, V0, B0, B0p) = E0 + (9*V0*B0/16) * ( (((V0/x)**(2./3.)-1)**3 * B0p ) + (((V0/x)**(2./3.)-1)**2 * (6-4*(V0/x)**(2./3.)) )
    3. Set initial parameters: Similar to Python code in Method 1, provide reasonable initial guess values (E0≈E_min, V0≈V@E_min, B0≈100-200, B0p≈4.0). 4. Execute fitting: Run iterative fitting, software will automatically find optimal parameters and their errors.

5. Results Analysis and Validation

  1. Graphical check: Always plot fitted curve and original data points! Ensure the curve passes through data points well, especially near the minimum (\(V_0\)).
  2. Parameter reasonableness:
  3. $V_0$ should be within your calculated volume range.
  4. $B_0$ value should match known physical properties of the material you're studying (compare with literature). For example, ordinary metals typically have \(B_0\) between 50-200 GPa.
  5. $B_0'$ is typically between 3.5 and 5.5.
  6. Error assessment: Check standard errors of parameters given by fitting software. If a parameter's error is very large (e.g., error value comparable to parameter value itself), it indicates insufficient data to accurately determine that parameter, or initial guess value is too poor.
  7. Convergence: Ensure your WIEN2k calculations have sufficiently converged (k-point grid, RKMAX, SCF iterations), otherwise data point noise will be large, affecting fitting accuracy.

6. Common Problems and Solutions

  • Fitting doesn't converge or results are unreasonable:
  • Problem: Initial guess values (p0) are too far from true values.
  • Solution: Provide better initial values. V0 and E0 can be estimated directly from data. B0 can try a reasonable range (50, 100, 150...).
  • B0 error is very large:
  • Problem: Too few data points or volume range too narrow to constrain \(B_0\) (bulk modulus is second derivative of energy, requires high data precision).
  • Solution: Add more volume points for calculation, especially near equilibrium volume \(V_0\).
  • Unit errors:
  • Problem: B0 fitting result value is 0.01 or millions, completely unreasonable.
  • Solution: Carefully check and unify all units (energy, volume), and perform correct unit conversion in equations. This is the most common problem!

By following these steps, you can reliably extract material equilibrium properties from WIEN2k calculation results.