-
Notifications
You must be signed in to change notification settings - Fork 838
/
explicit_permutations_plane_algorithm.py
153 lines (132 loc) · 5.82 KB
/
explicit_permutations_plane_algorithm.py
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
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
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
allcg = AllCoordinationGeometries()
while True:
cg_symbol = input("Enter symbol of the geometry for which you want to get the explicit permutations : ")
try:
cg = allcg[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 = []
ialgo = 1
for sep_plane_algo in cg._algorithms:
print(f"In {ialgo = :d}/{len(cg._algorithms):d}")
ialgo += 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):d}) "
print(perms_on_file)
print(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):d}) :")
print(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,
sepplane=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=1.0e-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):d}) : ")
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", "w") as f:
json.dump(cg_dict, f)