2 questions on `determine_topology`

[kaueltzen]

Hi,
I have two questions on the return value of determine_topology.
The first one concerns its length / the number of substructures with disjoint vertices.

Is it expected behaviour that this depends on the number of lattice points of the input structure provided?
Intuitively, I would expect different representations of the same crystal structure to yield the same determine_topology result but maybe I misunderstood the documentation there.
I appended 2 examples that compare the determine_topology output of (a) a structure and its 2x2x2 superstructure and (b) a centered structure and its primitive.

(CrystalNets.jl v1.1.1)

import juliacall
from pymatgen.core import Structure, Lattice
from pymatgen.analysis.local_env import NearNeighbors, JmolNN
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer


jl = juliacall.newmodule("Topo")
jl.seval("using CrystalNets")


options_input = jl.CrystalNets.Options(bonding=jl.CrystalNets.Bonding.Input,
                                       export_input=False, 
                                       export_subnets=False, 
                                       export_net=False)

def write_bonded_cif(nn: NearNeighbors, structure: Structure, file_name: str):
    bonds = []
    for site_id, site in enumerate(structure.sites):
        nbs = nn.get_nn_info(structure=structure, n=site_id)
        for nb in nbs:
            site_str = f"{site.species_string}{site_id}"
            site_to_str = f"{nb['site'].species_string}{nb['site_index']}"
            length = nb["site"].nn_distance
            bond_str = f"{site_str} {site_to_str} {length}"
            reverse_bond_str = f"{site_to_str} {site_str} {length}"
            if bond_str not in bonds and reverse_bond_str not in bonds:
                bonds.append(bond_str)
    
    topo_strings = ["loop_", "_geom_bond_atom_site_label_1", "_geom_bond_atom_site_label_2", "_geom_bond_distance"] +bonds
    structure.to_file(file_name)

    
    with open(file_name, "a") as file:
        for line in topo_strings:
            file.write(line + "\n")

# Example structure and 2x2x2 superstructure
structure = Structure.from_spacegroup(
                    "Pmmm", Lattice.orthorhombic(3.3, 7, 7), ["Po"], [[0, 0, 0],]
                )
super_structure = structure.make_supercell(scaling_matrix=[2, 2, 2], in_place=False)

file_name_simple = "structure.cif"
write_bonded_cif(nn=JmolNN(), structure=structure, file_name=file_name_simple)
res = jl.determine_topology(file_name_simple, options_input
                                             )
print("simple topo: ", len(res))
print(res)

file_name_super = "superstructure.cif"
write_bonded_cif(nn=JmolNN(), structure=super_structure, file_name=file_name_super)
super_res = jl.determine_topology(file_name_super, options_input
                                                           )
print("2x2x2 superstructure topo: ", len(super_res))
print(super_res)

Output:
simple topo: 1
1-rod (1 1 1 1)
2x2x2 superstructure topo: 4
4 interpenetrated substructures:
⋅ Subnet 1 → 1-rod (1 1 1 1)
⋅ Subnet 2 → 1-rod (1 1 1 1)
⋅ Subnet 3 → 1-rod (1 1 1 1)
⋅ Subnet 4 → 1-rod (1 1 1 1)

# Example with centered vs. primitive lattice
centered_structure = Structure.from_spacegroup(
                    "Cmmm", Lattice.orthorhombic(3.3, 9, 9), ["Po"], [[0, 0, 0]]
                )
sga_centered = SpacegroupAnalyzer(structure=centered_structure)
primitive_structure = sga_centered.get_primitive_standard_structure()

file_name_centered = "centered_structure.cif"
write_bonded_cif(nn=JmolNN(), structure=centered_structure, file_name=file_name_centered)
centered_res = jl.determine_topology(file_name_centered, options_input
                                             )
print("centered structure topo: ", len(centered_res))
print(centered_res)

file_name_primitive = "primitive_structure.cif"
write_bonded_cif(nn=JmolNN(), structure=primitive_structure, file_name=file_name_primitive)
primitive_res = jl.determine_topology(file_name_primitive, options_input
                                                           )
print("primitive structure topo: ", len(primitive_res))
print(primitive_res)

Output:
centered structure topo: 2
2 interpenetrated substructures:
⋅ Subnet 1 → 1-rod (1 1 1 1)
⋅ Subnet 2 → 1-rod (1 1 1 1)
primitive structure topo: 1
1-rod (1 1 1 1)

---

My second question is: is it possible to get the full topological genome of a crystal structure that has multiple substructures with disjoint vertices, not only the genomes of the substructures?

Thanks in advance!

[Liozou]

For the first question, the answer is that it is indeed expected, because we differentiate two kinds of interpenetration. That may not be very intuitive so I'll try to illustrate this with schematic pictures:

• First, there is the general case of two independent nets that overlap. Those are what we call the "substructures" in the manual and the topology of these substructures can be the same or not. This is represented below as the green and ochre nets, which happen to be of different topologies. If you removed any of these nets, it would not affect the others.
  
• Then, there is the case of catenated nets. We use this terminology to designate a situation in which each kind of vertex can be connected to each other kind (in the example below, it is possible to connect a circle to a triangle and to a square), but not all pairs of vertices can be connected (it is impossible to connect two vertices of different colors). One characteristic of this particular kind of interpenetration is that the interpenetrating nets are all the same, with a shift of an integer number of unit cells. Another is that the unit cell of a single component is larger than the initial one: if you only kept the black component, its repeating pattern does not fit in the red square that represents the unit cell of the overall topology.
  

The type returned by determine_topology is the list of pairs (substructure, n) where n is the number of catenation of that substructure (see the manual for more details). In the catenation example above (second picture), if the unit cell was twice as large, each component would be considered as an independent substructure of the whole. This is represented in the picture below:

Like in your examples, the last two pictures represent the same crystal but the difference of unit cell makes a difference in the nature of the interpenetration, hence the difference in the way CrystalNets reports them. Of course, each subtopology remains the same.

---

For the second question: we only define the topological genome of connected periodic graph, so there is no way to get a single genome for a crystal structure with interpenetration. The reason is that, in many cases, what is of interest in the case of interpenetration is not only the nature of the subtopologies, but also the way they are knotted together (see for example https://doi.org/10.1039/B102400K) which we cannot tackle with topological tools.
However, depending on what you are looking for, you may be interested in the function CrystalNets.total_interpenetration which takes the value returned by determine_topology and splits it by topology to access the number of interpenetrated components with that topology. If your structure only has a single subtopology, you may alternatively be interested in CrystalNets.one_topoloy, with the same interface, that returns the only topology in the value returned by determine_topology and fails if there is more than one.

[kaueltzen]

Thank you very much for your detailed answer @Liozou !