2.1. Si1−xGex alloy

For pure crystals such as Si or Ge, EOS is often sufficient because atoms remain on high-symmetry sites. In alloys like \(\mathrm{Si_{0.7}Ge_{0.3}}\), EOS must be combined with ionic relaxation, since Si–Si, Si–Ge, and Ge–Ge bonds have different equilibrium lengths, producing local strain.

In this tutorial we demonstrate how to perform EOS and atomic relaxation for a \(2\times 2\times 2\) SiGe alloy supercell using RESCU.

../../../../_images/sige.png

Fig. 2.1.1 Example of a \(\mathrm{Si_{0.7}Ge_{0.3}}\) supercell (\(2\times 2\times 2\) conventional cubic).

The calculation steps are given as follows:

  1. Generate alloy structures

    • Create a Special Quasirandom Structure (SQS) for the SiGe alloy.

    • Generate hydrostatically scaled versions of the SQS for the EOS scan.

  2. Run Equation of State (EOS) scan

    • Run SCF calculations for each scaled structure.

    • Fit the total energy vs. lattice constant data to find the equilibrium lattice constant \(a_0\).

  3. Perform geometry optimization

    • Fix the lattice to the equilibrium value \(a_0\).

    • Relax the internal atomic positions to minimize forces.

    • Analyze bond length distributions before and after relaxation.

2.1.1. Project Layout

We recommend organizing the project as:

eos_relaxation/
├── generate_sige_alloy.py
├── generate_sige_alloy_sqs.py
├── generate_sige_aEOS_alloy.py
├── xyz/
├── eos/
│   ├── scf_base.input
│   ├── run_eos_scan.py
│   └── fit_eos.py
└── relaxation/
    ├── relax_base.input
    ├── run_relax.py
    └── analyze_bonds.py

  • The xyz/ directory stores alloy structures.

  • EOS and relaxation each have their own subfolder.

  • Python helper scripts automate scans, fitting, and bond-length analysis.

2.1.2. Equation of State (EOS)

The EOS procedure consists of varying the lattice constant around an initial guess (e.g., Vegard’s law) and computing the total energy for each scaled volume. Plotting \(E(V)\) yields a nearly parabolic curve around the minimum. Fitting this curve allows extraction of the equilibrium lattice constant \(a_{0}\).

Step-by-step procedure

2.1.2.1. Generate the alloy structure

To model the disordered alloy Si1−xGex, we construct a 2×2×2 supercell (64 atoms) of the conventional cubic diamond cell with a Ge concentration of \(x = 0.3125\) (20 Ge atoms out of 64).

Rather than randomly assigning Ge atoms to lattice sites, we use a Special Quasirandom Structure (SQS) to statistically reproduce the local chemical environment of a random alloy, as defined by correlation functions of atom pairs, triplets, etc. This improves accuracy and reproducibility over simple random substitution.

The method is based on the seminal work of Wei, Ferreira, Bernard, and Zunger [WEI1990] and implemented using the sqsgenerator package [SQSGENERATOR2023]. This performs a Monte Carlo optimization to find the best atomic occupations that match the target composition and correlation shells.

Use the following script:

This generates:

  • xyz/SiGe_2x2x2_sqs_S90.xyz … _S110.xyz – per-scale hydrostatic variants for EOS (cell and atomic positions scaled consistently).

Optionally, for comparison or rapid prototyping, you may use random substitution:

This method simply shuffles atomic species using a fixed seed. However, it lacks the statistical realism and stability of the SQS method.

2.1.2.2. Prepare base input file

To begin the EOS scan, we create a base input file with placeholders for different lattice constants.

For this calculation, you will need the following pseudopotential files:

  • Si_DZP_PBE.mat

  • Ge_sp_SZP_PBE.mat

SCF input file

