From 839ab66c0d0a36c76fec4a1308413b1537c32451 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Wed, 29 May 2024 10:14:47 -0400 Subject: [PATCH] add db option and update documentation on substitute_1_2 --- README.md | 37 +++++++++++++++++- doc/Usage.rst | 8 ++-- pyxtal/db.py | 83 ++++++++++++++++++++++++++++------------ pyxtal/interface/vasp.py | 27 +++++++------ 4 files changed, 112 insertions(+), 43 deletions(-) diff --git a/README.md b/README.md index f430dc97..d4ea2d0e 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ ## Introduction -PyXtal is an open source Python package which was initiated by [Qiang Zhu](http://qzhu2017.github.io) and Scott Fredericks at department of Physics and Astronomy, University of Nevada Las Vegas. The goal of PyXtal is to develop a fundamental library to allow one to design the material structure with a certain symmetry constraint. These structures can exported to various structural formats for further study. See the [documentation](https://pyxtal.readthedocs.io/en/latest/) for more information. +PyXtal is an open source Python package which was initiated by [Qiang Zhu](http://qzhu2017.github.io) and Scott Fredericks. The goal of PyXtal is to develop a fundamental library to allow one to design the material structure with a certain symmetry constraint. These structures can exported to various structural formats for further study. See the [documentation](https://pyxtal.readthedocs.io/en/latest/) for more information. To contribute to this project, please check [How to contribute?](#how-to-contribute). @@ -58,6 +58,8 @@ pip install --upgrade git+https://github.com/qzhu2017/PyXtal.git@master ## Citation +### General PyXtal + Fredericks S, Parrish K, Sayre D, Zhu Q\*(2020) [PyXtal: a Python Library for Crystal Structure Generation and Symmetry Analysis](https://www.sciencedirect.com/science/article/pii/S0010465520304057). \[[arXiv link](https://arxiv.org/pdf/1911.11123.pdf)\] @@ -76,6 +78,39 @@ author = "Scott Fredericks and Kevin Parrish and Dean Sayre and Qiang Zhu", } ``` +### Symmetry relation +Zhu Q, Kang B, Parrish K (2022). +[Symmetry Relation Database and Its Application to Ferroelectric Materials Discovery](https://link.springer.com/article/10.1557/s43579-022-00268-4) + +```bib +@article{zhu2022symmetry, + title="Symmetry relation database and its application to ferroelectric materials discovery", + author="Zhu, Qiang and Kang, Byungkyun and Parrish, Kevin", + journal="MRS Communications", + volume="12", + number="5", + pages="686--691", + year="2022", + doi="https://doi.org/10.1557/s43579-022-00268-4", +``` + +### Organic crystal packing motif +Zhu Q, Tang W-L, Hattori S. (2022). +[Quantification of Crystal Packing Similarity from Spherical Harmonic Transform](https://pubs.acs.org/doi/10.1021/acs.cgd.2c00933) + +```bib +@article{zhu2022quantification, + title="Quantification of Crystal Packing Similarity from Spherical Harmonic Transform", + author="Zhu, Qiang and Tang, Weilun and Hattori, Shinnosuke", + journal="Crystal Growth \& Design", + volume="22", + number="12", + pages="7308--7316", + year="2022", + doi="https://doi.org/10.1021/acs.cgd.2c00933", +} +``` + ## How to contribute? This is an open-source project. Its growth depends on the community. To contribute to PyXtal, you don't necessarily have to write the code. Any contributions from the following list will be helpful. diff --git a/doc/Usage.rst b/doc/Usage.rst index 59ef306e..418b9765 100644 --- a/doc/Usage.rst +++ b/doc/Usage.rst @@ -594,9 +594,11 @@ Symmetry relation has been playing an important role in crystallography. PyXtal Chemical Substitution ---------------------- -In many cases, the crystal structures of mutlicompnent systems are strongly related to the structure of simple systems. For instance, the 1: 1 ratio boron nitrides, as an isoelectronic analogue to carbon, exihibit structural behaviors that are very similar to elemental carbon allotropes. Similarly, many of the known AlPO4 polymorphs are related to SiO2. Inspired by these known correlation, PyXtal offers the `substitue_1_2 `_ function to derive the BC compounds from A via subgroup relation (e.g., from C to BN or from SiO2 to AlPO4). The key idea is to split A's wyckoff site to B and C according to the BC composition constraints. Unlike the random substituion, the wyckoff position splitting strictly follow the group-subgroup relation. As such, the resulting compounds will retain a high space group symmetry from the parental structure. Below, we illustrate this function via a few examples. +In many cases, the crystal structures of mutlicompnent systems are strongly related to the structure of simple systems. For instance, the 1: 1 ratio boron nitrides, as an isoelectronic analogue to carbon, exihibit very similar structural behaviors as compared to elemental carbon allotropes. Similarly, many of the known AlPO4 polymorphs are related to SiO2. -Below is an example to make a BN compound from the given carbon allotrope diamond. +Inspired by these known correlation, PyXtal offers the `substitue_1_2 `_ function to derive the BC compounds from A via subgroup relation (e.g., from C to BN or from SiO2 to AlPO4). The key idea is to split A's Wyckoff sites to B and C according to the BC composition constraints. Unlike the random substitution, the Wyckoff position splitting strictly follows the group-subgroup relation. As such, the resulting compound retains a high space group symmetry from the parental structure. Below, we illustrate this function via a few examples. + +Below is a script to make a 1:1 ratio BN compound from the diamond carbon allotrope. .. code-block:: Python @@ -615,7 +617,7 @@ Below is an example to make a BN compound from the given carbon allotrope diamon Found 1 substitutions in total -If you want to generate more BN crystal, you can first generate the subgroup representation and then apply the `substitute_1_2` function. +If you want to generate more BN crystals, you can first generate the subgroup representation and then apply the ``substitute_1_2`` function. .. code-block:: Python diff --git a/pyxtal/db.py b/pyxtal/db.py index 5101671f..a4bfc1bb 100644 --- a/pyxtal/db.py +++ b/pyxtal/db.py @@ -479,17 +479,28 @@ def __init__(self, db_name, ltol=0.05, stol=0.05, atol=3): def vacuum(self): self.db.vacuum() - def get_pyxtal(self, id, use_ff=True): + def get_pyxtal(self, id, use_relaxed=None): """ - Get pyxtal based on row_id, if use_ff, get pyxtal from the ff_relaxed file + Get pyxtal based on row_id, if use_relaxed, get pyxtal from the ff_relaxed file + + Args: + id (int): row id + use_relaxed (str): 'ff_relaxed', 'vasp_relaxed' """ from pyxtal import pyxtal from pyxtal.util import ase2pymatgen from pymatgen.core import Structure row = self.db.get(id) #; print(id, row.id) - if use_ff and hasattr(row, 'ff_relaxed'): - pmg = Structure.from_str(row.ff_relaxed, fmt='cif') + if use_relaxed is not None: + + if hasattr(row, use_relaxed): + xtal_str = getattr(row, use_relaxed) + else: + raise ValueError('No ff or vasp relaxed attributes for structure', id) + + pmg = Structure.from_str(xtal_str, fmt='cif') + else: atom = self.db.get_atoms(id=id) pmg = ase2pymatgen(atom) @@ -808,7 +819,7 @@ def clean_structures_pmg(self, ids=(None, None), min_id=None, dtol=5e-2, criteri self.db.delete(to_delete) - def select_xtals(self, ids, overwrite=False, attribute=None, use_ff=False): + def select_xtals(self, ids, overwrite=False, attribute=None, use_relaxed=None): """ Extract xtals based on attribute name. Mostly called by update_row_ff_energy or update_row_dftb_energy. @@ -816,6 +827,8 @@ def select_xtals(self, ids, overwrite=False, attribute=None, use_ff=False): Args: ids: overwrite: + atttribute: + use_relaxed (str): 'ff_relaxed' or 'vasp_relaxed' """ (min_id, max_id) = ids if min_id is None: min_id = 1 @@ -825,7 +838,7 @@ def select_xtals(self, ids, overwrite=False, attribute=None, use_ff=False): for row in self.db.select(): if overwrite or attribute is None or not hasattr(row, attribute): if min_id <= row.id <= max_id: - xtal = self.get_pyxtal(row.id, use_ff) + xtal = self.get_pyxtal(row.id, use_relaxed) ids.append(row.id) xtals.append(xtal) if len(xtals) % 100 == 0: @@ -898,6 +911,7 @@ def update_row_dftb_energy(self, skf_dir, steps=500, skf_dir (str): GULP force field library (e.g., 'reaxff', 'tersoff') steps (int): relaxation steps ids (tuple): row ids e.g., (0, 100) + use_ff (bool): use the prerelaxed ff structure or not ncpu (int): number of parallel processes calc_folder (str): temporary folder for GULP calculations symmetrize (bool): impose symmetry in optimization @@ -905,7 +919,12 @@ def update_row_dftb_energy(self, skf_dir, steps=500, """ os.makedirs(calc_folder, exist_ok=True) - ids, xtals = self.select_xtals(ids, overwrite, 'dftb_energy', use_ff) + if use_ff: + use_relaxed = 'ff_relaxed' + else: + use_relaxed = None + + ids, xtals = self.select_xtals(ids, overwrite, 'dftb_energy', use_relaxed) dftb_results = [] os.chdir(calc_folder) @@ -949,13 +968,14 @@ def update_row_dftb_energy(self, skf_dir, steps=500, dftb_relaxed=xtal.to_file()) - def update_row_topology(self, StructureType='Auto', overwrite=True): + def update_row_topology(self, StructureType='Auto', overwrite=True, prefix=None): """ Update row topology base on the CrystalNets.jl Args: StructureType (str): 'Zeolite', 'MOF' or 'Auto' overwrite (bool): remove the existing attributes + prefix (str): prefix for tmp cif file """ try: import juliacall @@ -993,13 +1013,18 @@ def parse_topology(topology_info): else: option = jl.CrystalNets.Options(structure=jl.StructureType.Auto) + if prefix is not None: + cif_file = prefix + '.cif' + else: + cif_file = 'tmp.cif' + for row in self.db.select(): if overwrite or not hasattr(row, 'topology'): atoms = self.db.get_atoms(row.id) - atoms.write('tmp.cif', format='cif') + atoms.write(cif_file, format='cif') # Call crystalnet.jl - result = jl.determine_topology('tmp.cif', option) + result = jl.determine_topology(cif_file, option) #print(result) if len(result) > 1: results = [x for x in result] @@ -1052,7 +1077,7 @@ def update_db_description(self): def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, sort_by='similarity', overwrite=True, cutoff=None, - use_ff=True): + use_relaxed=None): """ export structures from database according to the given criterion @@ -1076,12 +1101,14 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, os.makedirs(folder) keys = [ + 'id', 'pearson_symbol', 'space_group_number', 'density', 'dof', 'similarity', 'ff_energy', + 'vasp_energy', 'topology', ] properties = [] @@ -1092,34 +1119,39 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, ps = row.pearson_symbol sim = float(row.similarity) if hasattr(row, 'similarity') and row.similarity is not None else None top = row.topology if hasattr(row, 'topology') else None - eng = float(row.ff_energy) if hasattr(row, 'ff_energy') else None - properties.append([row.id, ps, spg, den, dof, sim, eng, top, ]) + ff_eng = float(row.ff_energy) if hasattr(row, 'ff_energy') else None + vasp_eng = float(row.vasp_energy) if hasattr(row, 'vasp_energy') else None + properties.append([row.id, ps, spg, den, dof, sim, ff_eng, vasp_eng, top]) + + dicts = {} + for i, key in enumerate(keys): + if properties[0][i] is not None: + dicts[key] = [prop[i] for prop in properties] if sort_by in keys: - col = keys.index(sort_by) + 1 + col = keys.index(sort_by) #+ 1 else: print("supported attributes", keys) raise ValueError("Cannot sort by", sort_by) print("====Exporting {:} structures".format(len(properties))) - #properties = np.array(properties) - #mids = np.argsort(properties[:, col])[:cutoff] - #sorted_properties = [] properties = [prop for prop in properties if prop[col] is not None] sorted_properties = sorted(properties, key=lambda x: x[col]) - #for mid in mids: for entry in sorted_properties[:cutoff]: - [id, ps, spg, den, dof, sim, eng, top] = entry + [id, ps, spg, den, dof, sim, ff_eng, vasp_eng, top] = entry id = int(id) spg = int(spg) sim = float(sim) den = float(den) dof = int(dof) - if eng is not None: eng = float(eng) - try: + if vasp_eng is not None: + eng = float(vasp_eng) + elif ff_eng is not None: + eng = float(ff_eng) #if True: - xtal = self.get_pyxtal(id, use_ff) + try: + xtal = self.get_pyxtal(id, use_relaxed) number, symbol = xtal.group.number, xtal.group.symbol.replace('/','') # convert to the desired subgroup representation if needed if number != spg: @@ -1127,7 +1159,7 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, xtal = xtal.to_subgroup(paths) number, symbol = xtal.group.number, xtal.group.symbol.replace('/','') - label = os.path.join(folder, '{:s}-{:d}-{:d}-{:s}'.format(xtal.get_Pearson_Symbol(), id, number, symbol)) + label = os.path.join(folder, '{:d}-{:s}-{:d}-{:s}'.format(id, xtal.get_Pearson_Symbol(), number, symbol)) if criteria is not None: status = xtal.check_validity(criteria, True) @@ -1137,7 +1169,6 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, status = False if status: - #if top is not None: print(top) try: #if True: xtal.set_site_coordination() @@ -1152,7 +1183,7 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, if den is not None: label += '-D{:.2f}'.format(abs(den)) if eng is not None: label += '-E{:.3f}'.format(abs(eng)) if top is not None: label += '-T{:s}'.format(top) - if sim is not None: label += '-S{:.2f}'.format(sim) + #if sim is not None: label += '-S{:.2f}'.format(sim) print("====Exporting:", label) if fmt == 'vasp': @@ -1162,6 +1193,8 @@ def export_structures(self, fmt='vasp', folder='mof_out', criteria=None, else: print("====Skippng: ", label) + return dicts + def get_label(self, i): if i < 10: folder = f"cpu00{i}" diff --git a/pyxtal/interface/vasp.py b/pyxtal/interface/vasp.py index 0754d8fe..5d41ec0d 100644 --- a/pyxtal/interface/vasp.py +++ b/pyxtal/interface/vasp.py @@ -29,7 +29,7 @@ def __init__(self, struc, path='tmp', cmd='mpirun -np 16 vasp_std'): raise NotImplementedError("only support ASE atoms object") self.structure = struc - self.folder = path + self.folder = path if not os.path.exists(self.folder): os.makedirs(self.folder) self.pstress = 0.0 @@ -41,9 +41,9 @@ def __init__(self, struc, path='tmp', cmd='mpirun -np 16 vasp_std'): self.cputime = 0 self.error = True self.cmd = cmd - + def set_vasp(self, level=0, pstress=0.0000, setup=None): - self.pstress = pstress + self.pstress = pstress default0 = {'xc': 'pbe', 'npar': 8, 'kgamma': True, @@ -95,7 +95,7 @@ def set_vasp(self, level=0, pstress=0.0000, setup=None): 'ediff': 1e-4, 'nsw': 0, } - + dict_vasp = dict(default0, **default1) return Vasp(**dict_vasp) @@ -122,7 +122,7 @@ def read_OSZICAR(self, path='OSZICAR'): def read_bandgap(self, path='vasprun.xml'): from pyxtal.interface.vasprun import vasprun myrun = vasprun(path) - self.gap = myrun.values['gap'] + self.gap = myrun.values['gap'] def run(self, setup=None, pstress=0, level=0, clean=True, read_gap=False, walltime=None): if walltime is not None: @@ -182,11 +182,11 @@ def single_optimize(struc, level, pstress, setup, path, clean, """ single optmization - Args: + Args: struc: pyxtal structure level: vasp calc level pstress: external pressure - setup: vasp setup + setup: vasp setup path: calculation directory Returns: @@ -209,11 +209,11 @@ def single_point(struc, setup=None, path=None, clean=True): """ single optmization - Args: + Args: struc: pyxtal structure level: vasp calc level pstress: external pressure - setup: vasp setup + setup: vasp setup path: calculation directory Returns: @@ -223,9 +223,9 @@ def single_point(struc, setup=None, path=None, clean=True): calc.run(setup, level=4, clean=clean) return calc.energy, calc.forces, calc.error -def optimize(struc, path, levels=[0,2,3], pstress=0, setup=None, +def optimize(struc, path, levels=[0,2,3], pstress=0, setup=None, clean=True, cmd='mpirun -np 16 vasp_std', walltime="30m"): - + """ multi optimization @@ -234,7 +234,7 @@ def optimize(struc, path, levels=[0,2,3], pstress=0, setup=None, path: calculation directory levels: list of vasp calc levels pstress: external pressure - setup: vasp setup + setup: vasp setup Returns: list of structures, energies and time costs @@ -242,7 +242,7 @@ def optimize(struc, path, levels=[0,2,3], pstress=0, setup=None, time_total = 0 for i, level in enumerate(levels): - struc, eng, time, error = single_optimize(struc, level, pstress, setup, path, + struc, eng, time, error = single_optimize(struc, level, pstress, setup, path, clean, cmd, walltime) time_total += time @@ -263,7 +263,6 @@ def optimize(struc, path, levels=[0,2,3], pstress=0, setup=None, os.system("source /share/intel/mkl/bin/mklvars.sh intel64") cmd='mpirun -n 4 /share/apps/bin/vasp544-2019u2/vasp_std' - calc = VASP(struc, path='tmp', cmd=cmd) calc.run() print("Energy:", calc.energy)