File: explicit_permutations_plane_algorithm.py

package info (click to toggle)
pymatgen 2025.2.18%2Bdfsg1-5
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 85,888 kB
  • sloc: python: 176,173; javascript: 780; makefile: 221; sh: 78
file content (163 lines) | stat: -rw-r--r-- 6,202 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
"""
Development script of the ChemEnv utility to get the explicit permutations for coordination environments identified
with the separation plane algorithms (typically with coordination numbers >= 6).
"""

from __future__ import annotations

import itertools
import json

import numpy as np

from pymatgen.analysis.chemenv.coordination_environments.coordination_geometries import AllCoordinationGeometries
from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import (
    AbstractGeometry,
    LocalGeometryFinder,
)
from pymatgen.analysis.chemenv.utils.coordination_geometry_utils import Plane, collinear

if __name__ == "__main__":
    # Choose the geometry
    all_cg = AllCoordinationGeometries()
    while True:
        cg_symbol = input("Enter symbol of the geometry for which you want to get the explicit permutations : ")
        try:
            cg = all_cg[cg_symbol]
            break
        except LookupError:
            print("Wrong geometry, try again ...")
            continue

    # Check if the algorithm currently defined for this geometry corresponds to the explicit permutation algorithm
    for algo in cg.algorithms:
        if algo.algorithm_type != "SEPARATION_PLANE":
            raise ValueError("WRONG ALGORITHM !")

    new_algos = []

    idx = 1
    for sep_plane_algo in cg._algorithms:
        print(f"In {idx = }/{len(cg._algorithms)}")
        idx += 1
        if sep_plane_algo.algorithm_type != "SEPARATION_PLANE":
            raise ValueError("Should all be separation plane")

        perms_on_file = f"Permutations on file in this algorithm ({len(sep_plane_algo._permutations)}) "
        print(f"{perms_on_file}\n{sep_plane_algo._permutations}")
        permutations = sep_plane_algo.safe_separation_permutations(
            ordered_plane=sep_plane_algo.ordered_plane,
            ordered_point_groups=sep_plane_algo.ordered_point_groups,
        )

        sep_plane_algo._permutations = permutations

        print(f"Test permutations ({len(permutations)}):\n{permutations}")

        lgf = LocalGeometryFinder()
        lgf.setup_parameters(structure_refinement=lgf.STRUCTURE_REFINEMENT_NONE)
        lgf.setup_test_perfect_environment(
            cg_symbol,
            randomness=True,
            indices=range(cg.coordination_number),
            max_random_dist=0.05,
        )

        lgf.perfect_geometry = AbstractGeometry.from_cg(cg=cg)

        # Setting up the plane of separation
        local_plane = None
        found = False
        for n_points in range(
            sep_plane_algo.minimum_number_of_points,
            min(sep_plane_algo.maximum_number_of_points, 4) + 1,
        ):
            if found:
                break
            for ipoints in itertools.combinations(sep_plane_algo.plane_points, n_points):
                points_combination = [lgf.local_geometry.coords[ipoint] for ipoint in ipoints]
                if n_points == 2:
                    if collinear(
                        points_combination[0],
                        points_combination[1],
                        lgf.local_geometry.central_site,
                        tolerance=0.25,
                    ):
                        continue
                    local_plane = Plane.from_3points(
                        points_combination[0],
                        points_combination[1],
                        lgf.local_geometry.central_site,
                    )
                    found = True
                    break
                elif n_points == 3:
                    if collinear(
                        points_combination[0],
                        points_combination[1],
                        points_combination[2],
                        tolerance=0.25,
                    ):
                        continue
                    local_plane = Plane.from_3points(
                        points_combination[0],
                        points_combination[1],
                        points_combination[2],
                    )
                    found = True
                    break
                elif n_points > 3:
                    local_plane = Plane.from_npoints(points_combination, best_fit="least_square_distance")
                    found = True
                    break
                else:
                    raise ValueError("Wrong number of points to initialize separation plane")

        points_perfect = lgf.perfect_geometry.points_wocs_ctwocc()
        # Actual test of the permutations
        cgsm = lgf._cg_csm_separation_plane(
            coordination_geometry=cg,
            sep_plane=sep_plane_algo,
            local_plane=local_plane,
            plane_separations=[],
            dist_tolerances=[0.05, 0.1, 0.2, 0.3],
            testing=True,
            points_perfect=points_perfect,
        )

        print(cgsm)
        if cgsm[0] is None:
            print("IS NONE !")
            input()
            continue

        csms, perms, algos, sep_perms = cgsm[0], cgsm[1], cgsm[2], cgsm[3]

        print("Continuous symmetry measures")
        print(csms)
        csms_with_recorded_permutation: list[float] = []
        explicit_permutations = []
        for icsm, csm in enumerate(csms):
            found = False
            for csm2 in csms_with_recorded_permutation:
                if np.isclose(csm, csm2, rtol=0.0, atol=1e-6):
                    found = True
                    break
            if not found:
                print(perms[icsm], csm)
                csms_with_recorded_permutation.append(csm)
                explicit_permutations.append(sep_perms[icsm])

        print(perms_on_file)
        print(f"Permutations found ({len(explicit_permutations)}) : ")
        print(explicit_permutations)
        sep_plane_algo.explicit_permutations = explicit_permutations
        new_algos.append(sep_plane_algo)

    # Write update geometry file ?
    test = input('Save it ? ("y" to confirm)')
    if test == "y":
        cg._algorithms = new_algos
        cg_dict = cg.as_dict()
        with open(f"../coordination_geometries_files_new/{cg_symbol}.json", mode="w", encoding="utf-8") as file:
            json.dump(cg_dict, file)