[Python] How to correctly extract symmetry-ineqivalrnt atoms from Quantum Espresso output file...

I performed a vc-relax calculation in Quantum ESPRESSO on a structure that initially had 23 symmetry-inequivalent atoms. Due to the crystal symmetry, QE internally expanded this to 148 atoms. Unfortunately, the relaxation did not fully converge because the nstep limit was reached before the forces met the convergence threshold.

Now I want to continue the relaxation (or run a new calculation) using only the 23 inequivalent atoms relaxed out of 148 atoms, assuming they still represent a good approximation of the relaxed asymmetric unit. I used ASE to extract the first 23 atoms from vc-relax.out,Here is the error:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% task # 0 from check_atoms : error # 1 atoms # 97 and # 135 overlap! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
​

Here's the ASE code I used to extract the atoms from vc-relax.out and save the reduced structure for further DFT calculations:

import numpy as np
from ase.io import read
from ase.data import chemical_symbols
import spglib
from spglib import get_spacegroup_type

# --- Input: QE vc-relax output ---
pth = "vc-relax.out"
atoms = read(pth)

# --- Recover standard crystallographic lattice constants (Cmma, setting=2) ---
def get_standard_lattice_constants(atoms):
a_qe, b_qe, c_qe = atoms.cell.array
A_vec = a_qe - b_qe
B_vec = a_qe + b_qe
C_vec = c_qe
A = np.linalg.norm(A_vec)
B = np.linalg.norm(B_vec)
C = np.linalg.norm(C_vec)
return A, B, C

A, B, C = get_standard_lattice_constants(atoms)

# --- Use spglib to get symmetry information ---
cell = (atoms.cell.array, atoms.get_scaled_positions(), atoms.get_atomic_numbers())
sym_data = spglib.get_symmetry_dataset(cell, symprec=1e-5)

wyckoff_letters = sym_data['wyckoffs']
equiv_atoms = sym_data['equivalent_atoms']
unique_indices = sorted(set(equiv_atoms))

positions = atoms.get_scaled_positions()
atomic_numbers = atoms.get_atomic_numbers()

# Reduced site information
red_sites = []
red_Z = []
red_wyckoff = []
red_multiplicities = []

for i in unique_indices:
site = tuple(float(x) for x in positions)
red_sites.append(site)
red_Z.append(atomic_numbers)
wyck = wyckoff_letters
red_wyckoff.append(wyck)
multiplicity = list(wyckoff_letters).count(wyck)
red_multiplicities.append(multiplicity)

symbols = [chemical_symbols[z] for z in red_Z]

# --- Print summary to terminal ---
print("\nStandard lattice constants to use in QE:")
print(f"A = {A:.10f} Å")
print(f"B = {B:.10f} Å")
print(f"C = {C:.10f} Å")

print("\nSymmetry-inequivalent atoms:")
print("Element Wyckoff Multiplicity Fractional coordinates")
total_atoms = 0
for sym, wyck, mult, pos in zip(symbols, red_wyckoff, red_multiplicities, red_sites):
print(f"{sym:<8} {wyck:<8} {mult:<6} ({pos[0]:.6f}, {pos[1]:.6f}, {pos[2]:.6f})")
total_atoms += mult

print(f"\nTotal symmetry-inequivalent atoms: {len(red_sites)}")
print(f"Estimated total atoms in unit cell: {total_atoms}")

# --- Write to file in QE-style format ---
output_filename = "symmetry_inequivalent_atoms.txt"
with open(output_filename, "w") as f:
f.write("# Lattice constants (angstrom)\n")
f.write(f"A = {A:.10f}\n")
f.write(f"B = {B:.10f}\n")
f.write(f"C = {C:.10f}\n\n")
f.write("ATOMIC_POSITIONS crystal_sg\n")
for sym, pos in zip(symbols, red_sites):
f.write(f"{sym:<3} {pos[0]:.10f} {pos[1]:.10f} {pos[2]:.10f}\n")

print(f"\nOutput written to {output_filename}")



My question: Is there a correct and reliable way to extract only the relaxed symmetry-inequivalent atoms from Quantum ESPRESSO output for further calculations, without causing atom overlap or symmetry errors?

Continue reading...