diff --git a/0.py b/0.py new file mode 100644 index 00000000..42201417 --- /dev/null +++ b/0.py @@ -0,0 +1,8 @@ +from pyxtal import pyxtal +from pyxtal.lattice import Lattice +struc = pyxtal(molecular=True) +sites = [{"4e": [0.77, 0.57, 0.53]}] +lat = Lattice.from_para(11.43, 6.49, 11.19, 90, 83.31, 90, ltype="monoclinic") +for i in range(200): + struc.from_random(3, 14, ["aspirin"], [4], lattice=lat, sites=sites) + print(i, struc.lattice) diff --git a/CTMTNA.cif b/CTMTNA.cif new file mode 100644 index 00000000..ad02bb9c --- /dev/null +++ b/CTMTNA.cif @@ -0,0 +1,101 @@ +smiles: O=N(=O)N1CN(CN(C1)N(=O)=O)N(=O)=O +# Refcode: CTMTNA + +####################################################################### +# +# Cambridge Crystallographic Data Centre +# CCDC +# +####################################################################### +# +# If this CIF has been generated from an entry in the Cambridge +# Structural Database, then it will include bibliographic, chemical, +# crystal, experimental, refinement or atomic coordinate data resulting +# from the CCDC's data processing and validation procedures. +# +####################################################################### + +data_CTMTNA +_chemical_name_common 1,3,5-trinitro-1,3,5-triazinane +_chemical_formula_moiety 'C3 H6 N6 O6' +_chemical_name_systematic Cyclotrimethylene-trinitramine +_symmetry_cell_setting orthorhombic +_symmetry_space_group_name_H-M 'P b c a' +_symmetry_Int_Tables_number 61 +_space_group_name_Hall '-P 2ac 2ab' +loop_ +_symmetry_equiv_pos_site_id +_symmetry_equiv_pos_as_xyz +1 x,y,z +2 1/2-x,-y,1/2+z +3 1/2+x,1/2-y,-z +4 -x,1/2+y,1/2-z +5 -x,-y,-z +6 1/2+x,y,1/2-z +7 1/2-x,1/2+y,z +8 x,1/2-y,1/2+z +_cell_length_a 13.182(2) +_cell_length_b 11.574(2) +_cell_length_c 10.709(2) +_cell_angle_alpha 90 +_cell_angle_beta 90 +_cell_angle_gamma 90 +_cell_volume 1633.86 +_chemical_melting_point '478 K' +_cell_formula_units_Z 8 +_chemical_properties_physical explosive +loop_ +_atom_site_label +_atom_site_type_symbol +_atom_site_fract_x +_atom_site_fract_y +_atom_site_fract_z +C1 C 0.18390 0.35780 0.44000 +C2 C 0.05030 0.24400 0.33950 +C3 C 0.14870 0.38130 0.21590 +H1 H 0.24010 0.29440 0.42260 +H2 H 0.20130 0.40680 0.52440 +H3 H -0.02610 0.21170 0.35400 +H4 H 0.10150 0.17280 0.31920 +H5 H 0.20520 0.31660 0.19170 +H6 H 0.14400 0.44320 0.14130 +N1 N 0.17610 0.43600 0.33300 +N2 N 0.08770 0.29960 0.45370 +N3 N 0.05360 0.32220 0.23290 +N4 N 0.22600 0.53790 0.33460 +N5 N 0.01550 0.35250 0.52940 +N6 N -0.03330 0.38750 0.20780 +O1 O 0.22700 0.59320 0.23880 +O2 O 0.26490 0.56870 0.43380 +O3 O -0.06930 0.31420 0.52620 +O4 O 0.04540 0.42710 0.59900 +O5 O -0.11210 0.35340 0.25040 +O6 O -0.02360 0.46910 0.13900 +loop_ +_geom_bond_atom_site_label_1 +_geom_bond_atom_site_label_2 +_geom_bond_site_symmetry_1 +_geom_bond_site_symmetry_2 +C1 H1 1_555 1_555 +C2 H3 1_555 1_555 +C3 H5 1_555 1_555 +H2 C1 1_555 1_555 +H4 C2 1_555 1_555 +H6 C3 1_555 1_555 +N1 C1 1_555 1_555 +N2 C1 1_555 1_555 +N3 C2 1_555 1_555 +N4 N1 1_555 1_555 +N5 N2 1_555 1_555 +N6 N3 1_555 1_555 +O1 N4 1_555 1_555 +O2 N4 1_555 1_555 +O3 N5 1_555 1_555 +O4 N5 1_555 1_555 +O5 N6 1_555 1_555 +O6 N6 1_555 1_555 +C2 N2 1_555 1_555 +C3 N1 1_555 1_555 +C3 N3 1_555 1_555 + +#END diff --git a/CTMTNA.xyz b/CTMTNA.xyz new file mode 100644 index 00000000..818744bb --- /dev/null +++ b/CTMTNA.xyz @@ -0,0 +1,23 @@ +21 +Cyclotrimethylene-trinitramine +C 2.42417 4.14118 4.71196 +C 0.66305 2.82406 3.63571 +C 1.96016 4.41317 2.31207 +H 3.16500 3.40739 4.52562 +H 2.65354 4.70830 5.61580 +H -0.34405 2.45022 3.79099 +H 1.33797 1.99999 3.41831 +H 2.70495 3.66433 2.05292 +H 1.89821 5.12960 1.51318 +N 2.32135 5.04626 3.56610 +N 1.15606 3.46757 4.85867 +N 0.70656 3.72914 2.49413 +N 2.97913 6.22565 3.58323 +N 0.20432 4.07984 5.66934 +N -0.43896 4.48493 2.22533 +O 2.99231 6.86570 2.55731 +O 3.49191 6.58213 4.64556 +O -0.91351 3.63655 5.63508 +O 0.59846 4.94326 6.41469 +O -1.47770 4.09025 2.68153 +O -0.31110 5.42936 1.48855 diff --git a/csp_xyz.py b/csp_xyz.py new file mode 100644 index 00000000..f6f83562 --- /dev/null +++ b/csp_xyz.py @@ -0,0 +1,184 @@ +""" +This is an example to perform CSP based on a reference crystal. +The structures with good matches will be output to *-matched.cif. +Supports molecule input as .xyz or .smi. +""" +from pyxtal import pyxtal +from pyxtal.optimize import WFS, DFS, QRS +from pyxtal.molecule import pyxtal_molecule, compare_mol_connectivity +import argparse +import os + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("-g", "--gen", dest="gen", type=int, default=1, + help="Number of generation, default: 1") + parser.add_argument("-p", "--pop", dest="pop", type=int, default=10, + help="Population size, default: 10") + parser.add_argument("-n", "--ncpu", dest="ncpu", type=int, default=1, + help="cpu number, default: 1") + parser.add_argument("-a", "--algo", dest="algo", default="WFS", + help="algorithm, default: WFS") + parser.add_argument("--mlp", dest="mlp", default="MACE", + help="MLP backend for xyz-only relaxation (MACE/ANI/UMA), default: MACE") + parser.add_argument("--preopt", dest="preopt", action="store_true", + help="preoptimize the lattice and rotation") + parser.add_argument("--parameters", dest="parameters", default="parameters.xml", + help="forcefield parameter xml file, default: parameters.xml") + parser.add_argument("--ffstyle", dest="ffstyle", default="gaff", + help="forcefield style, default: gaff") + + # New: molecule + reference inputs + parser.add_argument("--mol", dest="mol", default="CTMTNA.xyz", + help="Molecule input (.xyz or .smi), default: CTMTNA.xyz") + parser.add_argument("--smiles", dest="smiles", default=None, + help="Optional SMILES for force-field mode; not needed for .xyz MLP-only runs") + parser.add_argument("--xyz-only", dest="xyz_only", action="store_true", + help="Force MLP-only mode (selected automatically for .xyz without --smiles)") + parser.add_argument("--nconf", dest="nconf", type=int, default=1, + help="Number of conformers generated from --smiles, default: 1") + parser.add_argument("--niter-conf", dest="niter_conf", type=int, default=5, + help="Embedding iterations for conformer generation, default: 5") + parser.add_argument("--conf-tol", dest="conf_tol", type=float, default=0.5, + help="RMSD tolerance for unique conformers, default: 0.5") + parser.add_argument("--uff", dest="use_uff", action="store_true", + help="Use UFF (instead of MMFF) in conformer generation") + parser.add_argument("--seed", dest="seed", default=None, + help="Optional reference crystal CIF used only for match checking") + parser.add_argument("--wdir", dest="wdir", default="ctmtna-simple", + help="Working directory, default: ctmtna-simple") + parser.add_argument("--sg", dest="sg", type=int, nargs="+", default=[61], + help="Space group list, default: 61") + # add CLI arg + parser.add_argument( + "--active-sites", + dest="active_sites", + type=int, + nargs=3, + metavar=("DONOR", "ACCEPTOR", "H"), + default=None, + help="Optional atom indices for active sites, e.g. --active-sites 11 12 20", + ) + + options = parser.parse_args() + + # build active_sites only if provided + active_sites = None + if options.active_sites is not None: + d, a, h = options.active_sites + active_sites = [[d], [a], [h]] + + mol_paths = [x.strip() for x in options.mol.split(",") if len(x.strip()) > 0] + if len(mol_paths) == 0: + raise ValueError("Please provide at least one molecule path via --mol") + mol_ext = os.path.splitext(mol_paths[0])[1].lower() + for mpath in mol_paths: + ext = os.path.splitext(mpath)[1].lower() + if ext != ".smi" and not os.path.exists(mpath): + raise FileNotFoundError(f"Cannot find molecule file: {mpath}") + if options.seed is not None and not os.path.exists(options.seed): + raise FileNotFoundError(f"Cannot find reference CIF: {options.seed}") + + xyz_only = options.xyz_only or ( + options.smiles is None + and all(os.path.splitext(path)[1].lower() == ".xyz" for path in mol_paths) + ) + + # Optimizer currently requires a valid SMILES string internally + if xyz_only: + if options.algo != "WFS": + raise ValueError("xyz-only mode currently supports --algo WFS only") + smiles_opt = options.smiles if options.smiles is not None else "xyz_only" + elif options.smiles is not None: + smiles_opt = options.smiles + elif mol_ext == ".smi": + smiles_opt = os.path.splitext(os.path.basename(mol_paths[0]))[0] + else: + raise ValueError( + "For .xyz input, please also provide --smiles '' " + "so optimizer torsion/FF initialization can proceed." + ) + + # Build molecule pool: + # - single geometry from --mol + # - or multiple conformers generated from --smiles + if xyz_only: + mol_objs = [pyxtal_molecule(mpath, active_sites=active_sites) for mpath in mol_paths] + molecules_for_opt = [mol_objs] + seed_molecules = [mol_objs[0]] + elif options.nconf > 1: + if options.smiles is None: + raise ValueError("Please provide --smiles when --nconf > 1") + confs = pyxtal_molecule.get_conformers_from_smiles( + options.smiles, + N_iter=options.niter_conf, + N_conf=options.nconf, + tol=options.conf_tol, + use_uff=options.use_uff, + ) + if len(confs) == 0: + raise RuntimeError("No conformers were generated from --smiles") + if active_sites is not None: + for m in confs: + m.active_sites = active_sites + molecules_for_opt = [confs] + seed_molecules = [confs[0]] + else: + m1 = pyxtal_molecule(mol_paths[0], active_sites=active_sites) + # For xyz input, verify topology matches the SMILES-based template used + # by optimizer internals (torsion/FF typing). + if mol_ext == ".xyz": + m_smi = pyxtal_molecule(smiles_opt + ".smi") + ok, _ = compare_mol_connectivity(m1.mol, m_smi.mol, ignore_HH=True) + if not ok: + raise ValueError( + "Input .xyz connectivity does not match --smiles template. " + "HTOCSP force-field/torsion setup is SMILES-driven, so xyz-only " + "input is not currently sufficient for robust optimization. " + "Use --mol '.smi' (or --nconf from SMILES), or provide " + "an xyz generated with the same atom/bond topology convention." + ) + molecules_for_opt = [[m1]] + seed_molecules = [m1] + + # A reference CIF is optional. It is used for match checking, not for + # defining the molecule; molecular geometry comes from --mol. + pmg = None + if options.seed is not None: + xtal = pyxtal(molecular=True) + xtal.from_seed(options.seed, molecules=seed_molecules) + pmg = xtal.to_pymatgen() + + # Sampling + fun = globals().get(options.algo) + if fun is None: + raise ValueError(f"Unknown algorithm: {options.algo}. Choose from WFS, DFS, QRS.") + + run_kwargs = dict( + tag="csp_run", + N_gen=options.gen, + N_pop=options.pop, + N_cpu=options.ncpu, + mlp=options.mlp, + skip_mlp=False if xyz_only else True, + ff_style=options.ffstyle, + ff_parameters=options.parameters, + molecules=molecules_for_opt, + xyz_only=xyz_only, + ) + if options.algo != "QRS": + run_kwargs["pre_opt"] = options.preopt + if xyz_only and options.algo == "WFS": + run_kwargs["fracs"] = [1.0, 0.0] + + go = fun( + smiles_opt, + options.wdir, + options.sg, + **run_kwargs, + ) + + go.run(ref_pmg=pmg) + if pmg is not None: + go.print_matches(header="Ref_match") + go.plot_results() diff --git a/pyxtal/optimize/WFS.py b/pyxtal/optimize/WFS.py index e38d6922..712c749d 100644 --- a/pyxtal/optimize/WFS.py +++ b/pyxtal/optimize/WFS.py @@ -92,6 +92,7 @@ def __init__( use_mpi: bool = False, pre_opt: bool = False, check: bool = True, + xyz_only: bool = False, ): if isinstance(random_state, Generator): self.random_state = random_state.spawn(1)[0] @@ -142,6 +143,7 @@ def __init__( check_stable, use_mpi, pre_opt, + xyz_only=xyz_only, ) if fracs is None: diff --git a/pyxtal/optimize/base.py b/pyxtal/optimize/base.py index 74c3d8a9..bc8d817f 100644 --- a/pyxtal/optimize/base.py +++ b/pyxtal/optimize/base.py @@ -176,6 +176,7 @@ def __init__( use_mpi: bool = False, pre_opt: bool = False, N_min_matches: int = 10, + xyz_only: bool = False, ): self.ncpu = N_cpu @@ -198,6 +199,7 @@ def __init__( # Molecular information self.smile = smiles self.smiles = self.smile.split(".") # list + self.xyz_only = xyz_only self.torsions = torsions self.molecules = molecules self.block = block @@ -205,9 +207,15 @@ def __init__( self.composition = [ 1] * len(self.smiles) if composition is None else composition self.N_torsion = 0 - for smi, comp in zip(self.smiles, self.composition): - self.N_torsion += len(find_rotor_from_smile(smi) - ) * int(max([comp, 1])) + if self.xyz_only and self.molecules is not None: + for pool, comp in zip(self.molecules, self.composition): + molecule = pool[0] if isinstance(pool, (list, tuple)) else pool + torsionlist = getattr(molecule, "torsionlist", None) or [] + self.N_torsion += len(torsionlist) * int(max([comp, 1])) + else: + for smi, comp in zip(self.smiles, self.composition): + self.N_torsion += len(find_rotor_from_smile(smi) + ) * int(max([comp, 1])) # Crystal information self.pre_opt = pre_opt @@ -259,7 +267,15 @@ def __init__( filename=self.log_file, level=logging.INFO) self.logging = logging - if info is not None: + if self.xyz_only: + if self.skip_mlp: + raise ValueError("xyz_only requires skip_mlp=False") + if self.molecules is None: + raise ValueError("xyz_only requires pre-built molecular geometries") + self.atom_info = {} + self.parameters = None + self.ff_opt = False + elif info is not None: self.atom_info = info self.parameters = None self.ff_opt = False @@ -461,8 +477,9 @@ def run(self, ref_pmg=None, ref_pxrd=None, max_rmsd=None): results = self._run(pool) except (EOFError, OSError) as e: print(f"Error in running the optimizer: {e}") - pool.terminate() - pool.join() + if pool is not None: + pool.terminate() + pool.join() return None if self.rank == 0: @@ -796,11 +813,22 @@ def _apply_gaussian(self, reps, engs, h1=0.1, h2=0.1, w1=0.2, w2=3): g1 = h1 * np.exp(-0.5 * diff1) # cell # Torsion g2 = 0 - if len(tor1) > 0: + # Relaxation may discover higher symmetry, so a Z'>1 + # candidate and a saved Z'=1 structure can have different + # numbers (or layouts) of molecular-site records. Their + # torsion vectors are not directly comparable. Keep the + # lattice Gaussian, but only compare torsions when the + # representation layouts match. + same_layout = ( + len(ref) == len(rep) + and all(len(ref[j]) == len(rep[j]) + for j in range(1, len(rep))) + ) + if len(tor1) > 0 and same_layout: tor2 = np.zeros(self.N_torsion) count = 0 - for j in range(1, len(rep)): - if len(rep[j]) > N_id: # for Cl- + for j in range(1, len(ref)): + if len(ref[j]) > N_id: # for Cl- tor2[count: count + len(ref[j]) - N_id - 1] = ref[j][N_id:-1] count += len(ref[j]) - N_id @@ -900,6 +928,7 @@ def _get_local_optimization_args(self): self.opt_lat, getattr(self, 'delta_length', 1.0), getattr(self, 'delta_angle', 15.0), + self.xyz_only, ] return args diff --git a/pyxtal/optimize/common.py b/pyxtal/optimize/common.py index cfffa1db..2f07365d 100644 --- a/pyxtal/optimize/common.py +++ b/pyxtal/optimize/common.py @@ -553,6 +553,7 @@ def optimizer( skip_mlp = False, output_mlp = True, pre_opt = False, + xyz_only = False, ): """ Structural relaxation for each individual pyxtal structure. @@ -577,6 +578,33 @@ def optimizer( if pre_opt: struc.optimize_lattice_and_rotation() + if xyz_only: + cwd = os.getcwd() + t0 = time() + os.makedirs(workdir, exist_ok=True) + os.chdir(workdir) + try: + s = ASE_relax( + struc, + mlp, + opt_lat=opt_lat, + step=200 if opt_lat else 50, + fmax=0.1, + logfile="ase.log", + ) + if s is None: + return None + eng = s.get_potential_energy() + stress = max(abs(s.get_stress())) / units.GPa + if stress > 30.0: + return None + + xtal = pyxtal(molecular=True) + xtal.from_seed(ase2pymatgen(s), molecules=struc.molecules) + return {"xtal": xtal, "energy": eng, "time": time() - t0} + finally: + os.chdir(cwd) + if calculators is None: calculators = ["CHARMM"] @@ -654,7 +682,7 @@ def optimizer( struc.energy < 9999 and struc.lattice.is_valid_matrix() # and struc.check_distance() - and 0.25 < struc.get_density() < 3.0 + and 1.25 < struc.get_density() < 3.0 ): s = struc.to_ase() step = 50 if mlp in ['MACE', 'ANI'] else 25 @@ -665,14 +693,13 @@ def optimizer( t = time() - t0 if t > max_time: - try: - print("!!!Long time in ani calculation", t) - print(struc.get_1D_representation().to_string()) - struc.optimize_lattice() - except: - print("Trouble in optLat") - return None - elif stress < stress_tol: + # The relaxation has already completed successfully. Do not + # discard it merely because it exceeded this advisory target; + # the worker-level SIGALRM enforces the actual timeout. + print(f"MLP relaxation exceeded advisory time: {t:.1f} s " + f"(target {max_time:.1f} s)") + + if stress < stress_tol: results = {} if output_mlp: xtal = pyxtal(molecular=True) @@ -732,6 +759,7 @@ def optimizer_par( opt_lat=True, delta_length=1.0, delta_angle=15.0, + xyz_only=False, ): """ A routine used for parallel structure optimization @@ -777,6 +805,7 @@ def optimizer_par( opt_lat=opt_lat, delta_length=delta_length, delta_angle=delta_angle, + xyz_only=xyz_only, label=labels[i] if labels is not None else None, ) results.append((id, xtal, match, stable)) @@ -814,6 +843,7 @@ def optimizer_single( opt_lat=True, delta_length=1.0, delta_angle=15.0, + xyz_only=False, label=None, ): """ @@ -859,7 +889,7 @@ def optimizer_single( else: res = optimizer(xtal, atom_info, workdir, job_tag, opt_lat, mlp=mlp, skip_mlp=skip_mlp, output_mlp=output_mlp, - pre_opt=pre_opt) + pre_opt=pre_opt, xyz_only=xyz_only) match = False # used for matching with reference stable = True # used for tagging if the structure is stable @@ -917,8 +947,11 @@ def optimizer_single( tag += 'Stable' else: tag += 'Shallow' - rep = xtal.get_1D_representation() - strs = rep.to_string(None, eng / N, tag) + try: + rep = xtal.get_1D_representation() + strs = rep.to_string(None, eng / N, tag) + except Exception: + strs = f"{tag:8s} E={eng / N:12.3f}" # 3. Check match w.r.t the reference if ref_pmg is not None: @@ -933,7 +966,7 @@ def optimizer_single( # Further refine the structure match = True str1 = f"Match {rmsd[0]:6.2f} {rmsd[1]:6.2f} {eng / N:12.3f} " - if not skip_mlp: + if not skip_mlp and not xyz_only: xtal, eng1 = refine_struc(xtal, smiles, ASE_relax, mlp) str1 += f"Full Relax -> {eng1 / N:12.3f}" eng = eng1