[Spglib-users] Can not find atoms related by space group symmetry

[Dorde D. <Dan...@ho...>]

Dear all,

I am having trouble finding atoms related by a spacegroup symmetry. I am not sure what I am doing wrong. I am sending a script that I am using to find symmetry-related atoms.

POSCAR:

  H
  1.00000000
  0.59175841  -0.59175841  1.52600353
  0.59175841  0.59175841  1.52600353
  -0.59175841  -0.59175841  1.52600353
  H
2
Direct
  0.00000000  0.00000000  0.00000000
  0.50000000  -0.25000000  0.25000000

test.py:

import numpy as np
import spglib
from phonopy.interface.calculator import read_crystal_structure

# Function to check if the vector is the translational lattice vector
def check_if_R(vec):
        rvec = False
        for i in range(-1,2):
                if not rvec:
                        for j in range(-1,2):
                                if not rvec:
                                        for k in range(-1,2):
                                                vec1 = vec + np.array([float(i), float(j), float(k)])
                                                if(np.linalg.norm(vec1) < 1.0e-5):
                                                        rvec = True
                                                        break
        return rvec

unitcell, _ = read_crystal_structure("POSCAR", interface_mode='vasp') # read in structure
cell = (unitcell.cell, unitcell.scaled_positions, unitcell.numbers)

dataset = spglib.get_symmetry_dataset(cell, symprec=1e-5, angle_tolerance=-1.0, hall_number=0) # get datasets
rotations = dataset['rotations'].copy()
translations = dataset['translations'].copy()

print(scaled_positions)

atom_image = np.zeros((len(rotations), len(numbers))) + len(numbers) # array to store to which atom we go with symmetry
for irot in range(len(rotations)):
        new_pos = np.einsum('ij,jk->ik', scaled_positions, rotations[irot]) + translations[irot] # symmetry operation on the scaled positions
        for iat in range(natom): # Try to match rotated atom with the original set
                found = False
                for jat in range(natom):
                        pos = new_pos[iat] - scaled_positions[jat]
                        pos -= np.rint(pos)
                        if(check_if_R(pos)):
                                atom_image[irot, iat] = jat
                                found = True
                                break
                if(not found): # Print rotated atom if could not match
                        print(irot, iat)
                        print(new_pos[iat])
# Print atom images
print(atom_image)

Could you please tell me what am I missing?

Thank you for your help.