add db option and update documentation on substitute_1_2

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.

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 <pyxtal.html#pyxtal.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 <pyxtal.html#pyxtal.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
 

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,14 +819,16 @@ 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.
 
  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,14 +911,20 @@ 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
  overwrite (bool): remove the existing attributes
  """
 
  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,42 +1119,47 @@ 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:
  paths = xtal.group.path_to_subgroup(spg)
  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}"