Listing 2.1.1 eos/scf_base.input
 units.length = 'Angstrom';
 LCAO.status = true;
 info.calculationType = 'self-consistent';

 info.savepath = './results/scf_TAG';
 info.outfile  = './results/scf_TAG.out';

 % Elements / pseudopotentials
 element(1).species = 'Si';
 element(1).path    = './Si_DZP_PBE.mat';
 element(2).species = 'Ge';
 element(2).path    = './Ge_sp_SZP_PBE.mat';

 % Lattice constant from Vegard's law (Å) and isotropic scaling
 a0 = A0;
 fac = SCALE / 100;
 a  = a0 * fac;

 % 2×2×2 conventional cubic supercell (64 atoms total)
 domain.latvec = [a*2, 0, 0; 0, a*2, 0; 0, 0, a*2];
 units.domain.latvec = 'Angstrom';

 % Atomic positions (XYZ, Å)
 atom.xyz       = 'XYZ_FILE';
 units.atom.xyz = 'Angstrom';

 % Exchange-correlation
 functional.list  = {'XC_GGA_X_PBE', 'XC_GGA_C_PBE'};
 functional.libxc = true;

 % k-points
 kpoint.gridn = [3, 3, 3];

 % Numerical settings
 domain.lowres = 0.14;

The keywords specify the following:

For a detailed description of most parameters used here, such as LCAO.status, info.calculationType, info.savepath, element(1).species, functional.list, kpoint.gridn, and domain.lowres, please refer to the Numerical Convergence tutorial.

Note

  • A0 and SCALE: The global volume is controlled by the cubic lattice parameter \(a = \mathrm{SCALE}\times A0/100\). A0 is obtained from Vegard’s law (linear interpolation of Si and Ge lattice constants at the chosen composition), a standard first guess for random alloys [VEGARD1921]. We then scan around it with several SCALE values (e.g., 90, 95, 100, 105, 110).

  • domain.latvec: Constructs a 2×2×2 conventional cubic supercell, i.e., the supercell edge is \(2a\). EOS varies only the volume (global, hydrostatic-like change), keeping the crystal shape cubic.

  • atom.xyz: Points to the per-scale XYZ file generated earlier (e.g., SiGe_2x2x2_sqs_S95.xyz), where both the cell and atomic positions are scaled consistently. Using per-scale XYZs ensures RESCU reads a geometry consistent with the volume setting in domain.latvec.

  • kpoint.gridn: A moderate 3×3×3 grid is adequate for an EOS scan. We recommend testing k-point and grid settings as shown in the Numerical Convergence tutorial and re-using those here. If your converged settings differ, update this accordingly.

  • domain.lowres: Real-space resolution (smaller is finer). Here we use 0.14 Å as a production-quality setting consistent with our convergence guidance. If you have not converged this yet, start with a slightly coarser value to map the EOS cheaply, then re-run at the converged value near the minimum.

  • info.savepath / info.outfile and TAG: Each EOS point writes into a unique directory/file pair labeled by TAG (set by the driver script), making post-processing straightforward.

Tip

EOS in RESCU uses total energies only (no stress tensor). We therefore fit \(E(a)\) near its minimum with a low-order polynomial to extract the equilibrium lattice constant \(a_{0}\) reliably, then perform the ionic relaxation at that fixed lattice (next section).

2.1.2.3. Run EOS scan

Run multiple scaled structures (e.g., SCALE = 90, 95, 100, 105, 110 around the Vegard’s law guess).

Download: run_eos_scan.py

cd eos
python run_eos_scan.py

This script:

  1. Computes A0 from Vegard’s law for \(x_{\rm Ge}=0.3125\) (20 Ge in 64 atoms).

  2. Checks that all per-scale XYZ files exist (SiGe_2x2x2_sqs_S90.xyz … _S110.xyz).

  3. Renders one input per SCALE by replacing A0, SCALE, XYZ_FILE and TAG.

  4. Launches RESCU for each case using mpiexec (see script header to adjust paths and MPI ranks).

2.1.2.4. Analyze EOS data and fit the minimum

Download: fit_eos.py

cd eos
python fit_eos.py

The script collects total energies from results/scf_*.mat (and/or *.out), extracts the cubic lattice constant \(a\) for each SCALE, and fits a cubic (default) \(E(a)\) near the minimum. The equilibrium lattice parameter \(a_{0}\) is taken from the fitted minimum. You may decrease the polynomial degree to 2 if the data points are very symmetric around the minimum.

Outputs:

  • CSV table of SCALE, \(a\), \(E\) (cell and per atom) and \(\Delta E\) vs the best point (placeholder: eos_data.csv).

  • Plot of \(E(a)\) with the cubic fit and the fitted minimum marker (placeholder: e_vs_a.png, Fig. 2.1.2).

  • Text report summarizing the fitted \(a_{0}\) and energy minimum (placeholder: eos_fit_report.txt).

../../../../_images/e_vs_a.png

