Skip to content

Commit 928ff76

Browse files
author
Musiha Mukta
committed
fix issue+add rmsd
1 parent 11157ec commit 928ff76

File tree

1 file changed

+127
-0
lines changed

1 file changed

+127
-0
lines changed

scripts/match.py

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
2+
import pymatgen as mg
3+
import numpy as np
4+
from pyxtal.util import parse_cif
5+
from optparse import OptionParser
6+
import pymatgen.analysis.structure_matcher as sm
7+
from pyxtal import pyxtal
8+
from pyxtal.XRD import Similarity
9+
from pyxtal.optimize.base import GlobalOptimize
10+
import warnings
11+
warnings.filterwarnings("ignore")
12+
13+
parser = OptionParser()
14+
parser.add_option("-f", dest="cif", default="WFS-gaff.cif", help="input cif file")
15+
parser.add_option("-r", dest="ref", help="reference")
16+
parser.add_option("-o", dest="out", default="Matched.cif", help="output")
17+
parser.add_option("--emin", dest="emin", type=float, default=0,
18+
help="minimum energy, default 0")
19+
parser.add_option("--emax", dest="emax", type=float, default=100,
20+
help="maximum energy, default 100")
21+
parser.add_option("--early_stop", dest="early", action="store_true", default=False,
22+
help="stop when the first match is found")
23+
parser.add_option("--XRD", dest="xrd",
24+
action="store_true", default=False,
25+
help="Compare XRD")
26+
parser.add_option("--smin", dest="smin", type=float, default=0.80,
27+
help="min similarity for XRD")
28+
29+
(options, args) = parser.parse_args()
30+
matcher = sm.StructureMatcher(ltol=0.3, stol=0.4, angle_tol=5.0)
31+
32+
with open(options.cif, 'r') as f:
33+
lines = f.readlines()
34+
smiles = []
35+
for l in lines:
36+
if 'smile' in l:
37+
smile_str = l.split(':')[1].strip()
38+
smiles = [s + '.smi' for s in smile_str.split('.')]
39+
break
40+
print(smiles)
41+
with open(options.out, 'w') as f: f.write(f'smiles: {smile_str}\n')
42+
43+
if options.ref is None:
44+
raise ValueError("Reference structure is required.")
45+
else:
46+
pmg_ref = mg.core.Structure.from_file(options.ref)
47+
xtal = pyxtal(molecular=True)
48+
xtal.from_seed(pmg_ref, molecules = smiles)
49+
print(f"Reference Structure loaded from {options.ref} {pmg_ref.density:.3f}")
50+
pmg_ref.remove_species("H")
51+
print(xtal)
52+
53+
if options.xrd:
54+
thetas = [0, 35.0]
55+
xrd = xtal.get_XRD(thetas=thetas)
56+
p_ref = xrd.get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
57+
58+
cifs, engs = parse_cif(options.cif, eng=True)
59+
print("Total Number of Structures:", len(cifs))
60+
engs = np.array(engs)
61+
ids = np.argsort(engs)
62+
cifs = [cifs[id] for id in ids]
63+
engs = engs[ids]
64+
print(f"Min energy in eV: {engs.min()+options.emin/96.485:.4f} {engs.min()+options.emax/96.485:.4f}")
65+
engs_norm = engs - engs.min() # Normalize energies to the lowest one
66+
engs_norm *= 96.485
67+
68+
# Find the id of energy that is between [options.emin, options.emax]
69+
n1 = np.searchsorted(engs_norm, options.emin, side='left')
70+
n2 = np.searchsorted(engs_norm, options.emax, side='right')
71+
engs = engs[n1:n2]
72+
engs_norm = engs_norm[n1:n2]
73+
cifs = [cifs[id] for id in range(n1, n2)]
74+
ids = ids[n1:n2]
75+
count = 0
76+
xtal = pyxtal(molecular=True)
77+
for id, cif in enumerate(cifs):
78+
pmg = mg.core.Structure.from_str(cif, fmt='cif')
79+
try:
80+
xtal.from_seed(pmg, molecules = smiles)
81+
#xtal.energy = engs_norm[id]
82+
spg = xtal.group.number
83+
den = xtal.get_density()
84+
raw_eng = engs[id]
85+
norm_eng = engs_norm[id]
86+
match = False
87+
strs = f"Struc {ids[id]:6d}: {spg:3d} {raw_eng:.3f} kJ/mol, {den:.3f} g/cm^3, {norm_eng:.3f}"
88+
if options.xrd:
89+
p1 = xtal.get_XRD(thetas=thetas).get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
90+
sim = Similarity(p1, p_ref, x_range=thetas).value
91+
if sim > options.smin: match = True
92+
strs += f'{sim:12.3f} in PXRD similarity'
93+
else:
94+
pmg.remove_species("H")
95+
if abs(pmg.density-pmg_ref.density) <= 0.35:
96+
strs += '****'
97+
if matcher.fit(pmg, pmg_ref):
98+
match = True
99+
100+
pmg_gen = xtal.to_pymatgen()
101+
pmg_gen.remove_species("H")
102+
103+
ref_copy = pmg_ref.copy()
104+
ref_copy.remove_species("H")
105+
print(f"Generated structure sites: {len(pmg_gen.sites)}, Reference sites: {len(ref_copy.sites)}")
106+
print(f"Species (gen): {[str(sp) for sp in pmg_gen.species]}")
107+
print(f"Species (ref): {[str(sp) for sp in ref_copy.species]}")
108+
rmsd = matcher.get_rms_dist( ref_copy, pmg_gen)
109+
if rmsd is not None:
110+
rms_lat, rms_cart = rmsd
111+
print(f"RMSD – Lattice: {rms_lat:.3f} Å, Cartesian: {rms_cart:.3f} Å")
112+
else:
113+
print("RMSD calculation failed: structures not comparable")
114+
115+
except:
116+
continue
117+
118+
if match:
119+
count += 1
120+
strs += '+++++++++++'
121+
label = f"{count}-d{den:.3f}-spg{spg}-e{norm_eng:.3f}"
122+
if options.xrd: label += f"-s{sim:.3f}"
123+
with open(options.out, 'a+') as f: f.writelines(xtal.to_file(header=label))
124+
if options.early: break
125+
print(strs)
126+
print(f"Found {count} matches")
127+

0 commit comments

Comments
 (0)