SQS Generator Error "RuntimeWarning: divide by zero encountered in long_scalars" with specific vacancy concentrations

Hi all! Relatively new icet user here. I’m using the generate_sqs to find a structure with disordered oxygen vacancies. I’ve run into the above error, “divide by zero encountered in long_scalars” causing the command to quit for specific vacancy concentrations. I can get my script, included below to run with 1/4, 1/3 and 1/15 oxygen vacancies. For 1/12 (the concentration needed) I get a reproducible divide by zero error. Note that 1/12 vacancies work for other structures… I’m really running into a bit of a wall! Any help and advice is much appreciated!

Below is my script I wrote based on icet examples provided, and then the oxygen sublattice cif.

from ase import Atom
from ase.build import bulk
from ase.io import write, read
from icet import ClusterSpace
from icet.tools.structure_generation import (generate_sqs,generate_sqs_from_supercells, generate_sqs_by_enumeration,generate_target_structure)
from icet.input_output.logging_tools import set_log_config
set_log_config(level='INFO')
#Build structure type.
#create a primitive cell with just oxygen, X is for vacancies
primitive = read('Pnma_matrix_for_SQS_primitive-oxygens.cif')
#Set number of different elements. Add or change if necessary
cs = ClusterSpace(primitive, [8.0, 4.0], ['O', 'X'])
#Set concentration.
#Decimals should be to 7th decimal to avoid rounding errors like 0.3333333
target_concentrations = {'O': 11/12, 'X': 1/12}
#Call generate_sqs to create sqs. Change size or sqs generator based on icet documentation
sqs = generate_sqs(cluster_space=cs,max_size=12,include_smaller_cells=False,target_concentrations=target_concentrations)
#create a cif for the sqs Atom object.
#Edit name if necessary, as this file will be overwritten
write('test-sqs.cif',sqs) #edit name

Oxygen sublattice cif.


#======================================================================
# CRYSTAL DATA
#----------------------------------------------------------------------
data_VESTA_phase_1

_chemical_name_common                  'cesium lead tribromide - gamma'
_cell_length_a                         8.210500
_cell_length_b                         8.282900
_cell_length_c                         11.803900
_cell_angle_alpha                      90.000000
_cell_angle_beta                       90.000000
_cell_angle_gamma                      90.000000
_cell_volume                           802.744829
_space_group_name_H-M_alt              'P 1'
_space_group_IT_number                 1

loop_
_space_group_symop_operation_xyz
   'x, y, z'

loop_
   _atom_site_label
   _atom_site_occupancy
   _atom_site_fract_x
   _atom_site_fract_y
   _atom_site_fract_z
   _atom_site_adp_type
   _atom_site_U_iso_or_equiv
   _atom_site_type_symbol
   O1         1.0     0.961400     0.473200     0.250000    Uiso  ? O
   O2         1.0     0.038600     0.526800     0.750000    Uiso  ? O
   O3         1.0     0.461400     0.026800     0.750000    Uiso  ? O
   O4         1.0     0.538600     0.973200     0.250000    Uiso  ? O
   O5         1.0     0.795100     0.206100     0.023600    Uiso  ? O
   O6         1.0     0.204900     0.793900     0.976400    Uiso  ? O
   O7         1.0     0.295100     0.293900     0.976400    Uiso  ? O
   O8         1.0     0.704900     0.706100     0.023600    Uiso  ? O
   O9         1.0     0.204900     0.793900     0.523600    Uiso  ? O
   O10        1.0     0.795100     0.206100     0.476400    Uiso  ? O
   O11        1.0     0.704900     0.706100     0.476400    Uiso  ? O
   O12        1.0     0.295100     0.293900     0.523600    Uiso  ? O

Full structure cif

# CIF file created by FINDSYM, version 7.1.3

data_findsym-output
_audit_creation_method FINDSYM

_cell_length_a    8.2829000000
_cell_length_b    11.8039000000
_cell_length_c    8.2105000000
_cell_angle_alpha 90.0000000000
_cell_angle_beta  90.0000000000
_cell_angle_gamma 90.0000000000
_cell_volume      802.7448816368

_symmetry_space_group_name_H-M "P 21/n 21/m 21/a"
_symmetry_Int_Tables_number 62
_space_group.reference_setting '062:-P 2ac 2n'
_space_group.transform_Pp_abc a,b,c;0,0,0

loop_
_space_group_symop_id
_space_group_symop_operation_xyz
1 x,y,z
2 x+1/2,-y+1/2,-z+1/2
3 -x,y+1/2,-z
4 -x+1/2,-y,z+1/2
5 -x,-y,-z
6 -x+1/2,y+1/2,z+1/2
7 x,-y+1/2,z
8 x+1/2,y,-z+1/2

loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_symmetry_multiplicity
_atom_site_Wyckoff_label
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
_atom_site_fract_symmform
La1 La   4 c 0.50610  0.25000 0.50400 1.00000 Dx,0,Dz
Ni1 Ni   4 b 0.00000  0.00000 0.50000 1.00000 0,0,0
O1  O    4 c -0.02680 0.25000 0.46140 1.00000 Dx,0,Dz
O2  O    8 d 0.70610  0.02360 0.29510 1.00000 Dx,Dy,Dz

Hi @Lauren_Walters,

Thank you for reporting! I think your system fails due to a bug that we have previously not encountered. Is it true that the code runs for a quite long time before failing? I didn’t have the patience to wait because the very long runtime in itself indicated that something was wrong. I have pushed a fix to a new branch. When I run that branch for your system it seems to work, but I need to test it a bit further before we merge with the master branch. For the time being you could try the new branch (or just replace icet/tools/structure_enumeration_support/normal_form_matrices.py with the corresponding file in that branch).

The worst case scenario if I did something stupid but the code still runs, is that not all potential supercells that are 12 times larger than the primitive cell are included in the search for SQS. That would only mean that the optimal SQS might not be found, but what you get would probably still be good given that you have many sites in your target SQS. The best SQS cells usually have roughly equally long axes, so a highly elongated cell would be a warning sign.