fix issue for rank and match

scripts/match.py

@@ -0,0 +1,127 @@
+
+import pymatgen as mg
+import numpy as np
+from pyxtal.util import parse_cif
+from optparse import OptionParser
+import pymatgen.analysis.structure_matcher as sm
+from pyxtal import pyxtal
+from pyxtal.XRD import Similarity
+from pyxtal.optimize.base import GlobalOptimize
+import warnings
+warnings.filterwarnings("ignore")
+
+parser = OptionParser()
+parser.add_option("-f", dest="cif", default="WFS-gaff.cif", help="input cif file")
+parser.add_option("-r", dest="ref", help="reference")
+parser.add_option("-o", dest="out", default="Matched.cif", help="output")
+parser.add_option("--emin", dest="emin", type=float, default=0,
+ help="minimum energy, default 0")
+parser.add_option("--emax", dest="emax", type=float, default=100,
+ help="maximum energy, default 100")
+parser.add_option("--early_stop", dest="early", action="store_true", default=False,
+ help="stop when the first match is found")
+parser.add_option("--XRD", dest="xrd",
+ action="store_true", default=False,
+ help="Compare XRD")
+parser.add_option("--smin", dest="smin", type=float, default=0.80,
+ help="min similarity for XRD")
+
+(options, args) = parser.parse_args()
+matcher = sm.StructureMatcher(ltol=0.3, stol=0.4, angle_tol=5.0)
+
+with open(options.cif, 'r') as f:
+ lines = f.readlines()
+ smiles = []
+ for l in lines:
+ if 'smile' in l:
+ smile_str = l.split(':')[1].strip()
+ smiles = [s + '.smi' for s in smile_str.split('.')]
+ break
+print(smiles)
+with open(options.out, 'w') as f: f.write(f'smiles: {smile_str}\n')
+
+if options.ref is None:
+ raise ValueError("Reference structure is required.")
+else:
+ pmg_ref = mg.core.Structure.from_file(options.ref)
+ xtal = pyxtal(molecular=True)
+ xtal.from_seed(pmg_ref, molecules = smiles)
+ print(f"Reference Structure loaded from {options.ref} {pmg_ref.density:.3f}")
+ pmg_ref.remove_species("H")
+ print(xtal)
+
+ if options.xrd:
+ thetas = [0, 35.0]
+ xrd = xtal.get_XRD(thetas=thetas)
+ p_ref = xrd.get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
+
+cifs, engs = parse_cif(options.cif, eng=True)
+print("Total Number of Structures:", len(cifs))
+engs = np.array(engs)
+ids = np.argsort(engs)
+cifs = [cifs[id] for id in ids]
+engs = engs[ids]
+print(f"Min energy in eV: {engs.min()+options.emin/96.485:.4f} {engs.min()+options.emax/96.485:.4f}")
+engs_norm = engs - engs.min() # Normalize energies to the lowest one
+engs_norm *= 96.485
+
+# Find the id of energy that is between [options.emin, options.emax]
+n1 = np.searchsorted(engs_norm, options.emin, side='left')
+n2 = np.searchsorted(engs_norm, options.emax, side='right')
+engs = engs[n1:n2]
+engs_norm = engs_norm[n1:n2]
+cifs = [cifs[id] for id in range(n1, n2)]
+ids = ids[n1:n2]
+count = 0
+xtal = pyxtal(molecular=True)
+for id, cif in enumerate(cifs):
+ pmg = mg.core.Structure.from_str(cif, fmt='cif')
+ try:
+ xtal.from_seed(pmg, molecules = smiles)
+ #xtal.energy = engs_norm[id]
+ spg = xtal.group.number
+ den = xtal.get_density()
+ raw_eng = engs[id]
+ norm_eng = engs_norm[id]
+ match = False
+ strs = f"Struc {ids[id]:6d}: {spg:3d} {raw_eng:.3f} kJ/mol, {den:.3f} g/cm^3, {norm_eng:.3f}"
+ if options.xrd:
+ p1 = xtal.get_XRD(thetas=thetas).get_profile(res=0.15, user_kwargs={"FWHM": 0.25})
+ sim = Similarity(p1, p_ref, x_range=thetas).value
+ if sim > options.smin: match = True
+ strs += f'{sim:12.3f} in PXRD similarity'
+ else:
+ pmg.remove_species("H")
+ if abs(pmg.density-pmg_ref.density) <= 0.35:
+ strs += '****'
+ if matcher.fit(pmg, pmg_ref):
+ match = True
+
+ pmg_gen = xtal.to_pymatgen()
+ pmg_gen.remove_species("H")
+
+ ref_copy = pmg_ref.copy()
+ ref_copy.remove_species("H")
+ print(f"Generated structure sites: {len(pmg_gen.sites)}, Reference sites: {len(ref_copy.sites)}")
+ print(f"Species (gen): {[str(sp) for sp in pmg_gen.species]}")
+ print(f"Species (ref): {[str(sp) for sp in ref_copy.species]}")
+ rmsd = matcher.get_rms_dist( ref_copy, pmg_gen)
+ if rmsd is not None:
+ rms_lat, rms_cart = rmsd
+ print(f"RMSD – Lattice: {rms_lat:.3f} Å, Cartesian: {rms_cart:.3f} Å")
+ else:
+ print("RMSD calculation failed: structures not comparable")
+
+ except:
+ continue
+
+ if match:
+ count += 1
+ strs += '+++++++++++'
+ label = f"{count}-d{den:.3f}-spg{spg}-e{norm_eng:.3f}"
+ if options.xrd: label += f"-s{sim:.3f}"
+ with open(options.out, 'a+') as f: f.writelines(xtal.to_file(header=label))
+ if options.early: break
+ print(strs)
+print(f"Found {count} matches")
+