Fig. 2.1.2 Placeholder: EOS energy vs. lattice constant with cubic fit and fitted \(a_{0}\).

Note

Ensure your SCALE range brackets the minimum (e.g., 90–110%). If curvature appears weak or the minimum sits at an edge of your scan, add intermediate points (e.g., 97, 98, 102, 103) and refit.

2.1.3. Geometry Optimization (Ionic Relaxation)

Once the equilibrium lattice constant is obtained from EOS, we fix the cubic cell and relax the internal atomic positions to minimize forces. In RESCU, this is done with info.calculationType = 'relaxation' and the sr.* controls for ionic optimization. See the RESCU relaxation tutorial for details.

Key sr.* options:

  • sr.method — optimization algorithm (e.g., 'cg')

  • sr.tol — force tolerance (eV/Å)

  • sr.maxStep — maximum ionic steps

  • sr.relaxedstructure — output XYZ file written when converged

  • sr.moveableAtomList — optional list of atoms allowed to move (not used here)

Step-by-step procedure

2.1.3.1. Prepare relaxation structure

Before running the relaxation, we must generate an XYZ file where the SQS structure is hydrostatically scaled to the equilibrium lattice constant \(a_{\rm EOS}\) found from the EOS fit.

Download the script:

Edit the script to set a_eq to the value you obtained from fit_eos.py (e.g., 5.538145). Then run it:

python generate_sige_aEOS_alloy.py

This will create xyz/SiGe_2x2x2_sqs_aEOS.xyz, which will be used as the input geometry for the relaxation calculation.

2.1.3.2. Prepare input file

We relax at the fixed EOS lattice constant \(a_{\rm EOS} = 5.538145\) Å using the SQS structure created earlier at the same lattice. The input below uses two placeholders that the helper script fills:

  • A0 — equilibrium lattice constant (Å), here 5.538145

  • XYZ_FILE — path to the pre-relaxation coordinates (here the SQS at \(a_{\rm EOS}\))

For this calculation, you will need the following pseudopotential files:

  • Si_DZP_PBE.mat

  • Ge_sp_SZP_PBE.mat

Relaxation input

Listing 2.1.2 relaxation/relax_base.input
units.length = 'Angstrom';
LCAO.status = true;
info.calculationType = 'relaxation';

info.savepath = './results/relax_sige_sqs_aEOS';
info.outfile  = './results/relax_sige_sqs_aEOS.out';

% Elements / pseudopotentials
element(1).species = 'Si';
element(1).path    = './Si_DZP_PBE.mat';
element(2).species = 'Ge';
element(2).path    = './Ge_sp_SZP_PBE.mat';

% Fixed 2×2×2 cubic supercell at the EOS lattice (Å)
a = A0;  % filled by run_relax.py (here: 5.538145)
domain.latvec = [2*a, 0, 0; 0, 2*a, 0; 0, 0, 2*a];
units.domain.latvec = 'Angstrom';

% Atomic positions at a = A0 (Å) — pre-relaxation structure
atom.xyz       = 'XYZ_FILE';  % ../xyz/SiGe_2x2x2_sqs_aEOS.xyz
units.atom.xyz = 'Angstrom';

% Exchange-correlation
functional.list  = {'XC_GGA_X_PBE', 'XC_GGA_C_PBE'};
functional.libxc = true;

% k-points (use your converged mesh from numerical convergence)
kpoint.gridn = [3, 3, 3];

% Numerical settings (use your converged real-space resolution)
domain.lowres = 0.14;

% SCF settings
spin.type = 'degenerate';
option.maxSCFiteration = 200;

% --- Structural relaxation (ionic) ---
sr.method           = 'cg';                % conjugate gradient
sr.tol              = 0.01;               % eV/Ang — force threshold
sr.maxStep          = 200;
sr.relaxedstructure = './results/SiGe_2x2x2_sqs_aEOS_relaxed.xyz';

2.1.3.3. Run relaxation

Download the run helper and execute it:

cd relaxation
python run_relax.py

What the helper does:

  1. Sets A0 = 5.538145 Å (from the EOS fit) and points XYZ_FILE to ../xyz/SiGe_2x2x2_sqs_aEOS.xyz.

  2. Renders relax.input from relax_base.input and runs RESCU with mpiexec.

  3. Writes the relaxed coordinates to results/SiGe_2x2x2_sqs_aEOS_relaxed.xyz (as specified by sr.relaxedstructure).

