The usage of maximum_distance_factor in LocalGeometryFinder.compute_structure_environments

[Jeff-oakley]

Something strange happens when I use LocalGeometryFinder.compute_structure_environments. The maximum_distance_factor is explained as "param maximum_distance_factor: If not set to None, neighbors beyond maximum_distance_factor*closest_neighbor_distance are not considered". It appears to be the maximum bond length when evaluating the bonding environment. Usually, I will think a larger enough number would at least give some results. However, for some structures, this will raise up error message due to the failure of calculating valid ce_dict. The following code can reproduce some frequently seen error messages. Could anyone let me know why this happens? Also, how do we choose maximum_distance_factor normally?

from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import LocalGeometryFinder
from pymatgen.ext.matproj import MPRester
from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import SimplestChemenvStrategy, MultiWeightsChemenvStrategy
from pymatgen.analysis.chemenv.coordination_environments.structure_environments import LightStructureEnvironments
import numpy as np

qe = MPRester('v8mU74QhZSIisN26')
probe_str = qe.get_entries('mp-22915',inc_structure=True)[0].structure
d_cut = np.min(probe_str.distance_matrix[~np.eye(probe_str.distance_matrix.shape[0],dtype=bool)]) * 1.4

for max_df in [1.0,1.1,1.2,2.0]:
    print(f'Test with maximum_distance_factor={max_df}')
    lgf = LocalGeometryFinder()
    lgf.setup_structure(structure=probe_str)
    se = lgf.compute_structure_environments(maximum_distance_factor=max_df,only_cations=False)
    strategy = SimplestChemenvStrategy(distance_cutoff=d_cut)
    lse = LightStructureEnvironments.from_structure_environments(strategy=strategy, structure_environments=se)

    print('{} neighbors_sets exist'.format(len(lse.neighbors_sets)))
    coord_lst = {}
    cation, anion = 'Ag', 'I'
    cation_inds = probe_str.indices_from_symbol(cation)
    anion_inds = probe_str.indices_from_symbol(anion)
    for ind, nb_set in enumerate(lse.neighbors_sets):
        if ind not in cation_inds: continue
        coord_lst[ind]=[]
        try:
            for dict_val in nb_set[0].neighb_sites_and_indices:
                if dict_val['index'] in anion_inds:
                    coord_lst[ind].append(dict_val['index'])
        except: continue

    cation_n_coords = [len(nn_ind_lst) for ind, nn_ind_lst in coord_lst.items()]
    print('Cation {} has neighbors of {}'.format(cation, cation_n_coords))

[JaGeo]

I am suspecting that this is related to stored reference polyhedra in ChemEnv. ChemEnv does not have any reference parameters for polyhedra with more than 13 neighbors. Starting from a max_distance_cutoff of 1.2, you get 14 however. I will further look into your case and hope to have a better answer soon.

[JaGeo]

To circumvent this, you can focus on cation-anion bonds or cation-cation bonds only. For example by providing valences:

import numpy as np

from pymatgen.analysis.bond_valence import BVAnalyzer
from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import SimplestChemenvStrategy
from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import LocalGeometryFinder
from pymatgen.analysis.chemenv.coordination_environments.structure_environments import LightStructureEnvironments
from pymatgen.analysis.chemenv.utils.defs_utils import AdditionalConditions
from pymatgen.ext.matproj import MPRester

qe = MPRester()
probe_str = qe.get_entries('mp-22915', inc_structure=True)[0].structure
valences = BVAnalyzer().get_valences(probe_str)
d_cut = np.min(probe_str.distance_matrix[~np.eye(probe_str.distance_matrix.shape[0], dtype=bool)]) * 1.4

for max_df in [1.0, 1.1, 1.2, 1.4, 2.0]:
    print(f'Test with maximum_distance_factor={max_df}')
    lgf = LocalGeometryFinder()
    lgf.setup_structure(structure=probe_str)
    se = lgf.compute_structure_environments(maximum_distance_factor=max_df, only_cations=False, valences=valences,
                                            additional_conditions=[AdditionalConditions.ONLY_ACB])
    strategy = SimplestChemenvStrategy(distance_cutoff=d_cut)
    lse = LightStructureEnvironments.from_structure_environments(strategy=strategy, structure_environments=se)

    print('{} neighbors_sets exist'.format(len(lse.neighbors_sets)))
    coord_lst = {}
    cation, anion = 'Ag', 'I'
    cation_inds = probe_str.indices_from_symbol(cation)
    anion_inds = probe_str.indices_from_symbol(anion)
    for ind, nb_set in enumerate(lse.neighbors_sets):
        if ind not in cation_inds: continue
        coord_lst[ind] = []
        try:
            for dict_val in nb_set[0].neighb_sites_and_indices:
                if dict_val['index'] in anion_inds:
                    coord_lst[ind].append(dict_val['index'])
        except:
            continue

    cation_n_coords = [len(nn_ind_lst) for ind, nn_ind_lst in coord_lst.items()]
    print('Cation {} has neighbors of {}'.format(cation, cation_n_coords))

With regard to choosing the parameter, I would use the default value but with valences.

[JaGeo]

It's on my list somewhere to update the tutorial in this regard but might take time. We should probably also add a clearer error-message.