scripts/rank.py

@@ -0,0 +1,148 @@
+import pymatgen as mg
+import numpy as np
+from pyxtal import pyxtal
+from pyxtal.util import parse_cif
+from optparse import OptionParser
+import pymatgen.analysis.structure_matcher as sm
+
+import warnings
+warnings.filterwarnings("ignore")
+
+def new_struc(xtal, xtals, max_num=100):
+ """
+ check if this is a new structure
+
+ Args:
+ xtal: input structure
+ xtals: list of reference structures
+
+ Return:
+ `None` or the id of matched structure
+ """
+ if len(xtals) > max_num:
+ start = len(xtals) - max_num
+ else:
+ start = 0
+ sg1 = xtal.group.number
+ pmg_s1 = xtal.to_pymatgen()
+ pmg_s1.remove_species("H")
+ vol1 = pmg_s1.lattice.volume
+
+ for xtal2 in xtals[start:]:
+ sg2 = xtal2.group.number
+ if sg1 == sg2:
+ pmg_s2 = xtal2.to_pymatgen()
+ vol2 = pmg_s2.lattice.volume
+ if abs(vol1-vol2)/vol1<5e-2:
+ pmg_s2.remove_species("H")
+ if sm.StructureMatcher().fit(pmg_s1, pmg_s2):
+ return False
+ return True
+
+
+parser = OptionParser()
+parser.add_option("-f", "--cif", dest="cif", default="WFS-gaff.cif",
+ help="cif file name, optional")
+parser.add_option("-r", "--rank", dest="rank", default='energy',
+ help="ranking criteria: default is energy")
+parser.add_option("-s", "--n1", dest="n1", type=int, default=0,
+ help="starting id, optional")
+parser.add_option("-e", "--n2", dest="n2", type=int, default=-1,
+ help="ennding id, optional")
+parser.add_option("-c", "--cut", dest="cut", type=int,
+ help="cutoff number, optional")
+parser.add_option("--dmax", dest="dmax", type=float, default=10.0,
+ help="maximum density in g/cm^3, optional")
+
+(options, args) = parser.parse_args()
+rank = options.rank
+n1 = options.n1
+n2 = options.n2
+output1 = 'Ranked.cif'
+
+"""
+Read the smile from the following contents
+-------Global Crystal Structure Prediction------
+smile : CC(=O)OC1=CC=CC=C1C(=O)O
+"""
+
+with open(options.cif, 'r') as f:
+ lines = f.readlines()
+ smiles = []
+ for l in lines:
+ if 'smile' in l:
+ smile_str = l.split(':')[1].strip()
+ smiles = [smile_str + '.smi']
+ break
+print(smiles)
+
+cifs, engs = parse_cif(options.cif, eng=True)
+print("Total Number of Structures:", len(cifs))
+if options.cut is None:
+ cut = len(cifs)
+else:
+ cut = options.cut
+
+if options.rank == 'energy':
+ engs = np.array(engs)
+ ids = np.argsort(engs)
+else:
+ cifs, sims = parse_cif(options.cif, sim=True)
+ sims = np.array(sims)
+ ids = np.argsort(-1*sims)
+ sims = [sims[id] for id in ids]
+
+cifs = [cifs[id] for id in ids]
+engs = engs[ids]
+eng0 = engs.min()
+
+if n2 == -1:
+ cifs = cifs[n1:]
+ engs = engs[n1:]
+ ids = ids[n1:]
+else:
+ n2 = min(n2, len(cifs))
+ cifs = cifs[n1:n2]
+ engs = engs[n1:n2]
+ ids = ids[n1:n2]
+print("Index", n1, n2, cut, len(cifs))
+with open(output1, 'w') as f: f.write(l)
+
+xtals = []
+count = 0
+with open(output1, 'a+') as f:
+ for id, cif in enumerate(cifs):
+ pmg = mg.core.Structure.from_str(cif, fmt='cif')
+ try:
+ xtal = pyxtal(molecular=True)
+ xtal.from_seed(pmg, molecules = smiles)
+ xtal.energy = engs[id]
+ if new_struc(xtal, xtals, 100):
+ xtals.append(xtal)
+ spg = xtal.group.number
+ den = xtal.get_density()
+ eng = engs[id]
+ label = f"{count}-d{den:.3f}-spg{spg}-e{eng:.3f}"
+ if den < options.dmax:
+ cif_lines = xtal.to_file(header=label).splitlines()
+ # Extract energy from label
+ energy_from_label = float(label.split('-e')[1])
+ energy_line = f"#Energy: {energy_from_label:.3f} eV/cell\n"
+ output_lines = []
+ energy_replaced = False
+ for line in cif_lines:
+ stripped = line.strip()
+ # Replace any existing energy line
+ if stripped.startswith("#Energy:"):
+ output_lines.append(energy_line)
+ energy_replaced = True
+ continue
+ output_lines.append(line + '\n')
+ f.writelines(output_lines)
+ print(f"{ids[id]:6d} {label} {(eng - eng0)*96.485:6.2f}")
+ count += 1
+ if count == cut:
+ print("Stop", count)
+ break
+ except:
+ print("Problem in reading")