Tip

During visualization, atoms can appear to “jump” because of periodic boundary conditions (wrapping across cell faces). To quantify changes robustly, analyze bond-length distributions rather than raw displacements.

2.1.3.4. Analyze results

Use the analysis script to compare first-neighbor bond lengths before vs after relaxation.

cd relaxation
python analyze_bonds.py

The script reads:

  • Before: ../xyz/SiGe_2x2x2_sqs_aEOS.xyz

  • After: results/SiGe_2x2x2_sqs_aEOS_relaxed.xyz

and produces:

A text summary to stdout with a table of mean bond lengths per species pair, before and after relaxation.

Expected outcome:

  • Before relaxation, all bonds cluster near the uniform EOS value (~2.398 Å).

  • After relaxation, the distributions separate as Si–Si < Si–Ge < Ge–Ge, e.g. means around ~2.37, ~2.41, ~2.47 Å, reflecting the different preferred bond lengths of Si and Ge environments.

Note

We use a 3.0 Å cutoff to isolate first neighbors in diamond-like SiGe (second neighbors are ~3.9 Å). If you change the lattice or chemistry, adjust the cutoff accordingly.

2.1.4. Results

The equilibrium lattice constant for the Si0.6875Ge0.3125 alloy obtained from the EOS fit is 5.538 Å. This value is in good agreement with experimental and other theoretical results, as shown in the table below.

Table 2.1.1 Comparison of Si0.6875Ge0.3125 lattice constant.

Method

Lattice Constant (Å)

This work (RESCU, PBE)

5.538

Vegard’s Law [VEGARD1921]

5.501

Experiment [DISMUKES1964]

5.493

The analysis of bond lengths before and after ionic relaxation demonstrates the local structural adjustments within the alloy. Before relaxation, all bond lengths are identical, corresponding to a uniform scaling of the ideal diamond lattice. After relaxation, the bond lengths diverge to accommodate the different atomic sizes of Si and Ge.

Table 2.1.2 Average bond lengths before and after relaxation.

Bond Type

Before Relaxation (Å)

After Relaxation (Å)

Experiment (Å) [MENONI1992], [AUBRY2001]

Si-Si

2.398

2.373

2.36

Si-Ge

2.398

2.411

2.40

Ge-Ge

2.398

2.470

2.45

The calculated bond lengths are in good agreement with experimental values obtained from EXAFS (Extended X-ray Absorption Fine Structure) and XAFS measurements. The trend where Si-Si bonds are shortest, Ge-Ge bonds are longest, and Si-Ge bonds are intermediate is correctly reproduced. This confirms that the geometry optimization has successfully captured the local strain relaxation in the SiGe alloy.

2.1.5. References

[PAYNE1992]

M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev. Mod. Phys. 64, 1045 (1992).

[WEI1990]

S.-H. Wei, L. G. Ferreira, J. E. Bernard, and A. Zunger, Phys. Rev. B 42, 9622–9649 (1990), https://doi.org/10.1103/PhysRevB.42.9622

[SQSGENERATOR2023]

D. Gehringer et al., Computer Physics Communications 289, 108791 (2023), https://doi.org/10.1016/j.cpc.2023.108791

[VEGARD1921] (1,2)

L. Vegard, Die Konstitution der Mischkristalle und die Raumfüllung der Atome. Z. Phys. 5, 17 (1921).

[DISMUKES1964]

J. P. Dismukes, L. Ekstrom, and R. J. Paff. Lattice Parameter and Density in Germanium-Silicon Alloys. J. Phys. Chem. 68, 3021 (1964).

[MENONI1992]

C. S. Menoni, J. B. Boyce, B. M. Clemens, J. F. Morar, R. H. Kappauf, and R. J. Nemanich. Bond-length relaxation in crystalline Si:sub:`1−x`Ge:sub:`x` alloys: An EXAFS study. Phys. Rev. B 45, 14005 (1992). https://doi.org/10.1103/PhysRevB.45.14005

[AUBRY2001]

S. Aubry, D. A. Drabold, and P. Ordejón. Strain relaxation mechanisms and local geometric/electronic structure in Si:sub:`1−x`Ge:sub:`x` alloys. Phys. Rev. B 64, 165205 (2001). https://doi.org/10.1103/PhysRevB.64.165205