From 3d70a8ff83e573eac5e37f56c4f5792e538be5a4 Mon Sep 17 00:00:00 2001 From: qzhu2017 Date: Sat, 17 Aug 2024 23:01:54 -0400 Subject: [PATCH] migrate mof_builder to pyxtal --- doc/Usage.rst | 4 + doc/conf.py | 4 +- doc/index.rst | 2 +- doc/pyxtal.interface.rst | 1 - doc/pyxtal.lego.SO3.rst | 7 + ...sprun.rst => pyxtal.lego.basinhopping.rst} | 4 +- doc/pyxtal.lego.builder.rst | 7 + doc/pyxtal.lego.rst | 17 + doc/pyxtal.optimize.basin_hopping.rst | 7 - doc/pyxtal.optimize.rst | 1 - doc/pyxtal.rst | 1 + pyxtal/db.py | 6 +- pyxtal/interface/LJ.py | 71 +- pyxtal/interface/vasp.py | 2 +- pyxtal/interface/vasprun.py | 899 ------------ pyxtal/lattice.py | 8 +- pyxtal/lego/SO3.py | 648 +++++++++ pyxtal/lego/__init__.py | 0 pyxtal/lego/basinhopping.py | 786 +++++++++++ pyxtal/lego/builder.py | 1212 +++++++++++++++++ requirements.txt | 6 + setup.py | 5 +- 22 files changed, 2741 insertions(+), 957 deletions(-) create mode 100644 doc/pyxtal.lego.SO3.rst rename doc/{pyxtal.interface.vasprun.rst => pyxtal.lego.basinhopping.rst} (54%) create mode 100644 doc/pyxtal.lego.builder.rst create mode 100644 doc/pyxtal.lego.rst delete mode 100644 doc/pyxtal.optimize.basin_hopping.rst delete mode 100644 pyxtal/interface/vasprun.py create mode 100644 pyxtal/lego/SO3.py create mode 100644 pyxtal/lego/__init__.py create mode 100644 pyxtal/lego/basinhopping.py create mode 100644 pyxtal/lego/builder.py diff --git a/doc/Usage.rst b/doc/Usage.rst index d11151ae..5d66fda4 100644 --- a/doc/Usage.rst +++ b/doc/Usage.rst @@ -630,6 +630,7 @@ Executing this above scripts will lead to the following output: This way, you can easily find derivative crystals in the suboptimal representations. Conversely, it is also possible to identify the likely supergroup xtal. The following snippet codes can be used to design illustrate pyxtal functionalities. .. code-block:: Python + from pyxtal import pyxtal # load a graphite crystal and make the subgroup representation @@ -957,6 +958,7 @@ between Pyxtal and its 1D representation. With this module, one can represent th c1.from_seed('pyxtal/database/cifs/aspirin.cif', ['CC(=O)OC1=CC=CC=C1C(=O)O.smi']) rep = c1.get_1D_representation() print(rep.to_string()) + :: 81 11.23 6.54 11.23 95.9 1 0.23 0.59 0.03 44.1 -25.2 32.5 82.9 2.8 -178.3 1 @@ -1003,6 +1005,7 @@ To create a new database file (e.g., `test.db`), print("Initial list of codes", db.codes) db.add_from_code('NAPHTA') print("Updated list of codes", db.codes) + :: 0 ACSALA @@ -1017,6 +1020,7 @@ To view the database file, .. code-block:: Python $ ase db test.db + :: csd_code|space_group|mol_smi diff --git a/doc/conf.py b/doc/conf.py index 256c105e..abdc9c39 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -23,9 +23,9 @@ author = "Qiang Zhu, Scott Fredericks, Kevin Parrish" # The short X.Y version -version = "1.0.1" +version = "1.0.2" # The full version, including alpha/beta/rc tags -release = "1.0.1" +release = "1.0.2" # -- General configuration --------------------------------------------------- diff --git a/doc/index.rst b/doc/index.rst index f200e187..5a4f6adb 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -24,7 +24,7 @@ ab-initio generation of random crystal structures. It has the following features - Structural manipulation via symmetry relation (both subgroup and supergroup) - Geometry optimization from built-in and external optimization methods -The current version is ``1.0.1`` at `GitHub `_. +The current version is ``1.0.2`` at `GitHub `_. It is available for use under the MIT license. Expect updates upon request by `Qiang Zhu\'s group `_ at the University of North Carolina at Charlotte. diff --git a/doc/pyxtal.interface.rst b/doc/pyxtal.interface.rst index c2c58377..f8ac534f 100644 --- a/doc/pyxtal.interface.rst +++ b/doc/pyxtal.interface.rst @@ -19,4 +19,3 @@ Submodules pyxtal.interface.gulp pyxtal.interface.lammpslib pyxtal.interface.vasp - pyxtal.interface.vasprun diff --git a/doc/pyxtal.lego.SO3.rst b/doc/pyxtal.lego.SO3.rst new file mode 100644 index 00000000..c93622e2 --- /dev/null +++ b/doc/pyxtal.lego.SO3.rst @@ -0,0 +1,7 @@ +pyxtal.lego.SO3 module +========================== + +.. automodule:: pyxtal.lego.SO3 + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/pyxtal.interface.vasprun.rst b/doc/pyxtal.lego.basinhopping.rst similarity index 54% rename from doc/pyxtal.interface.vasprun.rst rename to doc/pyxtal.lego.basinhopping.rst index 50f7eff7..2f83b898 100644 --- a/doc/pyxtal.interface.vasprun.rst +++ b/doc/pyxtal.lego.basinhopping.rst @@ -1,7 +1,7 @@ -pyxtal.interface.vasprun module +pyxtal.lego.basinhopping module =============================== -.. automodule:: pyxtal.interface.vasprun +.. automodule:: pyxtal.lego.basinhopping :members: :undoc-members: :show-inheritance: diff --git a/doc/pyxtal.lego.builder.rst b/doc/pyxtal.lego.builder.rst new file mode 100644 index 00000000..2072f161 --- /dev/null +++ b/doc/pyxtal.lego.builder.rst @@ -0,0 +1,7 @@ +pyxtal.lego.builder module +========================== + +.. automodule:: pyxtal.lego.builder + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/pyxtal.lego.rst b/doc/pyxtal.lego.rst new file mode 100644 index 00000000..71b22289 --- /dev/null +++ b/doc/pyxtal.lego.rst @@ -0,0 +1,17 @@ +pyxtal.lego package +=================== + +.. automodule:: pyxtal.lego + :members: + :undoc-members: + :show-inheritance: + +Submodules +---------- + +.. toctree:: + :maxdepth: 4 + + pyxtal.lego.builder + pyxtal.lego.SO3 + pyxtal.lego.basinhopping diff --git a/doc/pyxtal.optimize.basin_hopping.rst b/doc/pyxtal.optimize.basin_hopping.rst deleted file mode 100644 index 1210a596..00000000 --- a/doc/pyxtal.optimize.basin_hopping.rst +++ /dev/null @@ -1,7 +0,0 @@ -pyxtal.optimize.basin\_hopping module -===================================== - -.. automodule:: pyxtal.optimize.basin_hopping - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pyxtal.optimize.rst b/doc/pyxtal.optimize.rst index 01c91acd..50f87a55 100644 --- a/doc/pyxtal.optimize.rst +++ b/doc/pyxtal.optimize.rst @@ -16,6 +16,5 @@ Submodules pyxtal.optimize.QRS pyxtal.optimize.WFS pyxtal.optimize.base - pyxtal.optimize.basin_hopping pyxtal.optimize.benchmark pyxtal.optimize.common diff --git a/doc/pyxtal.rst b/doc/pyxtal.rst index bd37a3b3..a1681a30 100644 --- a/doc/pyxtal.rst +++ b/doc/pyxtal.rst @@ -14,6 +14,7 @@ Subpackages pyxtal.database pyxtal.interface + pyxtal.lego pyxtal.optimize pyxtal.potentials diff --git a/pyxtal/db.py b/pyxtal/db.py index 0454dfd9..e72f0c5b 100644 --- a/pyxtal/db.py +++ b/pyxtal/db.py @@ -700,9 +700,9 @@ def clean_structures(self, ids=(None, None), dtol=2e-3, etol=1e-3, criteria=None """ Clean up the db by removing the duplicate structures Here we check the follow criteria - - same number of atoms - - same density - - same energy + - same number of atoms + - same density + - same energy Args: dtol (float): tolerance of density diff --git a/pyxtal/interface/LJ.py b/pyxtal/interface/LJ.py index c06c659b..784c912d 100644 --- a/pyxtal/interface/LJ.py +++ b/pyxtal/interface/LJ.py @@ -311,38 +311,39 @@ def symmetrized_stress(self, stress): return m2 -from spglib import get_symmetry_dataset - -from pyxtal.crystal import random_crystal - -for _i in range(10): - crystal = random_crystal(11, ["C"], [4], 1.0) - if crystal.valid: - crystal1 = deepcopy(crystal) - test = LJ(epsilon=0.01, sigma=3.40, rcut=8.0) - struc = (crystal1.lattice_matrix, crystal1.frac_coords, [6] * 4) - eng, enth, force, stress = test.calc(crystal1) - sg = get_symmetry_dataset(struc)["number"] - print(f"\nBefore relaxation Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") - - dyn1 = FIRE(crystal1, test, f_tol=1e-5, dt=0.2, maxmove=0.2) # , symmetrize=True) - dyn1.run(500) - eng, enth, force, stress = test.calc(crystal1) - struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6] * 4) - sg = get_symmetry_dataset(struc, symprec=0.1)["number"] - print(f"After relaxation without symm Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") - - dyn1 = FIRE(crystal, test, f_tol=1e-5, dt=0.2, maxmove=0.2, symmetrize=True) - dyn1.run(500) - eng, enth, force, stress = test.calc(crystal) - struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6] * 4) - sg = get_symmetry_dataset(struc, symprec=0.02)["number"] - if sg is None: - sg = 0 - print(f"After relaxation with symm Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") - # dyn1 = FIRE(crystal, test, f_tol=1e-5, dt=0.2, maxmove=0.2) - # struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6]*4) - # sg = get_symmetry_dataset(struc, symprec=0.02)['number'] - # print('After relaxation without symm Space group: {:4d} Energy: {:12.4} Enthalpy: {:12.4}'.format(sg, eng, enth)) - # right now, it seems structures goes to either HCP of FCC after relaxation, which is expected for 3D LJ system - # need to compare with other code to see if the energy is correct +if __name__ == "__main__": + from spglib import get_symmetry_dataset + from pyxtal import pyxtal + + for _i in range(10): + crystal = pyxtal() + crystal.from_random(3, 11, ["C"], [4]) + if crystal.valid: + crystal1 = deepcopy(crystal) + test = LJ(epsilon=0.01, sigma=3.40, rcut=8.0) + struc = (crystal1.lattice_matrix, crystal1.frac_coords, [6] * 4) + eng, enth, force, stress = test.calc(crystal1) + sg = get_symmetry_dataset(struc)["number"] + print(f"\nBefore relaxation Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") + + dyn1 = FIRE(crystal1, test, f_tol=1e-5, dt=0.2, maxmove=0.2) # , symmetrize=True) + dyn1.run(500) + eng, enth, force, stress = test.calc(crystal1) + struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6] * 4) + sg = get_symmetry_dataset(struc, symprec=0.1)["number"] + print(f"After relaxation without symm Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") + + dyn1 = FIRE(crystal, test, f_tol=1e-5, dt=0.2, maxmove=0.2, symmetrize=True) + dyn1.run(500) + eng, enth, force, stress = test.calc(crystal) + struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6] * 4) + sg = get_symmetry_dataset(struc, symprec=0.02)["number"] + if sg is None: + sg = 0 + print(f"After relaxation with symm Space group: {sg:4d} Energy: {eng:12.4} Enthalpy: {enth:12.4}") + # dyn1 = FIRE(crystal, test, f_tol=1e-5, dt=0.2, maxmove=0.2) + # struc = (dyn1.struc.lattice_matrix, dyn1.struc.frac_coords, [6]*4) + # sg = get_symmetry_dataset(struc, symprec=0.02)['number'] + # print('After relaxation without symm Space group: {:4d} Energy: {:12.4} Enthalpy: {:12.4}'.format(sg, eng, enth)) + # right now, it seems structures goes to either HCP of FCC after relaxation, which is expected for 3D LJ system + # need to compare with other code to see if the energy is correct diff --git a/pyxtal/interface/vasp.py b/pyxtal/interface/vasp.py index 5e9f112e..27419bdd 100644 --- a/pyxtal/interface/vasp.py +++ b/pyxtal/interface/vasp.py @@ -131,7 +131,7 @@ def read_OSZICAR(self, path="OSZICAR"): self.energy = energy # this is actually enthalpy def read_bandgap(self, path="vasprun.xml"): - from pyxtal.interface.vasprun import vasprun + from vasprun import vasprun myrun = vasprun(path) self.gap = myrun.values["gap"] diff --git a/pyxtal/interface/vasprun.py b/pyxtal/interface/vasprun.py deleted file mode 100644 index 05b21a0d..00000000 --- a/pyxtal/interface/vasprun.py +++ /dev/null @@ -1,899 +0,0 @@ -import matplotlib as mpl -import numpy as np -import pandas as pd -from lxml import etree - -mpl.use("Agg") -import contextlib - -import matplotlib.pyplot as plt -from matplotlib import rcParams - -rcParams.update({"figure.autolayout": True}) -plt.style.use("bmh") - - -def smear_data(data, sigma): - """ - Apply Gaussian smearing to spectrum y value. - - Args: - sigma: Std dev for Gaussian smear function - """ - from scipy.ndimage.filters import gaussian_filter1d - - diff = [data[i + 1, 0] - data[i, 0] for i in range(len(data) - 1)] - avg_x_per_step = np.sum(diff) / len(diff) - data[:, 1] = gaussian_filter1d(data[:, 1], sigma / avg_x_per_step) - return data - - -class units: - am2kg = 1.6605402e-27 - ev2j = 1.60217733e-19 - plank = 6.626075e-34 - c = 2.99792458e10 - pi = 3.1415926 - proton_mass = 1836 - ev2hartree = 27.211386245988 - a2bohr = 0.529 - ev2cm = 8065.6 - - -class vasprun: - """ - parse vasprun.xml and return all useful info to self.values - - Args: - vasp_file: the path of vasprun.xml - verbosity: output error msgs or not - """ - - def __init__(self, vasp_file="vasprun.xml", verbosity=0): - self.error = False - self.errormsg = "" - self.values = {} - try: - doc = etree.parse(vasp_file) - doc = doc.getroot() - self.parse_vaspxml(doc) - self.get_band_gap() - N_atom = len(self.values["finalpos"]["positions"]) - self.values["calculation"]["energy_per_atom"] = self.values["calculation"]["energy"] / N_atom - except etree.XMLSyntaxError: - self.error = True - self.errormsg = "corrupted file found" - - if verbosity > 0 and self.error is True: - print("----------Warning---------------") - print(self.errormsg) - print("--------------------------------") - - def parse_vaspxml(self, xml_object): - """ - parse the following tags - - - incar - - kpoints - - atominfo - composition, elements, formula - - calculation - eigenvalues, energy, dos, fermi energy, stress, force - - finalpos - lattice, rec_lat, positions - """ - - for child in xml_object.iterchildren(): - self.values.setdefault(child.tag, {}) - if child.tag == "incar": - self.values[child.tag] = self.parse_i_tag_collection(child) - elif child.tag == "kpoints": - self.values[child.tag] = self.parse_kpoints(child) - elif child.tag == "parameters": - self.values[child.tag] = self.parse_parameters(child) - elif child.tag == "atominfo": - self.values["name_array"] = self.parse_name_array(child) - self.values["composition"] = self.parse_composition(child) - self.values["elements"] = self.get_element(self.values["composition"]) - self.values["formula"] = self.get_formula(self.values["composition"]) - ( - self.values["pseudo_potential"], - self.values["potcar_symbols"], - self.values["valence"], - self.values["mass"], - ) = self.get_potcar(child) - elif child.tag == "calculation": - self.values["calculation"], scf_count = self.parse_calculation(child) - if self.values["parameters"]["electronic"]["electronic convergence"]["NELM"] == scf_count: - self.error = True - self.errormsg = "SCF is not converged" - - if ( - self.values["parameters"]["electronic"]["electronic spin"]["LSORBIT"] - or self.values["parameters"]["electronic"]["electronic spin"]["ISPIN"] == 2 - ): - self.spin = True - else: - self.spin = False - elif child.tag == "structure" and child.attrib.get("name") == "finalpos": - self.values["finalpos"] = self.parse_finalpos(child) - elif child.tag not in ( - "i", - "r", - "v", - "incar", - "kpoints", - "atominfo", - "calculation", - ): - self.values[child.tag] = self.parse_vaspxml(child) - # else: - # return 1 - self.dict_clean(self.values) - - @staticmethod - def dict_clean(dict_del): - """ - Delete the keys that is {} and None - Loops recursively over nested dictionaries. - """ - dict_foo = dict_del.copy() # Used as iterator to avoid the 'DictionaryHasChanged' error - for key in dict_foo: - if isinstance(dict_foo[key], dict): - vasprun.dict_clean(dict_del[key]) - - # if len(dict_foo[key])==0: - if dict_foo[key] == {} or dict_foo[key] is None: - with contextlib.suppress(KeyError): - del dict_del[key] - - return dict_del - - def parse_finalpos(self, finalpos): - """ - obtain final configuration - """ - d = {} - for i in finalpos.iter("varray"): - name = i.attrib.get("name") - d[name] = self.parse_varray(i) - return d - - def parse_dynmat(self, dynmat): - hessian, eigenvalues, eigenvectors = [], [], [] - for va in dynmat.findall("varray"): - if va.attrib.get("name") == "hessian": - hessian = self.parse_varray(va) - elif va.attrib.get("name") == "eigenvectors": - eigenvectors = self.parse_varray(va) - factor = np.sqrt(units.ev2j / 1e-20 / units.am2kg) - for v in dynmat.findall("v"): - for i in v.text.split(): - eigenvalues.append(np.sqrt(abs(float(i))) * factor * units.plank / units.ev2j / 2 / units.pi) - # eigenvalues.append(np.sqrt(abs(float(i)))) - return hessian, eigenvalues, eigenvectors - - def parse_born_chg(self, charge): - """ - obtain the born charge - """ - chg = [] - for info in charge.findall("set"): - chg.append(self.parse_varray(info)) - return chg - - def parse_i_tag_collection(self, itags_collection): - d = {} - for info in itags_collection.findall("i"): - name = info.attrib.get("name") - type = info.attrib.get("type") - content = info.text - d[name] = self.assign_type(type, content) - return d - - @staticmethod - def parse_varray_pymatgen(elem): - def _vasprun_float(f): - """ - Large numbers are often represented as ********* in the vasprun. - This function parses these values as np.nan - """ - try: - return float(f) - except ValueError as e: - f = f.strip() - if f == "*" * len(f): - warnings.warn("Float overflow (*******) encountered in vasprun") - return np.nan - raise e - - if elem.get("type", None) == "logical": - m = [[i == "T" for i in v.text.split()] for v in elem] - else: - m = [[_vasprun_float(i) for i in v.text.split()] for v in elem] - - return m - - @staticmethod - def parse_varray(varray): - if varray.get("type") == "int": - m = [[int(number) for number in v.text.split()] for v in varray.findall("v")] - else: - try: - m = [[float(number) for number in v.text.split()] for v in varray.findall("v")] - except ValueError: - m = [[0 for number in v.text.split()] for v in varray.findall("v")] - return m - - @staticmethod - def parse_array(array): - array_dictionary = {} - values = [] - dimension_list = {} - field_list = [] - - for dimension in array.findall("dimension"): - dimension_list[dimension.attrib.get("dim")] = dimension.text - - for field in array.findall("field"): - field_list.append(field.text) - - for r in array.iter("r"): - values.append([float(number) for number in r.text.split()]) - - array_dictionary["value"] = values - array_dictionary["dimensions"] = dimension_list - array_dictionary["fields"] = field_list - - return array_dictionary - - @staticmethod - def assign_type(type, content): - if type == "logical": - content = content.replace(" ", "") - if content in ("T", "True", "true"): - return True - elif content in ("F", "False", "false"): - return False - else: - Warning("logical text " + content + " not T, True, true, F, False, false, set to False") - return False - elif type == "int": - return int(content) if len(content.split()) == 1 else [int(number) for number in content.split()] - elif type == "string": - return content - elif type is None: - return float(content) if len(content.split()) == 1 else [float(number) for number in content.split()] - else: - Warning("New type: " + type + ", set to string") - return content - - @staticmethod - def parse_composition(atom_info): - atom_names = {} - for set in atom_info.findall("array"): - if set.attrib.get("name") == "atoms": - for rc in set.iter("rc"): - atom_name = rc.find("c").text.replace(" ", "") - if atom_name in atom_names: - atom_names[atom_name] += 1 - else: - atom_names[atom_name] = 1 - - break - return atom_names - - @staticmethod - def get_element(atom_names_dictionary): - elements = [] - for atom_name in atom_names_dictionary: - elements.append(atom_name.replace(" ", "")) - return elements - - @staticmethod - def get_formula(atom_names_dictionary): - formula = "" - for atom_name in atom_names_dictionary: - formula += atom_name.replace(" ", "") + str(atom_names_dictionary[atom_name]) - return formula - - def get_potcar(self, child): - """ - parse the potcar information - - - {'labels': ['O', 'Sr_sv'], 'pot_type': 'paw', 'functional': 'pbe'} - - ['PAW_PBE', 'N', '08Apr2002'] - """ - pseudo = {"labels": [], "pot_type": [], "functional": []} - potcar_symbol = [] - valence = [] - mass = [] - for i in child.iterchildren(): - if i.tag == "array" and i.attrib.get("name") == "atomtypes": - ll = list(i.iter("c")) - for i in range(3, len(ll), 5): - valence.append(float(ll[i].text)) - - for i in range(2, len(ll), 5): - mass.append(float(ll[i].text)) - - for i in range(4, len(ll), 5): - text = ll[i].text.split() - label = text[1] - pot = text[0].split("_")[0] - try: - xc = text[0].split("_")[1] - except: - xc = "unknown" - pseudo["labels"].append(label) - pseudo["pot_type"].append(pot) - pseudo["functional"].append(xc) - potcar_symbol.append(xc + " " + label) - return pseudo, potcar_symbol, valence, mass - - @staticmethod - def parse_name_array(atominfo): - atom_names = [] - for array in atominfo.findall("array"): - if array.attrib["name"] == "atoms": - atom_names = [rc.find("c").text.strip() for rc in array.find("set")] - - if atom_names == []: - ValueError("No atomname found in file") - - return atom_names - - # def parse_eigenvalue(self, eigenvalue): - # eigenvalue = eigenvalue.find("array") - # eigenvalues = self.parse_array(eigenvalue) - # return eigenvalues - def parse_parameters(self, child): - parameters = {} - for i in child: - if i.tag == "separator": - name = i.attrib.get("name") - d = self.parse_i_tag_collection(i) - parameters[name] = d - for ii in i: - if ii.tag == "separator": - name2 = ii.attrib.get("name") - d2 = self.parse_i_tag_collection(ii) - parameters[name][name2] = d2 - return parameters - - def parse_eigenvalue(self, eigenvalue): - eigenvalues = [] - for s in eigenvalue.find("array").find("set").findall("set"): - for ss in s.findall("set"): - eigenvalues.append(self.parse_varray_pymatgen(ss)) - return eigenvalues - - def parse_dos(self, dos): - t_dos = [] - p_dos = [] - for s in dos.find("total").find("array").findall("set"): - for ss in s.findall("set"): - t_dos.append(self.parse_varray_pymatgen(ss)) - if dos.find("partial") is not None and len(dos.find("partial")) > 0: - for s in dos.find("partial").find("array").findall("set"): - for _i, ss in enumerate(s.findall("set")): - p = [] - for sss in ss.findall("set"): - p.append(self.parse_varray_pymatgen(sss)) - p_dos.append(p) - - return t_dos, p_dos - - def parse_projected(self, proj): - projected = [] - for s in proj.find("array").find("set").findall("set"): - for ss in s.findall("set"): - p = [] - for sss in ss.findall("set"): - p.append(self.parse_varray_pymatgen(sss)) - projected.append(p) - - return projected # [N_kpts, N_bands] - - def parse_calculation(self, calculation): - stress = [] - force = [] - efermi = 0.0 - eigenvalues = [] - energy = 0.0 - scf_count = 0 - tdos = [] - pdos = [] - born_chgs = [] - hessian = [] - dyn_eigenvalues = [] - dyn_eigenvectors = [] - epsilon_ion = [] - proj = [] - for i in calculation.iterchildren(): - if i.attrib.get("name") == "stress": - stress = self.parse_varray(i) - elif i.attrib.get("name") == "forces": - force = self.parse_varray(i) - elif i.tag == "dos": - for j in i.findall("i"): - if j.attrib.get("name") == "efermi": - efermi = float(j.text) - break - tdos, pdos = self.parse_dos(i) - elif i.tag == "projected": - proj = self.parse_projected(i) - elif i.tag == "eigenvalues": - eigenvalues = self.parse_eigenvalue(i) - elif i.tag == "scstep": - for j in i.iterchildren(): - if j.tag == "energy": - for e in j.findall("i"): - if e.attrib.get("name") == "e_fr_energy": - scf_count += 1 - - elif i.tag == "energy": - for e in i.findall("i"): - if e.attrib.get("name") == "e_fr_energy": - try: - energy = float(e.text) - except ValueError: - energy = 100000000 - else: - Warning("No e_fr_energy found in tag, energy set to 0.0") - elif i.tag == "array" and i.attrib.get("name") == "born_charges": - born_chgs = self.parse_born_chg(i) - elif i.tag == "varray" and i.attrib.get("name") == "epsilon_ion": - epsilon_ion = self.parse_varray(i) - - elif i.tag == "dynmat": - hessian, dyn_eigenvalues, dyn_eigenvectors = self.parse_dynmat(i) - - calculation = {} - calculation["stress"] = stress - calculation["efermi"] = efermi - calculation["force"] = force - calculation["eband_eigenvalues"] = eigenvalues - calculation["energy"] = energy - calculation["tdos"] = tdos - calculation["pdos"] = pdos - calculation["projected"] = proj - calculation["born_charges"] = born_chgs - calculation["hessian"] = hessian - calculation["normal_modes_eigenvalues"] = dyn_eigenvalues - calculation["normal_modes_eigenvectors"] = dyn_eigenvectors - calculation["epsilon_ion"] = epsilon_ion - return calculation, scf_count - - def parse_kpoints(self, kpoints): - kpoints_dict = {"list": [], "weights": [], "divisions": [], "mesh_scheme": ""} - - for i in kpoints.iterchildren(): - if i.tag == "generation": - kpoints_dict["mesh_scheme"] = i.attrib.get("param") - for j in i.iterchildren(): - if j.attrib.get("name") == "divisions": - kpoints_dict["divisions"] = [int(number) for number in j.text.split()] - break - - for va in kpoints.findall("varray"): - name = va.attrib["name"] - if name == "kpointlist": - kpoints_dict["list"] = self.parse_varray(va) - elif name == "weights": - kpoints_dict["weights"] = self.parse_varray(va) - - return kpoints_dict - - def get_bands(self): - """ - Function for computing the valence band index from the count of electrons - - Args: - None - - Returns: - bands: an integer number - occupy: bool number - """ - self.values["valence"] - self.values["composition"] - total = int(self.values["parameters"]["electronic"]["NELECT"]) - - # if self.spin: - # fac = 1 - # else: - # fac = 2 - fac = 2 - - if total % 2 == 0: - IBAND = int(total / fac) - occupy = True - else: - IBAND = int(total / fac) + 1 - occupy = False - - self.values["bands"] = IBAND - self.values["occupy"] = occupy - - @staticmethod - def get_cbm(kpoints, efermi, eigens, IBAND): - ind = np.argmin(eigens[:, IBAND, 0]) - pos = kpoints[ind] - value = eigens[ind, IBAND, 0] - efermi - return {"kpoint": pos, "value": value} - - @staticmethod - def get_vbm(kpoints, efermi, eigens, IBAND): - ind = np.argmax(eigens[:, IBAND, 0]) - pos = kpoints[ind] - value = eigens[ind, IBAND, 0] - efermi - return {"kpoint": pos, "value": value} - - def get_band_gap(self): - self.get_bands() - IBAND = self.values["bands"] - occupy = self.values["occupy"] - self.values["metal"] = False - self.values["gap"] = None - self.values["cbm"] = None - self.values["vbm"] = None - if occupy is True: - efermi = self.values["calculation"]["efermi"] - eigens = np.array(self.values["calculation"]["eband_eigenvalues"]) - kpoints = np.array(self.values["kpoints"]["list"]) - if np.shape(eigens)[0] > np.shape(kpoints)[0]: - kpoints = np.tile(kpoints, [2, 1]) - - cbm = self.get_cbm(kpoints, efermi, eigens, IBAND) - vbm = self.get_vbm(kpoints, efermi, eigens, IBAND - 1) - self.values["gap"] = cbm["value"] - vbm["value"] - self.values["cbm"] = cbm - self.values["vbm"] = vbm - if self.values["gap"] < 0: - self.values["metal"] = True - self.values["gap"] = 0 - else: - self.values["metal"] = True - self.values["gap"] = 0 - - def eigenvalues_by_band(self, band=0): - efermi = self.values["calculation"]["efermi"] - eigens = np.array(self.values["calculation"]["eband_eigenvalues"]) - return eigens[:, band, 0] - efermi - - def show_eigenvalues_by_band(self, bands=None): - if bands is None: - bands = [0] - kpts = self.values["kpoints"]["list"] - col_name = {"K-points": kpts} - for band in bands: - eigen = self.eigenvalues_by_band(band) - if self.spin: - eigens = np.reshape(eigen, [int(len(eigen) / 2), 2]) - name1 = "band" + str(band) + "up" - name2 = "band" + str(band) + "down" - col_name[name1] = eigens[:, 0] - col_name[name2] = eigens[:, 1] - else: - name = "band" + str(band) - col_name[name] = eigen - df = pd.DataFrame(col_name) - print(df) - - def export_incar(self, filename=None, print_incar=True): - """export incar""" - contents = [] - for key in self.values["incar"]: - content = key + " = " + str(self.values["incar"][key]) - content += "\n" - contents.append(str(content)) - if filename is not None: - with open(filename, "w") as f: - f.writelines(contents) - elif print_incar: - print(contents) - self.incar = contents - - def export_kpoints(self, filename=None): - """export kpoints""" - contents = ["KPOINTS\n"] - contents += str(len(self.values["kpoints"]["list"])) + "\n" - contents += ["Cartesian\n"] - for kpt, wt in zip(self.values["kpoints"]["list"], self.values["kpoints"]["weights"]): - content = f"{kpt[0]:10.4f} {kpt[1]:10.4f} {kpt[2]:10.4f} {wt[0]:10.4f}" - if filename is None: - print(content) - else: - content += "\n" - contents.append(str(content)) - if filename is not None: - with open(filename, "w") as f: - f.writelines(contents) - - def export_poscar(self, filename): - """ - export poscar - - Args: - filename: string - Returns: - a POSCAR file - """ - - comp = self.values["composition"] - self.values["name_array"] - latt = self.values["finalpos"]["basis"] - pos = self.values["finalpos"]["positions"] - - with open(filename, "w") as f: - string = "" - for key in comp: - string += key - string += str(comp[key]) - string += "\n" - f.write(string) - f.write("1.0\n") - f.write(f"{latt[0][0]:12.6f} {latt[0][1]:12.6f} {latt[0][2]:12.6f}\n") - f.write(f"{latt[1][0]:12.6f} {latt[1][1]:12.6f} {latt[1][2]:12.6f}\n") - f.write(f"{latt[2][0]:12.6f} {latt[2][1]:12.6f} {latt[2][2]:12.6f}\n") - for key in comp: - f.write(f"{key:4s}") - f.write("\n") - for key in comp: - f.write(f"{comp[key]:4d}") - f.write("\n") - f.write("Direct\n") - for coor in pos: - f.write(f"{coor[0]:12.6f} {coor[1]:12.6f} {coor[2]:12.6f}\n") - - def parse_bandpath(self): - kpts = self.values["kpoints"]["list"] - rec_basis = np.array(self.values["finalpos"]["rec_basis"]) - - def inline(kpt, path): - if len(path) < 2: - return True - else: - v1 = np.array(kpt) - np.array(path[-1]) - v2 = np.array(path[-1]) - np.array(path[-2]) - v1_norm = np.linalg.norm(v1) - v2_norm = np.linalg.norm(v2) - if v1_norm < 1e-3: - return False - else: - cos = np.dot(v1, v2) / v1_norm / v2_norm - if abs(cos - 1) < 1e-2: - return True - else: - # print(kpt, path[-2], path[-1], angle, np.dot(v1, v2)/v1_norm/v2_norm) - return False - - paths = [] - path = [] - for i, kpt in enumerate(kpts): - if inline(kpt, path): - path.append(kpt) - if i == len(kpts) - 1: - paths.append(path) - else: - paths.append(path) - path = [] - path.append(kpt) - - band_points = [] - pointer = 0 - for i, path in enumerate(paths): - path = np.array(path) - dist = np.linalg.norm(np.dot(path[0, :] - path[-1, :], rec_basis)) - x = np.linspace(pointer, pointer + dist, len(path[:, 0])) - band_paths = x if i == 0 else np.hstack((band_paths, x)) - pointer += dist - band_points.append(pointer) - - self.values["band_paths"] = band_paths - self.values["band_points"] = band_points - - def plot_band( - self, - filename=None, - styles="normal", - ylim=None, - plim=None, - saveBands=False, - dpi=300, - ): - """ - plot the bandstructure - - Args: - filename: string - styles: string (`normal` or `projected`) - ylim: list, the range of energy values on the y-axis, e.g. [-5, 3] - p_max: float (the ratio of color plot in the `projected` mode) - - Returns: - A figure with band structure - """ - if plim is None: - plim = [0.0, 0.5] - if ylim is None: - ylim = [-20, 3] - self.parse_bandpath() - efermi = self.values["calculation"]["efermi"] - eigens = np.array(self.values["calculation"]["eband_eigenvalues"]) - paths = self.values["band_paths"] - band_pts = self.values["band_points"] - proj = np.array(self.values["calculation"]["projected"]) # [N_kpts, N_band, Ions, 9] - cm = plt.cm.get_cmap("RdYlBu") - nkpt, nband, nocc = np.shape(eigens) - for i in range(nband): - band = eigens[:, i, 0] - efermi - if np.all(band < ylim[0]) or np.all(band > ylim[1]): - continue - p = np.empty([len(paths)]) - for kpt, _ in enumerate(paths): - p[kpt] = np.sum(proj[kpt, i, :, :]) - if len(band) / len(paths) == 2: - plt.plot(paths, band[: len(paths)], c="black", lw=1.0) - plt.plot(paths, band[len(paths) :], c="red", lw=1.0) - else: - plt.plot(paths, band, c="black", lw=1.0) - if styles == "projected": - p[p < plim[0]] = plim[0] - p[p > plim[1]] = plim[1] - plt.scatter(paths, band, c=p, vmin=plim[0], vmax=plim[1], cmap=cm, s=10) - if saveBands: - np.savetxt("band%04d.dat" % i, np.transpose([band, p])) - else: - if saveBands: - np.savetxt("band%04d.dat" % i, band) - - for pt in band_pts: - plt.axvline(x=pt, ls="-", color="k", alpha=0.5) - - if styles == "projected": - cbar = plt.colorbar( - orientation="horizontal", - ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], - ) - cbar.set_label("ratio of projected DOS") # , fontsize=12) - plt.ylabel("Energy (eV)") - plt.ylim(ylim) - plt.xlim([0, paths[-1]]) - plt.xticks([]) - if filename is None: - plt.show() - else: - plt.savefig(filename, dpi=dpi) - plt.close() - - def get_dos(self, rows, style="t"): - mydos = [] - labels = [] - a_array = self.values["name_array"] - N_atom = len(a_array) - tdos = np.array(self.values["calculation"]["tdos"]) - pdos = np.array(self.values["calculation"]["pdos"]) - a, b, c, d = np.shape(pdos) - pdos = np.reshape(pdos, [b, a, c, d]) - if style == "t": - for spin in tdos: - mydos.append(spin[rows, 1]) - labels.append("total") - - elif style in ["s", "p", "d"]: - for spin in pdos: - spd = spin[0, :, :] - for i in range(1, N_atom): - spd += spin[i, :, :] - - if style == "s": - mydos.append(spd[rows, 1]) - elif style == "p": - mydos.append(spd[rows, 2] + spd[rows, 3] + spd[rows, 4]) - else: - mydos.append(spd[rows, 5] + spd[rows, 6] + spd[rows, 7] + spd[rows, 8] + spd[rows, 9]) - labels.append(style) - - elif style[0] == "a": - if style[1].isdigit(): - ids = style[1].split("-") - start, end = int(ids[0]), int(ids[1]) - ids = range(start, end + 1) - else: - ele = style[1:] - ids = [] - for i in range(N_atom): - if ele == a_array[i]: - ids.append(i) - for spin in pdos: - spd = spin[ids[0], :, :] - for i in ids[1:]: - spd += spin[i, :, :] - mydos.append( - spd[rows, 1] - + spd[rows, 2] - + spd[rows, 3] - + spd[rows, 4] - + spd[rows, 5] - + spd[rows, 6] - + spd[rows, 7] - + spd[rows, 8] - + spd[rows, 9] - ) - labels.append(style[1:]) - - if len(labels) == 2: - labels[0] += "-up" - labels[1] += "-down" - mydos[1] *= -1 - return mydos, labels - - def plot_dos(self, filename=None, smear=None, styles="t", xlim=None, dpi=300): - """ - plot the DOS - - Args: - filename: string - styles: string (`t` or `s` or `t+spd`) - xlim: list, the range of energy values on the x-axis, e.g. [-5, 3] - smear: float (the width of smearing, defult: None) - - Returns: - A figure with band structure - """ - if xlim is None: - xlim = [-3, 3] - efermi = self.values["calculation"]["efermi"] - tdos = np.array(self.values["calculation"]["tdos"][0]) - tdos[:, 0] -= efermi - e = tdos[:, 0] - rows = (e > xlim[0]) & (e < xlim[1]) - e = e[rows] - plt_obj = {} - for option in styles.split("+"): - option = ["s", "p", "d"] if option == "spd" else [option] - for style in option: - mydos, labels = self.get_dos(rows, style) - for data, label in zip(mydos, labels): - plt_obj[label] = data - - fig, ax = plt.subplots() - lines1 = [] - lines2 = [] - labels1 = [] - labels2 = [] - for label in plt_obj: - e = np.reshape(e, [len(e), 1]) - data = np.reshape(plt_obj[label], [len(e), 1]) - if smear is not None: - data = np.hstack((e, data)) - data = smear_data(data, smear) - data = data[:, 1] - if label.find("down") > 0: - lines2 += ax.plot(e, data) - labels2.append(label) - else: - lines1 += ax.plot(e, data) - labels1.append(label) - leg1 = ax.legend(lines1, list(labels1), fancybox=True, loc="upper right") - leg1.get_frame().set_alpha(0.5) - if len(lines2) > 0: - from matplotlib.legend import Legend - - leg2 = Legend( - ax, - lines2, - list(labels2), - fancybox=True, - loc="lower right", - ) - ax.add_artist(leg2) - leg2.get_frame().set_alpha(0.5) - - plt.xlabel("Energy (eV)") - plt.ylabel("DOS") - plt.xlim(xlim) - if filename is None: - plt.show() - else: - plt.savefig(filename, dpi=dpi) - plt.close() diff --git a/pyxtal/lattice.py b/pyxtal/lattice.py index 7e42f031..90ad5fb1 100644 --- a/pyxtal/lattice.py +++ b/pyxtal/lattice.py @@ -28,13 +28,13 @@ class Lattice: kwargs: various values which may be defined. If none are defined, random ones will be generated. Values will be passed to generate_lattice. Options include: - area: The cross-sectional area (in Ang^2). Only for 1D crystals - thickness: The cell's thickness (in Angstroms) for 2D crystals - unique_axis: The unique axis for certain symmetry (and especially + 'area': The cross-sectional area (in Ang^2). Only for 1D crystals + 'thickness': The cell's thickness (in Angstroms) for 2D crystals + 'unique_axis': The unique axis for certain symmetry (and especially layer) groups. Because the symmetry operations are not also transformed, you should use the default values for random crystal generation - random: If False, keeps the stored values for the lattice geometry + 'random': If False, keeps the stored values for the lattice geometry even upon applying reset_matrix. To alter the matrix, use set_matrix() or set_para 'unique_axis': the axis ('a', 'b', or 'c') which is not symmetrically diff --git a/pyxtal/lego/SO3.py b/pyxtal/lego/SO3.py new file mode 100644 index 00000000..1cc84845 --- /dev/null +++ b/pyxtal/lego/SO3.py @@ -0,0 +1,648 @@ +from __future__ import division +import numpy as np +from scipy.special import sph_harm, spherical_in +from ase import Atoms + +class SO3: + ''' + A class to generate the SO3 power spectrum components + based off of the Gaussian atomic neighbor density function + defined in "On Representing Atomic Environments". + + args: + nmax: int, degree of radial expansion + lmax: int, degree of spherical harmonic expansion + rcut: float, cutoff radius for neighbor calculation + alpha: float, gaussian width parameter + weight_on: bool, if True, the neighbors with different type will be counted as negative + ''' + + def __init__(self, nmax=3, lmax=3, rcut=3.5, alpha=2.0, + weight_on=False, neighborlist='ase'): + # populate attributes + self.nmax = nmax + self.lmax = lmax + self.rcut = rcut + self.alpha = alpha + self._type = "SO3" + self.cutoff_function = 'cosine' + self.weight_on = weight_on + self.neighborcalc = neighborlist + #return + + def __str__(self): + s = "SO3 descriptor with Cutoff: {:6.3f}".format(self.rcut) + s += " lmax: {:d}, nmax: {:d}, alpha: {:.3f}\n".format(self.lmax, self.nmax, self.alpha) + s += "neighborlist: {:s}\n".format(self.neighborcalc) + return s + + def __repr__(self): + return str(self) + + def load_from_dict(self, dict0): + self.nmax = dict0["nmax"] + self.lmax = dict0["lmax"] + self.rcut = dict0["rcut"] + self.alpha = dict0["alpha"] + self.derivative = dict0["derivative"] + + def save_dict(self): + """ + save the model as a dictionary in json + """ + dict = {"nmax": self.nmax, + "lmax": self.lmax, + "rcut": self.rcut, + "alpha": self.alpha, + "derivative": self.derivative, + "_type": "SO3", + } + return dict + + @property + def nmax(self): + return self._nmax + + @nmax.setter + def nmax(self, nmax): + if isinstance(nmax, int) is True: + if nmax < 1: + raise ValueError('nmax must be greater than or equal to 1') + if nmax > 11: + raise ValueError('nmax > 11 yields complex eigenvalues which will mess up the calculation') + self._nmax = nmax + else: + raise ValueError('nmax must be an integer') + + @property + def lmax(self): + return self._lmax + + @lmax.setter + def lmax(self, lmax): + if isinstance(lmax, int) is True: + if lmax < 0: + raise ValueError('lmax must be greater than or equal to zero') + elif lmax > 32: + raise NotImplementedError('''Currently we only support Wigner-D matrices and spherical harmonics + for arguments up to l=32. If you need higher functionality, raise an issue + in our Github and we will expand the set of supported functions''') + self._lmax = lmax + else: + raise ValueError('lmax must be an integer') + + @property + def rcut(self): + return self._rcut + + @rcut.setter + def rcut(self, rcut): + if isinstance(rcut, float) or isinstance(rcut, int): + if rcut <= 0: + raise ValueError('rcut must be greater than zero') + self._rcut = rcut + else: + raise ValueError('rcut must be a float') + + @property + def alpha(self): + return self._alpha + + @alpha.setter + def alpha(self, alpha): + if isinstance(alpha, float) or isinstance(alpha, int): + if alpha <= 0: + raise ValueError('alpha must be greater than zero') + self._alpha = alpha + else: + raise ValueError('alpha must be a float') + + @property + def derivative(self): + return self._derivative + + @derivative.setter + def derivative(self, derivative): + if isinstance(derivative, bool) is True: + self._derivative = derivative + else: + raise ValueError('derivative must be a boolean value') + + @property + def cutoff_function(self): + return self._cutoff_function + + @cutoff_function.setter + def cutoff_function(self, cutoff_function): + self._cutoff_function = Cosine + + def clear_memory(self): + ''' + Clears all non essential attributes for the calculator + ''' + attrs = list(vars(self).keys()) + for attr in attrs: + if attr not in {'_nmax', '_lmax', '_rcut', '_alpha', '_derivative', '_cutoff_function', 'weight_on', 'neighborcalc'}: + delattr(self, attr) + return + + def calculate(self, atoms, atom_ids=None, derivative=False): + ''' + Calculates the SO(3) power spectrum components of the + smoothened atomic neighbor density function + for given nmax, lmax, rcut, and alpha. + + Args: + atoms: an ASE atoms object corresponding to the desired + atomic arrangement + atom_ids: + derivative: bool, whether to calculate the gradient of not + ''' + self._atoms = atoms + self.build_neighbor_list(atom_ids) + self.initialize_arrays() + + ncoefs = self.nmax*(self.nmax+1)//2*(self.lmax+1) + tril_indices = np.tril_indices(self.nmax, k=0) + + ls = np.arange(self.lmax+1) + norm = np.sqrt(2*np.sqrt(2)*np.pi/np.sqrt(2*ls+1)) + + if derivative: + # get expansion coefficients and derivatives + cs, dcs = compute_dcs(self.neighborlist, self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) + + # weight cs and dcs + cs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis] + dcs *= self.atomic_weights[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis] + cs = np.einsum('inlm,l->inlm', cs, norm) + dcs = np.einsum('inlmj,l->inlmj', dcs, norm) + #print('cs, dcs', self.neighbor_indices, cs.shape, dcs.shape) + + # Assign cs and dcs to P and dP + # cs: (N_ij, n, l, m) => P (N_i, N_des) + # dcs: (N_ij, n, l, m, 3) => dP (N_i, N_j, N_des, 3) + # (n, l, m) needs to be merged to 1 dimension + + for i in range(len(atoms)): + # find atoms for which i is the center + centers = self.neighbor_indices[:, 0] == i + + # total up the c array for the center atom + ctot = cs[centers].sum(axis=0) #(n, l, m) + + # power spectrum P = c*c_conj + # eq_3 (n, n', l) eliminate m + P = np.einsum('ijk, ljk->ilj', ctot, np.conj(ctot)).real + + # merge (n, n', l) to 1 dimension + self._plist[i] = P[tril_indices].flatten() + + # gradient of P for each neighbor, eq_26 + # (N_ijs, n, n', l, 3) + # dc * c_conj + c * dc_conj + dP = np.einsum('wijkn,ljk->wiljn', dcs[centers], np.conj(ctot)) + dP += np.conj(np.transpose(dP, axes=[0,2,1,3,4])) + dP = dP.real + + #print("shape of P/dP", P.shape, dP.shape)#; import sys; sys.exit() + + #ijs = self.neighbor_indices[centers] + #for _id, j in enumerate(ijs[:, 1]): + # self._dplist[i, j, :, :] += dP[_id][tril_indices].flatten().reshape(ncoefs, 3) + # # QZ: to check + # self._dplist[i, i, :, :] += dP[_id][tril_indices].flatten().reshape(ncoefs, 3) + + ijs = self.neighbor_indices[centers] + for _id, (i_idx, j_idx) in enumerate(ijs):#(ijs[:, 1]): + Rij = atoms.positions[j_idx] - atoms.positions[i_idx] + norm_Rij = np.linalg.norm(Rij) + for m in range(len(atoms)): + if m != i_idx and m != j_idx: + normalization_factor = 0 + self._dplist[i, m, :, :] += dP[_id][tril_indices].flatten().reshape(ncoefs, 3) * normalization_factor + elif m == i_idx: + normalization_factor = -1 / norm_Rij + self._dplist[i, m, :, :] += dP[_id][tril_indices].flatten().reshape(ncoefs, 3) * normalization_factor + elif m == j_idx: + normalization_factor = 1 / norm_Rij + self._dplist[i, m, :, :] += dP[_id][tril_indices].flatten().reshape(ncoefs, 3) * normalization_factor + + x = {'x':self._plist, + 'dxdr':self._dplist, + 'elements':list(atoms.symbols)} + else: + if len(self.neighborlist) > 0: + cs = compute_cs(self.neighborlist, self.nmax, self.lmax, self.rcut, self.alpha, self._cutoff_function) + cs *= self.atomic_weights[:,np.newaxis,np.newaxis,np.newaxis] + cs = np.einsum('inlm,l->inlm', cs, norm) + # everything good up to here + for i in range(len(atoms)): + centers = self.neighbor_indices[:,0] == i + ctot = cs[centers].sum(axis=0) + P = np.einsum('ijk,ljk->ilj', ctot, np.conj(ctot)).real + self._plist[i] = P[tril_indices].flatten() + x = {'x': self._plist, + 'dxdr': None, + 'elements': list(atoms.symbols)} + + self.clear_memory() + return x + + def initialize_arrays(self): + # number of atoms + natoms = len(self._atoms) #self._atoms) + + # degree of spherical harmonic expansion + lmax = self.lmax + + # degree of radial expansion + nmax = self.nmax + + # number of unique power spectrum components + # this is given by the triangular elements of + # the radial expansion multiplied by the degree + # of spherical harmonic expansion (including 0) + ncoefs = nmax*(nmax+1)//2*(lmax+1) + + self._plist = np.zeros((natoms, ncoefs), dtype=np.float64) + self._dplist = np.zeros((natoms, natoms, ncoefs, 3), dtype=np.float64) + + return + + def build_neighbor_list(self, atom_ids=None): + ''' + Builds a neighborlist for the calculation of bispectrum components for + a given ASE atoms object given in the calculate method. + ''' + atoms = self._atoms + if atom_ids is None: + atom_ids = range(len(atoms)) + + cutoffs = [self.rcut/2]*len(atoms) + if self.neighborcalc == 'ase': + from ase.neighborlist import NeighborList + nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0) + nl.update(atoms) + else: + from neighborlist import NeighborList + nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0) + nl.update(atoms, atom_ids) + #print(atoms, atom_ids) + #print(atoms.get_scaled_positions()) + + center_atoms = [] + neighbors = [] + neighbor_indices = [] + atomic_weights = [] + temp_indices = [] + + for i in atom_ids: + # get center atom position vector + center_atom = atoms.positions[i] + # get indices and cell offsets for each neighbor + indices, offsets = nl.get_neighbors(i) + #print(indices); import sys; sys.exit() + temp_indices.append(indices) + for j, offset in zip(indices, offsets): + pos = atoms.positions[j] + np.dot(offset,atoms.get_cell()) - center_atom + # to prevent division by zero + if np.sum(np.abs(pos)) < 1e-3: pos += 0.001 + center_atoms.append(center_atom) + neighbors.append(pos) + if self.weight_on and atoms[j].number != atoms[i].number: + factor = -1 + else: + factor = 1 + atomic_weights.append(factor*atoms[j].number) + neighbor_indices.append([i,j]) + + neighbor_indices = np.array(neighbor_indices, dtype=np.int64) + + self.center_atoms = np.array(center_atoms, dtype=np.float64) + self.neighborlist = np.array(neighbors, dtype=np.float64) + self.atomic_weights = np.array(atomic_weights, dtype=np.int64) + self.neighbor_indices = neighbor_indices + return + +def Cosine(Rij, Rc, derivative=False): + # Rij is the norm + if derivative is False: + result = 0.5 * (np.cos(np.pi * Rij / Rc) + 1.) + else: + result = -0.5 * np.pi / Rc * np.sin(np.pi * Rij / Rc) + return result + +def W(nmax): + arr = np.zeros((nmax,nmax), np.float64) + for alpha in range(1, nmax+1, 1): + temp1 = (2*alpha+5)*(2*alpha+6)*(2*alpha+7) + for beta in range(1, alpha+1, 1): + temp2 = (2*beta+5)*(2*beta+6)*(2*beta+7) + arr[alpha-1, beta-1] = np.sqrt(temp1*temp2)/(5+alpha+beta)/(6+alpha+beta)/(7+alpha+beta) + arr[beta-1, alpha-1] = arr[alpha-1, beta-1] + + sinv = np.linalg.inv(arr) + eigvals, V = np.linalg.eig(sinv) + sqrtD = np.diag(np.sqrt(eigvals)) + arr[:,:] = np.dot(np.dot(V, sqrtD), np.linalg.inv(V)) + return arr + +def phi(r, alpha, rcut): + ''' + See g below + ''' + return (rcut-r)**(alpha+2)/np.sqrt(2*rcut**(2*alpha+7)/(2*alpha+5)/(2*alpha+6)/(2*alpha+7)) + +def g(r, n, nmax, rcut, w): + + Sum = 0.0 + for alpha in range(1, nmax+1): + Sum += w[n-1, alpha-1]*phi(r, alpha, rcut) + + return Sum + +def GaussChebyshevQuadrature(nmax, lmax): + NQuad = (nmax+lmax+1) * 10 # + quad_array = np.zeros(NQuad, dtype=np.float64) + for i in range(1, NQuad+1): + # roots of Chebyshev polynomial of degree N + x = np.cos((2*i-1)*np.pi/2/NQuad) + quad_array[i-1] = x + return quad_array, np.pi/NQuad + +def compute_cs(pos, nmax, lmax, rcut, alpha, cutoff): + """ + Compute exapnsion coefficients + + Args: + pos: + nmax (int): + lmax (int): + rcut (float): + alpha (float): + cutoff (callable): + + Returns: + C(N_ij, nmax, lmax+1, 2lmax+1) + """ + # compute the overlap matrix + w = W(nmax) + + # get the norm of the position vectors + Ris = np.linalg.norm(pos, axis=1) # (Nneighbors) + + # initialize Gauss Chebyshev Quadrature + GCQuadrature, weight = GaussChebyshevQuadrature(nmax, lmax) #(Nquad) + weight *= rcut/2 + # transform the quadrature from (-1,1) to (0, rcut) + Quadrature = rcut/2*(GCQuadrature+1) + + # compute the arguments for the bessel functions + BesselArgs = 2*alpha*np.outer(Ris,Quadrature)#(Nneighbors x Nquad) + + # initalize the arrays for the bessel function values + # and the G function values + Bessels = np.zeros((len(Ris), len(Quadrature), lmax+1), dtype=np.float64) #(Nneighbors x Nquad x lmax+1) + Gs = np.zeros((nmax, len(Quadrature)), dtype=np.float64) # (nmax, nquad) + + # compute the g values + for n in range(1,nmax+1,1): + Gs[n-1,:] = g(Quadrature, n, nmax, rcut, w) + + # compute the bessel values + for l in range(lmax+1): + Bessels[:,:,l] = spherical_in(l, BesselArgs) + + # mutliply the terms in the integral separate from the Bessels + Quad_Squared = Quadrature**2 + Gs *= Quad_Squared * np.exp(-alpha*Quad_Squared) * np.sqrt(1-GCQuadrature**2) * weight + + # perform the integration with the Bessels + integral_array = np.einsum('ij,kjl->kil', Gs, Bessels) # (Nneighbors x nmax x lmax+1) + + # compute the gaussian for each atom and multiply with 4*pi + # to minimize floating point operations + # weight can also go here since the Chebyshev gauss quadrature weights are uniform + exparray = 4*np.pi*np.exp(-alpha*Ris**2) # (Nneighbors) + + cutoff_array = cutoff(Ris, rcut) + + exparray *= cutoff_array + + # get the spherical coordinates of each atom + thetas = np.arccos(pos[:,2]/Ris[:]) + phis = np.arctan2(pos[:,1], pos[:,0]) + + # determine the size of the m axis + msize = 2*lmax+1 + # initialize an array for the spherical harmonics + ylms = np.zeros((len(Ris), lmax+1, msize), dtype=np.complex128) + + # compute the spherical harmonics + for l in range(lmax+1): + for m in range(-l,l+1,1): + midx = msize//2 + m + ylms[:,l,midx] = sph_harm(m, l, phis, thetas) + + # multiply the spherical harmonics and the radial inner product + Y_mul_innerprod = np.einsum('ijk,ilj->iljk', ylms, integral_array) + + # multiply the gaussians into the expression + C = np.einsum('i,ijkl->ijkl', exparray, Y_mul_innerprod) + return C + +def compute_dcs(pos, nmax, lmax, rcut, alpha, cutoff): + """ + Compute exapnsion coefficients + + Args: + pos: + nmax (int): + lmax (int): + rcut (float): + alpha (float): + cutoff (callable): + + Returns: + c(N_ij, nmax, lmax+1, 2lmax+1) + dc(N_ij, nmax, lmax+1, 2lmax+1, 3) for each x,y,z + """ + # compute the overlap matrix + w = W(nmax) + + # get the norm of the position vectors + Ris = np.linalg.norm(pos, axis=1) # (Nneighbors) + + # get unit vectors + upos = pos/Ris[:,np.newaxis] + + # initialize Gauss Chebyshev Quadrature + GCQuadrature, weight = GaussChebyshevQuadrature(nmax,lmax) #(Nquad) + weight *= rcut/2 + # transform from (-1,1) to (0, rcut) + Quadrature = rcut/2*(GCQuadrature+1) + + # compute the arguments for the bessel functions + BesselArgs = 2*alpha*np.outer(Ris,Quadrature)#(Nneighbors x Nquad) + + # initalize the arrays for the bessel function values + # and the G function values + Bessels = np.zeros((len(Ris), len(Quadrature), lmax+1), dtype=np.float64) #(Nneighbors x Nquad x lmax+1) + Gs = np.zeros((nmax, len(Quadrature)), dtype=np.float64) # (nmax, nquad) + dBessels = np.zeros((len(Ris), len(Quadrature), lmax+1), dtype=np.float64) #(Nneighbors x Nquad x lmax+1) + + # compute the g values + for n in range(1,nmax+1,1): + Gs[n-1,:] = g(Quadrature, n, nmax, rcut,w)*weight + + # compute the bessel values + for l in range(lmax+1): + Bessels[:,:,l] = spherical_in(l, BesselArgs) + dBessels[:,:,l] = spherical_in(l, BesselArgs, derivative=True) + + #(Nneighbors x Nquad x lmax+1) unit vector here + gradBessels = np.einsum('ijk,in->ijkn',dBessels,upos) + gradBessels *= 2*alpha + + # multiply with r for the integral + gradBessels = np.einsum('ijkn,j->ijkn',gradBessels,Quadrature) + + # mutliply the terms in the integral separate from the Bessels + Quad_Squared = Quadrature**2 + Gs *= Quad_Squared * np.exp(-alpha*Quad_Squared) * np.sqrt(1-GCQuadrature**2) + + # perform the integration with the Bessels + integral_array = np.einsum('ij,kjl->kil', Gs, Bessels) # (Nneighbors x nmax x lmax+1) + + grad_integral_array = np.einsum('ij,kjlm->kilm', Gs, gradBessels)# (Nneighbors x nmax x lmax+1, 3) + + # compute the gaussian for each atom + exparray = 4*np.pi*np.exp(-alpha*Ris**2) # (Nneighbors) + + gradexparray = (-2*alpha*Ris*exparray)[:,np.newaxis]*upos + + cutoff_array = cutoff(Ris, rcut) + + grad_cutoff_array = np.einsum('i,in->in',cutoff(Ris, rcut, True), upos) + + # get the spherical coordinates of each atom + thetas = np.arccos(pos[:,2]/Ris[:]) + phis = np.arctan2(pos[:,1], pos[:,0]) + + # the size changes temporarily for the derivative + # determine the size of the m axis + Msize = 2*(lmax+1)+1 + msize = 2*lmax + 1 + # initialize an array for the spherical harmonics and gradients + #(Nneighbors, l, m, *3*) + ylms = np.zeros((len(Ris), lmax+1+1, Msize), dtype=np.complex128) + gradylms = np.zeros((len(Ris), lmax+1, msize, 3), dtype=np.complex128) + # compute the spherical harmonics + for l in range(lmax+1+1): + for m in range(-l,l+1,1): + midx = Msize//2 + m + ylms[:,l,midx] = sph_harm(m, l, phis, thetas) + + + for l in range(1, lmax+1): + for m in range(-l, l+1, 1): + midx = msize//2 + m + Midx = Msize//2 + m + # get gradient with recpect to spherical covariant components + xcov0 = -np.sqrt(((l+1)**2-m**2)/(2*l+1)/(2*l+3))*l*ylms[:,l+1,Midx]/Ris + + if abs(m) <= l-1: + xcov0 += np.sqrt((l**2-m**2)/(2*l-1)/(2*l+1))*(l+1)*ylms[:,l-1,Midx]/Ris + + + xcovpl1 = -np.sqrt((l+m+1)*(l+m+2)/2/(2*l+1)/(2*l+3))*l*ylms[:,l+1,Midx+1]/Ris + + if abs(m+1) <= l-1: + xcovpl1 -= np.sqrt((l-m-1)*(l-m)/2/(2*l-1)/(2*l+1))*(l+1)*ylms[:,l-1,Midx+1]/Ris + + + xcovm1 = -np.sqrt((l-m+1)*(l-m+2)/2/(2*l+1)/(2*l+3))*l*ylms[:,l+1,Midx-1]/Ris + + if abs(m-1) <= l-1: + xcovm1 -= np.sqrt((l+m-1)*(l+m)/2/(2*l-1)/(2*l+1))*(l+1)*ylms[:,l-1,Midx-1]/Ris + + #transform the gradient to cartesian + gradylms[:,l,midx,0] = 1/np.sqrt(2)*(xcovm1-xcovpl1) + gradylms[:,l,midx,1] = 1j/np.sqrt(2)*(xcovm1+xcovpl1) + gradylms[:,l,midx,2] = xcov0 + + # index ylms to get rid of extra terms for derivative + ylms = ylms[:, 0:lmax+1, 1:1+2*lmax+1] + # multiply the spherical harmonics and the radial inner product + Y_mul_innerprod = np.einsum('ijk,ilj->iljk', ylms, integral_array) + # multiply the gradient of the spherical harmonics with the radial inner get_radial_inner_product + dY_mul_innerprod = np.einsum('ijkn,ilj->iljkn', gradylms, integral_array) + # multiply the spherical harmonics with the gradient of the radial inner get_radial_inner_product + Y_mul_dinnerprod = np.einsum('ijk,iljn->iljkn', ylms, grad_integral_array) + # multiply the gaussians into the expression with 4pi + C = np.einsum('i,ijkl->ijkl', exparray, Y_mul_innerprod) + # multiply the gradient of the gaussian with the other terms + gradexp_mul_y_inner = np.einsum('in,ijkl->ijkln', gradexparray, Y_mul_innerprod) + # add gradient of inner product and spherical harmonic terms + gradHarmonics_mul_gaussian = np.einsum('ijkln,i->ijkln', dY_mul_innerprod+Y_mul_dinnerprod, exparray) + dC = gradexp_mul_y_inner + gradHarmonics_mul_gaussian + dC *= cutoff_array[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis] + dC += np.einsum('in,ijkl->ijkln', grad_cutoff_array, C) + C *= cutoff_array[:, np.newaxis, np.newaxis, np.newaxis] + return C, dC + +if __name__ == "__main__": + from optparse import OptionParser + from ase.io import read + import time + # ---------------------- Options ------------------------ + parser = OptionParser() + parser.add_option("-c", "--crystal", dest="structure", + help="crystal from file, cif or poscar, REQUIRED", + metavar="crystal") + + parser.add_option("-r", "--rcut", dest="rcut", default=3.0, type=float, + help="cutoff for neighbor calcs, default: 3.0" + ) + + parser.add_option("-l", "--lmax", dest="lmax", default=2, type=int, + help="lmax, default: 1" + ) + + parser.add_option("-n", "--nmax", dest="nmax", default=1, type=int, + help="nmax, default: 1" + ) + + parser.add_option("-a", "--alpha", dest="alpha", default=2.0, type=float, + help="cutoff for neighbor calcs, default: 2.0" + ) + + parser.add_option("-f", dest="der", default=True, + action='store_false',help='derivative flag') + + (options, args) = parser.parse_args() + + if options.structure is None: + from ase.build import bulk + test = bulk('Si', 'diamond', a=5.459, cubic=True) + else: + test = read(options.structure, format='vasp') + + cell = test.get_cell() + cell[0,1] += 0.5 + test.set_cell(cell) + + lmax = options.lmax + nmax = options.nmax + rcut = options.rcut + alpha = options.alpha + der = options.der + + start1 = time.time() + f = SO3(nmax=nmax, lmax=lmax, rcut=rcut, alpha=alpha, cutoff_function='cosine') + x = f.calculate(test, derivative=True) + start2 = time.time() + print('x', x['x']) + print('dxdr', x['dxdr']) + print('calculation time {}'.format(start2-start1)) diff --git a/pyxtal/lego/__init__.py b/pyxtal/lego/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pyxtal/lego/basinhopping.py b/pyxtal/lego/basinhopping.py new file mode 100644 index 00000000..0f3f5700 --- /dev/null +++ b/pyxtal/lego/basinhopping.py @@ -0,0 +1,786 @@ +""" +basinhopping: The basinhopping global optimization algorithm +""" +import numpy as np +import math +import inspect +import scipy.optimize +from scipy._lib._util import check_random_state +from copy import deepcopy + +__all__ = ['basinhopping'] + + +_params = (inspect.Parameter('res_new', kind=inspect.Parameter.KEYWORD_ONLY), + inspect.Parameter('res_old', kind=inspect.Parameter.KEYWORD_ONLY)) +_new_accept_test_signature = inspect.Signature(parameters=_params) + + +class Storage: + """ + Class used to store the lowest energy structure + """ + def __init__(self, minres): + self._add(minres) + + def _add(self, minres): + self.minres = minres + self.minres.x = np.copy(minres.x) + + def update(self, minres): + #if minres.success and (minres.fun < self.minres.fun + # or not self.minres.success): + if minres.fun < self.minres.fun: + self._add(minres) + return True + else: + return False + + def get_lowest(self): + return self.minres + + +class BasinHoppingRunner: + """This class implements the core of the basinhopping algorithm. + + x0 : ndarray + The starting coordinates. + minimizer : callable + The local minimizer, with signature ``result = minimizer(x)``. + The return value is an `optimize.OptimizeResult` object. + step_taking : callable + This function displaces the coordinates randomly. Signature should + be ``x_new = step_taking(x)``. Note that `x` may be modified in-place. + accept_tests : list of callables + Each test is passed the kwargs `f_new`, `x_new`, `f_old` and + `x_old`. These tests will be used to judge whether or not to accept + the step. The acceptable return values are True, False, or ``"force + accept"``. If any of the tests return False then the step is rejected. + If ``"force accept"``, then this will override any other tests in + order to accept the step. This can be used, for example, to forcefully + escape from a local minimum that ``basinhopping`` is trapped in. + disp : bool, optional + Display status messages. + + """ + def __init__(self, x0, minimizer, step_taking, accept_tests, disp=False): + self.x = np.copy(x0) + self.minimizer = minimizer + self.step_taking = step_taking + self.accept_tests = accept_tests + self.disp = disp + self.xtrials = [] + self.energy_trials = [] + self.nstep = 0 + + # initialize return object + self.res = scipy.optimize.OptimizeResult() + self.res.minimization_failures = 0 + + # do initial minimization + minres = minimizer(self.x) + if not minres.success: + self.res.minimization_failures += 1 + if self.disp: + print("warning: basinhopping: local minimization failure") + self.x = np.copy(minres.x) + self.energy = minres.fun + self.incumbent_minres = minres # best minimize result found so far + if self.disp: + print("basinhopping step %d: f %g" % (self.nstep, self.energy)) + + self.xtrials.append(self.x) + self.energy_trials.append(self.energy) + + # initialize storage class + self.storage = Storage(minres) + + if hasattr(minres, "nfev"): + self.res.nfev = minres.nfev + if hasattr(minres, "njev"): + self.res.njev = minres.njev + if hasattr(minres, "nhev"): + self.res.nhev = minres.nhev + + def _monte_carlo_step(self): + """Do one Monte Carlo iteration + + Randomly displace the coordinates, minimize, and decide whether + or not to accept the new coordinates. + """ + # Take a random step. Make a copy of x because the step_taking + # algorithm might change x in place + x_after_step = np.copy(self.x) + x_after_step = self.step_taking(x_after_step) + + # do a local minimization + minres = self.minimizer(x_after_step) + #print('====================================', minres.fun) + x_after_quench = minres.x + energy_after_quench = minres.fun + if not minres.success: + self.res.minimization_failures += 1 + if self.disp: + print("warning: basinhopping: local minimization failure") + if hasattr(minres, "nfev"): + self.res.nfev += minres.nfev + if hasattr(minres, "njev"): + self.res.njev += minres.njev + if hasattr(minres, "nhev"): + self.res.nhev += minres.nhev + + # accept the move based on self.accept_tests. If any test is False, + # then reject the step. If any test returns the special string + # 'force accept', then accept the step regardless. This can be used + # to forcefully escape from a local minimum if normal basin hopping + # steps are not sufficient. + accept = True + for test in self.accept_tests: + if inspect.signature(test) == _new_accept_test_signature: + testres = test(res_new=minres, res_old=self.incumbent_minres) + else: + testres = test(f_new=energy_after_quench, x_new=x_after_quench, + f_old=self.energy, x_old=self.x) + + if testres == 'force accept': + accept = True + break + elif testres is None: + raise ValueError("accept_tests must return True, False, or " + "'force accept'") + elif not testres: + accept = False + + + # Report the result of the acceptance test to the take step class. + # This is for adaptive step taking + if hasattr(self.step_taking, "report"): + self.step_taking.report(accept, f_new=energy_after_quench, + x_new=x_after_quench, f_old=self.energy, + x_old=self.x) + + if minres.fun < self.incumbent_minres.fun and not accept: import sys; sys.exit() + return accept, minres + + def one_cycle(self): + """Do one cycle of the basinhopping algorithm + """ + self.nstep += 1 + new_global_min = False + + accept, minres = self._monte_carlo_step() + + if accept: + self.energy = minres.fun + self.x = np.copy(minres.x) + self.incumbent_minres = minres # best minimize result found so far + new_global_min = self.storage.update(minres) + #print('========================================================update', minres.fun, self.storage.minres.fun) + + # print some information + if self.disp: + self.print_report(minres.fun, accept) + if new_global_min: + print("found new global minimum on step %d with function" + " value %g" % (self.nstep, self.energy)) + + # save some variables as BasinHoppingRunner attributes + self.xtrial = minres.x + self.energy_trial = minres.fun + self.accept = accept + self.xtrials.append(minres.x) + self.energy_trials.append(minres.fun) + + return new_global_min + + def print_report(self, energy_trial, accept): + """print a status update""" + minres = self.storage.get_lowest() + print("basinhopping step %d: f %g trial_f %g accepted %d " + " lowest_f %g" % (self.nstep, self.energy, energy_trial, + accept, minres.fun)) + + +class AdaptiveStepsize: + """ + Class to implement adaptive stepsize. + + This class wraps the step taking class and modifies the stepsize to + ensure the true acceptance rate is as close as possible to the target. + + Parameters + ---------- + takestep : callable + The step taking routine. Must contain modifiable attribute + takestep.stepsize + accept_rate : float, optional + The target step acceptance rate + interval : int, optional + Interval for how often to update the stepsize + factor : float, optional + The step size is multiplied or divided by this factor upon each + update. + verbose : bool, optional + Print information about each update + + """ + def __init__(self, takestep, accept_rate=0.5, interval=50, factor=0.9, + verbose=True): + self.takestep = takestep + self.target_accept_rate = accept_rate + self.interval = interval + self.factor = factor + self.verbose = verbose + + self.nstep = 0 + self.nstep_tot = 0 + self.naccept = 0 + + def __call__(self, x): + return self.take_step(x) + + def _adjust_step_size(self): + old_stepsize = self.takestep.stepsize + accept_rate = float(self.naccept) / self.nstep + if accept_rate > self.target_accept_rate: + # We're accepting too many steps. This generally means we're + # trapped in a basin. Take bigger steps. + self.takestep.stepsize /= self.factor + else: + # We're not accepting enough steps. Take smaller steps. + self.takestep.stepsize *= self.factor + if self.verbose: + print("adaptive stepsize: acceptance rate {:f} target {:f} new " + "stepsize {:g} old stepsize {:g}".format(accept_rate, + self.target_accept_rate, self.takestep.stepsize, + old_stepsize)) + + def take_step(self, x): + self.nstep += 1 + self.nstep_tot += 1 + if self.nstep % self.interval == 0: + self._adjust_step_size() + return self.takestep(x) + + def report(self, accept, **kwargs): + "called by basinhopping to report the result of the step" + if accept: + self.naccept += 1 + + +class RandomDisplacement: + """Add a random displacement of maximum size `stepsize` to each coordinate. + + Calling this updates `x` in-place. + + Parameters + ---------- + stepsize : float, optional + Maximum stepsize in any dimension + random_gen : {None, int, `numpy.random.Generator`, + `numpy.random.RandomState`}, optional + + If `seed` is None (or `np.random`), the `numpy.random.RandomState` + singleton is used. + If `seed` is an int, a new ``RandomState`` instance is used, + seeded with `seed`. + If `seed` is already a ``Generator`` or ``RandomState`` instance then + that instance is used. + + """ + + def __init__(self, stepsize=0.5, random_gen=None): + self.stepsize = stepsize + self.random_gen = check_random_state(random_gen) + + def __call__(self, x): + x += self.random_gen.uniform(-self.stepsize, self.stepsize, + np.shape(x)) + return x + +class MinimizerWrapper: + """ + wrap a minimizer function as a minimizer class + """ + + def __init__(self, minimizer, func=None, **kwargs): + self.minimizer = minimizer + self.func = func + self.kwargs = kwargs + + def __call__(self, x0): + if self.func is None: + return self.minimizer(x0, **self.kwargs) + else: + if 'method' in self.kwargs and type(self.kwargs['method']) is list: + x = x0 + for i, method in enumerate(self.kwargs['method']): + kwargs = deepcopy(self.kwargs) + kwargs['method'] = method + if method == 'Nelder-Mead': + kwargs['options'].pop('ftol') + else: + kwargs['options'].pop('fatol') + res = self.minimizer(self.func, x, **kwargs) + if i + 1 < len(self.kwargs['method']): + x = res.x + #print(method, res.fun) + return res + else: + return self.minimizer(self.func, x0, **self.kwargs) + + +class Metropolis: + """Metropolis acceptance criterion. + + Parameters + ---------- + T : float + The "temperature" parameter for the accept or reject criterion. + random_gen : {None, int, `numpy.random.Generator`, + `numpy.random.RandomState`}, optional + + If `seed` is None (or `np.random`), the `numpy.random.RandomState` + singleton is used. + If `seed` is an int, a new ``RandomState`` instance is used, + seeded with `seed`. + If `seed` is already a ``Generator`` or ``RandomState`` instance then + that instance is used. + Random number generator used for acceptance test. + + """ + + def __init__(self, T, random_gen=None): + # Avoid ZeroDivisionError since "MBH can be regarded as a special case + # of the BH framework with the Metropolis criterion, where temperature + # T = 0." (Reject all steps that increase energy.) + self.beta = 1.0 / T if T != 0 else float('inf') + self.random_gen = check_random_state(random_gen) + + def accept_reject(self, res_new, res_old): + """ + Assuming the local search underlying res_new was successful: + If new energy is lower than old, it will always be accepted. + If new is higher than old, there is a chance it will be accepted, + less likely for larger differences. + """ + with np.errstate(invalid='ignore'): + # The energy values being fed to Metropolis are 1-length arrays, and if + # they are equal, their difference is 0, which gets multiplied by beta, + # which is inf, and array([0]) * float('inf') causes + # + # RuntimeWarning: invalid value encountered in multiply + # + # Ignore this warning so when the algorithm is on a flat plane, it always + # accepts the step, to try to move off the plane. + prod = -(res_new.fun - res_old.fun) * self.beta + w = math.exp(min(0, prod)) + #print("=================================================================", res_new.fun, res_old.fun, w) + + rand = self.random_gen.uniform() + return w >= rand #and (res_new.success or not res_old.success) + + def __call__(self, *, res_new, res_old): + """ + f_new and f_old are mandatory in kwargs + """ + return bool(self.accept_reject(res_new, res_old)) + + +def basinhopping(func, x0, niter=100, T=1.0, stepsize=0.5, + minimizer_kwargs=None, take_step=None, accept_test=None, + callback=None, interval=50, disp=False, niter_success=None, + seed=None, *, target_accept_rate=0.5, stepwise_factor=0.9): + """Find the global minimum of a function using the basin-hopping algorithm. + + Basin-hopping is a two-phase method that combines a global stepping + algorithm with local minimization at each step. Designed to mimic + the natural process of energy minimization of clusters of atoms, it works + well for similar problems with "funnel-like, but rugged" energy landscapes + [5]_. + + As the step-taking, step acceptance, and minimization methods are all + customizable, this function can also be used to implement other two-phase + methods. + + Parameters + ---------- + func : callable ``f(x, *args)`` + Function to be optimized. ``args`` can be passed as an optional item + in the dict `minimizer_kwargs` + x0 : array_like + Initial guess. + niter : integer, optional + The number of basin-hopping iterations. There will be a total of + ``niter + 1`` runs of the local minimizer. + T : float, optional + The "temperature" parameter for the acceptance or rejection criterion. + Higher "temperatures" mean that larger jumps in function value will be + accepted. For best results `T` should be comparable to the + separation (in function value) between local minima. + stepsize : float, optional + Maximum step size for use in the random displacement. + minimizer_kwargs : dict, optional + Extra keyword arguments to be passed to the local minimizer + `scipy.optimize.minimize` Some important options could be: + + method : str + The minimization method (e.g. ``"L-BFGS-B"``) + args : tuple + Extra arguments passed to the objective function (`func`) and + its derivatives (Jacobian, Hessian). + + take_step : callable ``take_step(x)``, optional + Replace the default step-taking routine with this routine. The default + step-taking routine is a random displacement of the coordinates, but + other step-taking algorithms may be better for some systems. + `take_step` can optionally have the attribute ``take_step.stepsize``. + If this attribute exists, then `basinhopping` will adjust + ``take_step.stepsize`` in order to try to optimize the global minimum + search. + accept_test : callable, ``accept_test(f_new=f_new, x_new=x_new, f_old=fold, x_old=x_old)``, optional + Define a test which will be used to judge whether to accept the + step. This will be used in addition to the Metropolis test based on + "temperature" `T`. The acceptable return values are True, + False, or ``"force accept"``. If any of the tests return False + then the step is rejected. If the latter, then this will override any + other tests in order to accept the step. This can be used, for example, + to forcefully escape from a local minimum that `basinhopping` is + trapped in. + callback : callable, ``callback(x, f, accept)``, optional + A callback function which will be called for all minima found. ``x`` + and ``f`` are the coordinates and function value of the trial minimum, + and ``accept`` is whether that minimum was accepted. This can + be used, for example, to save the lowest N minima found. Also, + `callback` can be used to specify a user defined stop criterion by + optionally returning True to stop the `basinhopping` routine. + interval : integer, optional + interval for how often to update the `stepsize` + disp : bool, optional + Set to True to print status messages + niter_success : integer, optional + Stop the run if the global minimum candidate remains the same for this + number of iterations. + seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional + + If `seed` is None (or `np.random`), the `numpy.random.RandomState` + singleton is used. + If `seed` is an int, a new ``RandomState`` instance is used, + seeded with `seed`. + If `seed` is already a ``Generator`` or ``RandomState`` instance then + that instance is used. + Specify `seed` for repeatable minimizations. The random numbers + generated with this seed only affect the default Metropolis + `accept_test` and the default `take_step`. If you supply your own + `take_step` and `accept_test`, and these functions use random + number generation, then those functions are responsible for the state + of their random number generator. + target_accept_rate : float, optional + The target acceptance rate that is used to adjust the `stepsize`. + If the current acceptance rate is greater than the target, + then the `stepsize` is increased. Otherwise, it is decreased. + Range is (0, 1). Default is 0.5. + + .. versionadded:: 1.8.0 + + stepwise_factor : float, optional + The `stepsize` is multiplied or divided by this stepwise factor upon + each update. Range is (0, 1). Default is 0.9. + + .. versionadded:: 1.8.0 + + Returns + ------- + res : OptimizeResult + The optimization result represented as a `OptimizeResult` object. + Important attributes are: ``x`` the solution array, ``fun`` the value + of the function at the solution, and ``message`` which describes the + cause of the termination. The ``OptimizeResult`` object returned by the + selected minimizer at the lowest minimum is also contained within this + object and can be accessed through the ``lowest_optimization_result`` + attribute. See `OptimizeResult` for a description of other attributes. + + See Also + -------- + minimize : + The local minimization function called once for each basinhopping step. + `minimizer_kwargs` is passed to this routine. + + Notes + ----- + Basin-hopping is a stochastic algorithm which attempts to find the global + minimum of a smooth scalar function of one or more variables [1]_ [2]_ [3]_ + [4]_. The algorithm in its current form was described by David Wales and + Jonathan Doye [2]_ http://www-wales.ch.cam.ac.uk/. + + The algorithm is iterative with each cycle composed of the following + features + + 1) random perturbation of the coordinates + + 2) local minimization + + 3) accept or reject the new coordinates based on the minimized function + value + + The acceptance test used here is the Metropolis criterion of standard Monte + Carlo algorithms, although there are many other possibilities [3]_. + + This global minimization method has been shown to be extremely efficient + for a wide variety of problems in physics and chemistry. It is + particularly useful when the function has many minima separated by large + barriers. See the `Cambridge Cluster Database + `_ for databases of molecular + systems that have been optimized primarily using basin-hopping. This + database includes minimization problems exceeding 300 degrees of freedom. + + See the free software program `GMIN `_ + for a Fortran implementation of basin-hopping. This implementation has many + variations of the procedure described above, including more + advanced step taking algorithms and alternate acceptance criterion. + + For stochastic global optimization there is no way to determine if the true + global minimum has actually been found. Instead, as a consistency check, + the algorithm can be run from a number of different random starting points + to ensure the lowest minimum found in each example has converged to the + global minimum. For this reason, `basinhopping` will by default simply + run for the number of iterations `niter` and return the lowest minimum + found. It is left to the user to ensure that this is in fact the global + minimum. + + Choosing `stepsize`: This is a crucial parameter in `basinhopping` and + depends on the problem being solved. The step is chosen uniformly in the + region from x0-stepsize to x0+stepsize, in each dimension. Ideally, it + should be comparable to the typical separation (in argument values) between + local minima of the function being optimized. `basinhopping` will, by + default, adjust `stepsize` to find an optimal value, but this may take + many iterations. You will get quicker results if you set a sensible + initial value for ``stepsize``. + + Choosing `T`: The parameter `T` is the "temperature" used in the + Metropolis criterion. Basinhopping steps are always accepted if + ``func(xnew) < func(xold)``. Otherwise, they are accepted with + probability:: + + exp( -(func(xnew) - func(xold)) / T ) + + So, for best results, `T` should to be comparable to the typical + difference (in function values) between local minima. (The height of + "walls" between local minima is irrelevant.) + + If `T` is 0, the algorithm becomes Monotonic Basin-Hopping, in which all + steps that increase energy are rejected. + + .. versionadded:: 0.12.0 + + References + ---------- + .. [1] Wales, David J. 2003, Energy Landscapes, Cambridge University Press, + Cambridge, UK. + .. [2] Wales, D J, and Doye J P K, Global Optimization by Basin-Hopping and + the Lowest Energy Structures of Lennard-Jones Clusters Containing up to + 110 Atoms. Journal of Physical Chemistry A, 1997, 101, 5111. + .. [3] Li, Z. and Scheraga, H. A., Monte Carlo-minimization approach to the + multiple-minima problem in protein folding, Proc. Natl. Acad. Sci. USA, + 1987, 84, 6611. + .. [4] Wales, D. J. and Scheraga, H. A., Global optimization of clusters, + crystals, and biomolecules, Science, 1999, 285, 1368. + .. [5] Olson, B., Hashmi, I., Molloy, K., and Shehu1, A., Basin Hopping as + a General and Versatile Optimization Framework for the Characterization + of Biological Macromolecules, Advances in Artificial Intelligence, + Volume 2012 (2012), Article ID 674832, :doi:`10.1155/2012/674832` + + Examples + -------- + The following example is a 1-D minimization problem, with many + local minima superimposed on a parabola. + + >>> import numpy as np + >>> from scipy.optimize import basinhopping + >>> func = lambda x: np.cos(14.5 * x - 0.3) + (x + 0.2) * x + >>> x0 = [1.] + + Basinhopping, internally, uses a local minimization algorithm. We will use + the parameter `minimizer_kwargs` to tell basinhopping which algorithm to + use and how to set up that minimizer. This parameter will be passed to + `scipy.optimize.minimize`. + + >>> minimizer_kwargs = {"method": "BFGS"} + >>> ret = basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs, + ... niter=200) + >>> print("global minimum: x = %.4f, f(x) = %.4f" % (ret.x, ret.fun)) + global minimum: x = -0.1951, f(x) = -1.0009 + + Next consider a 2-D minimization problem. Also, this time, we + will use gradient information to significantly speed up the search. + + >>> def func2d(x): + ... f = np.cos(14.5 * x[0] - 0.3) + (x[1] + 0.2) * x[1] + (x[0] + + ... 0.2) * x[0] + ... df = np.zeros(2) + ... df[0] = -14.5 * np.sin(14.5 * x[0] - 0.3) + 2. * x[0] + 0.2 + ... df[1] = 2. * x[1] + 0.2 + ... return f, df + + We'll also use a different local minimization algorithm. Also, we must tell + the minimizer that our function returns both energy and gradient (Jacobian). + + >>> minimizer_kwargs = {"method":"L-BFGS-B", "jac":True} + >>> x0 = [1.0, 1.0] + >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs, + ... niter=200) + >>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0], + ... ret.x[1], + ... ret.fun)) + global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109 + + Here is an example using a custom step-taking routine. Imagine you want + the first coordinate to take larger steps than the rest of the coordinates. + This can be implemented like so: + + >>> class MyTakeStep: + ... def __init__(self, stepsize=0.5): + ... self.stepsize = stepsize + ... self.rng = np.random.default_rng() + ... def __call__(self, x): + ... s = self.stepsize + ... x[0] += self.rng.uniform(-2.*s, 2.*s) + ... x[1:] += self.rng.uniform(-s, s, x[1:].shape) + ... return x + + Since ``MyTakeStep.stepsize`` exists basinhopping will adjust the magnitude + of `stepsize` to optimize the search. We'll use the same 2-D function as + before + + >>> mytakestep = MyTakeStep() + >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs, + ... niter=200, take_step=mytakestep) + >>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0], + ... ret.x[1], + ... ret.fun)) + global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109 + + Now, let's do an example using a custom callback function which prints the + value of every minimum found + + >>> def print_fun(x, f, accepted): + ... print("at minimum %.4f accepted %d" % (f, int(accepted))) + + We'll run it for only 10 basinhopping steps this time. + + >>> rng = np.random.default_rng() + >>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs, + ... niter=10, callback=print_fun, seed=rng) + at minimum 0.4159 accepted 1 + at minimum -0.4317 accepted 1 + at minimum -1.0109 accepted 1 + at minimum -0.9073 accepted 1 + at minimum -0.4317 accepted 0 + at minimum -0.1021 accepted 1 + at minimum -0.7425 accepted 1 + at minimum -0.9073 accepted 1 + at minimum -0.4317 accepted 0 + at minimum -0.7425 accepted 1 + at minimum -0.9073 accepted 1 + + The minimum at -1.0109 is actually the global minimum, found already on the + 8th iteration. + + """ + if target_accept_rate <= 0. or target_accept_rate >= 1.: + raise ValueError('target_accept_rate has to be in range (0, 1)') + if stepwise_factor <= 0. or stepwise_factor >= 1.: + raise ValueError('stepwise_factor has to be in range (0, 1)') + + x0 = np.array(x0) + + # set up the np.random generator + rng = check_random_state(seed) + + # set up minimizer + if minimizer_kwargs is None: + minimizer_kwargs = dict() + wrapped_minimizer = MinimizerWrapper(scipy.optimize.minimize, func, + **minimizer_kwargs) + + # set up step-taking algorithm + if take_step is not None: + if not callable(take_step): + raise TypeError("take_step must be callable") + # if take_step.stepsize exists then use AdaptiveStepsize to control + # take_step.stepsize + if hasattr(take_step, "stepsize"): + take_step_wrapped = AdaptiveStepsize( + take_step, interval=interval, + accept_rate=target_accept_rate, + factor=stepwise_factor, + verbose=disp) + else: + take_step_wrapped = take_step + else: + # use default + displace = RandomDisplacement(stepsize=stepsize, random_gen=rng) + take_step_wrapped = AdaptiveStepsize(displace, interval=interval, + accept_rate=target_accept_rate, + factor=stepwise_factor, + verbose=disp) + + # set up accept tests + accept_tests = [] + if accept_test is not None: + if not callable(accept_test): + raise TypeError("accept_test must be callable") + accept_tests = [accept_test] + + # use default + metropolis = Metropolis(T, random_gen=rng) + accept_tests.append(metropolis) + + if niter_success is None: + niter_success = niter + 2 + + bh = BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped, + accept_tests, disp=disp) + + # The wrapped minimizer is called once during construction of + # BasinHoppingRunner, so run the callback + if callable(callback): + val = callback(bh.storage.minres.x, bh.storage.minres.fun, True) + # no bh search if the target is already reached + if val is not None: + if val: + niter = 0 + # start main iteration loop + count, i = 0, 0 + message = ["requested number of basinhopping iterations completed" + " successfully"] + for i in range(niter): + new_global_min = bh.one_cycle() + + if callable(callback): + # should we pass a copy of x? + val = callback(bh.xtrial, bh.energy_trial, bh.accept) + if val is not None: + if val: + message = ["callback function requested stop early by" + "returning True"] + break + + count += 1 + if new_global_min: + count = 0 + elif count > niter_success: + message = ["success condition satisfied"] + break + + # prepare return object + res = bh.res + res.lowest_optimization_result = bh.storage.get_lowest() + res.x = np.copy(res.lowest_optimization_result.x) + res.fun = res.lowest_optimization_result.fun + res.message = message + res.nit = i + 1 + res.success = res.lowest_optimization_result.success + res.xtrials = bh.xtrials + res.energy_trials = bh.energy_trials + return res diff --git a/pyxtal/lego/builder.py b/pyxtal/lego/builder.py new file mode 100644 index 00000000..a116ea88 --- /dev/null +++ b/pyxtal/lego/builder.py @@ -0,0 +1,1212 @@ +""" +mof_builder aims to generate crystal structure from the defined +building blocks (e.g., SiO2 with 4-coordined Si and 2-coordinated O) +1. Generate all possible wp combinations +2. For each wp combination, + 2.1 generate structure randomly + 2.2 optimize the geomtry + 2.3 parse the coordination + 2.4 save the qualified structure to ase.db + +Todo: +1. Symmetry support for dSdX +2. add parallel for gulp optimiza +""" + +# Standard Libraries +import os +import numpy as np +import random +from scipy.optimize import minimize +from pyxtal.lego.basinhopping import basinhopping +from pyxtal.lego.SO3 import SO3 + +import multiprocessing as mp +from concurrent.futures import ProcessPoolExecutor + +# Material science Libraries +from ase.io import write +from ase.db import connect +from pyxtal import pyxtal +from pyxtal.util import generate_wp_lib +from pyxtal.symmetry import Group +from pyxtal.lattice import Lattice +from pyxtal.db import database_topology + +# Logging and debugging +import logging +from time import time +#np.set_printoptions(precision=3, suppress=True) + +def generate_wp_lib_par(spgs, composition, num_wp, num_fu, num_dof): + """ + A wrapper to generate the wp list in parallel + """ + + my_spgs, wp_libs = [], [] + for spg in spgs: + wp_lib = generate_wp_lib([spg], composition, num_wp, num_fu, num_dof), + if len(wp_lib) > 0: + wp_libs.append(wp_lib) + my_spgs.append(spg) + return (my_spgs, wp_libs) + +def generate_xtal_par(wp_libs, niter, dim, elements, calculator, ref_environments, + criteria, T, N_max, early_quit): + """ + A wrapper to call generate_xtal function in parallel + """ + xtals, sims = [], [] + for wp_lib in wp_libs: + (number, spg, wps, dof) = wp_lib + xtal, sim = generate_xtal(dim, spg, wps, niter*dof, elements, calculator, + ref_environments, criteria, T, N_max, early_quit) + if xtal is not None: + xtals.append(xtal) + sims.append(sim) + + return (xtals, sims) + +def minimize_from_x_par(dim, wp_libs, elements, calculator, ref_environments, + opt_type, T, niter, early_quit, minimizers): + """ + A wrapper to call minimize_from_x function in parallel + """ + xtals = [] + xs = [] + for wp_lib in wp_libs: + (x, spg, wps) = wp_lib + res = minimize_from_x(x, dim, spg, wps, elements, calculator, + ref_environments, T, niter, + early_quit, opt_type, minimizers, + filename=None) + if res is not None: + xtals.append(res[0]) + xs.append(res[1]) + return xtals, xs + +def generate_xtal(dim, spg, wps, niter, elements, calculator, + ref_environments, criteria, T, N_max, early_quit, + dump=False, random_state=None, verbose=False): + """ + Generate a crystal with the desired local environment + + Args: + dim (int): 1, 2, 3 + spg: pyxtal.symmetry.Group object + wps: list of wps for the disired crystal (e.g., [wp1, wp2]) + ref_env: reference enviroment + f: callable function to compute env + n_iter (int): + T (float): for basinhopping + + Returns: + xtal and its similarity + """ + + # Here we start to optimize the xtal based on similarity + + print("\n", dim, spg, wps, T, N_max, early_quit) + count = 0 + while True: + + filename = 'opt.txt' if dump else None + result = minimize_from_x(None, dim, spg, wps, elements, calculator, + ref_environments, T, niter, early_quit, + 'global', filename=filename, + random_state=random_state) + + if result is not None: + (xtal, xs) = result + if xtal.check_validity(criteria, verbose=verbose): + x = xtal.get_1d_rep_x() + sim1 = calculate_S(x, xtal, ref_environments, calculator) + print(xtal.get_xtal_string({'sim': sim1})) + else: + xtal = None + sim1 = None + return xtal, sim1 + + count += 1 + if count == N_max: + break + return None, None + + +def minimize_from_x(x, dim, spg, wps, elements, calculator, ref_environments, + T=0.2, niter=20, early_quit=0.02, opt_type='local', + minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], + filename='local_opt_data.txt', random_state=None): + """ + Generate xtal from the 1d representation + + Args: + x: list of 1D array + spg (int): space group number (1-230) + wps (string): e.g. [['4a', '8b']] + elements (string): e.g., ['Si', 'O'] + """ + g, wps, dof = get_input_from_letters(spg, wps, dim) + l_type = g.lattice_type + sites_wp = [] + sites = [] + numIons = [] + ref_envs = None + for i, wp in enumerate(wps): + site = [] + numIon = 0 + for w in wp: + sites.append((elements[i], w))#.get_label()) + site.append(w.get_label()) + numIon += w.multiplicity + if ref_envs is None: + ref_envs = ref_environments[i] + else: + ref_envs = np.vstack((ref_envs, ref_environments[i])) + sites_wp.append(site) + numIons.append(numIon) + + if len(ref_envs.shape) == 1: + ref_envs = ref_envs.reshape((1, len(ref_envs))) + + xtal = pyxtal() + if x is None: + count = 0 + while True: + count += 1 + try: + xtal.from_random(dim, g, elements, numIons, sites=sites_wp, factor=1.0, random_state=random_state) + except RuntimeError: + print(g.number, numIons, sites) + print("Trouble in generating random xtals from pyxtal, try again") + if xtal.valid: + atoms = xtal.to_ase(resort=False, add_vaccum=False) + try: + des = calculator.calculate(atoms)['x'] + except: + if filename is not None: print('Not a good structure, skip') + continue + x = xtal.get_1d_rep_x() + break + elif count == 5: + return None + else: + sites = [] + for ele, _wps in zip(elements, wps): + for wp in _wps: + sites.append((ele, wp)) + try: + xtal.from_1d_rep(x, sites, dim=dim) + except: + return None + + x0 = x.copy() + # Extract variables, call from Pyxtal + [N_abc, N_ang] = Lattice.get_dofs(xtal.lattice.ltype) + rep = xtal.get_1D_representation() + xyzs = rep.x[1:] + + # Set constraints and minimization + bounds = [(1.5, 50)] * (N_abc) + [(30, 150)] * (N_ang) + + # Special treatment in case the random lattice is small + for i in range(N_abc): + if x[i] < 1.5: x[i] = 1.5 + if x[i] > 50.0: x[i] = 50.0 + + for i in range(N_abc, N_abc + N_ang): + if x[i] < 30.0: x[i] = 30.0 + if x[i] > 150.0: x[i] = 150.0 + + for xyz in xyzs: + if len(xyz) > 2: bounds += [(0.0, 1.0)] * len(xyz[2:]) + + if len(x) != len(bounds): print('debug before min', xtal, x, bounds, len(x), len(bounds)) + + sim0 = calculate_S(x, xtal, ref_envs, calculator) + if filename is not None: + with open(filename, 'a+') as f0: + f0.write('\nSpace Group: {:d}\n'.format(xtal.group.number)) + for element, numIon, site in zip(elements, numIons, sites_wp): + strs = 'Element: {:2s} {:4d} '.format(element, numIon) + for s in site: strs += '{:s} '.format(s) + strs += '\n' + f0.write(strs) + # Initial value + strs = 'Init: {:9.3f} '.format(sim0) + for x0 in x: strs += '{:8.4f} '.format(x0) + strs += '\n' + print(strs) + f0.write(strs) + + # Run local minimization + if opt_type == 'local': + # set call back function for debugging + def print_local_fun(x): + f = calculate_S(x, xtal, ref_envs, calculator) + #if filename is not None: + print("{:.4f} ".format(f), x) + with open(filename, 'a+') as f0: + strs = 'Iter: {:9.3f} '.format(f) + for x0 in x: strs += '{:8.4f} '.format(x0) + strs += '\n' + f0.write(strs) + callback = print_local_fun if filename is not None else None + + + for minimizer in minimizers: + #print("Starting", xtal.lattice, xtal.lattice.is_valid_lattice()) + (method, step) = minimizer + if len(x) != len(bounds): print('debug min', xtal, x, bounds, len(x), len(bounds)) + res = minimize(calculate_S, x, + method=method, + args=(xtal, ref_envs, calculator), + bounds=bounds, + options={'maxiter': step}, #'disp': True}, + callback=callback) + x = res.x + if xtal.lattice is None: + return None + + if filename is not None: + with open(filename, 'a+') as f0: f0.write('END\n') + else: + # set call back function for debugging + def print_fun_local(x): + f = calculate_S(x, xtal, ref_envs, calculator) + #print("{:.4f} ".format(f), x) + if f < early_quit: + return True + else: + return None + callback = print_fun_local #if verbose else None + + minimizer_kwargs = {'method': ['Nelder-Mead', 'l-bfgs-b', + 'Nelder-Mead', 'l-bfgs-b'], + 'args': (xtal, ref_envs, calculator), + 'bounds': bounds, + 'callback': callback, + 'options': {'maxiter': 100, + 'fatol': 1e-6, + 'ftol': 1e-6}} + + bounded_step = RandomDisplacementBounds(np.array([b[0] for b in bounds]), + np.array([b[1] for b in bounds]), + id1 = N_abc + N_ang, + id2 = N_abc) + + # set call back function for debugging + def print_fun(x, f, accepted): + if filename is not None: + print("minimum {:.4f}[{:.4f}] accepted {:d} ".format(\ + f, early_quit, int(accepted)), x[:N_abc]) + if f < early_quit: + #print("Return True", True is not None) + return True + else: + return None + callback = print_fun #if verbose else None + + + # Run BH optimization + res = basinhopping(calculate_S, x, T = T, + minimizer_kwargs = minimizer_kwargs, + niter = niter, + take_step = bounded_step, + callback = callback) + if xtal.lattice is None: + return None + + # Extract the optimized xtal + xtal = pyxtal() + xtal.from_1d_rep(res.x, sites, dim=dim) + return xtal, (x0, res.x) + + +def calculate_dSdx(x, xtal, des_ref, f, eps=1e-4, symmetry=True, verbose=False): + """ + Compute dSdx via the chain rule of dSdr*drdx. + This routine will serve as the jacobian for scipy.minimize + + Args: + x (float): input variable array to describe a crystal structure + xtal (instance): pyxtal object + des_ref (array): reference environment + f (callable): function to compute the environment + symmetry (bool): whether or not turn on symmetry + verbose (bool): output more information + """ + xtal.update_from_1d_rep(x) + atoms = xtal.to_ase(resort=False, add_vaccum=False) #* 2 + atoms.set_positions(atoms.positions + 1e-3 * np.random.random([len(atoms), 3])) + ref_pos = atoms.positions + + #results = f.calculate(atoms, ids, derivative=True) + results = f.calculate(atoms, derivative=True) + dPdr = results['dxdr']#; print(dpdr>0) + P = results['x'] + + # Compute dSdr + dSdr = np.einsum("ik, ijkl -> jl", 2*(P - des_ref), dPdr) + + # Compute drdx via numerical func + drdx = np.zeros([len(atoms), 3, len(x)]) + xtal0 = xtal.copy() + for i in range(len(x)): + x0 = x.copy() + x0[i] += eps + xtal0.update_from_1d_rep(x0) + #print(xtal0) + pos = xtal0.to_ase(resort=False, add_vaccum=False).positions + drdx[:, :, i] = (pos - ref_pos)/eps + #print("drdx", i, x0, '\n', drdx[:, :, i]) + + #dPdx = np.einsum('ijkl, klm->ijm', dPdr, drdx) + #dSdx = np.einsum('ij, ijk->k', 2*(P-des_ref), dPdx) + # dSdx = dSdr * drdx + dSdx = np.einsum("ij, ijk -> k", dSdr, drdx) + + return dSdx, dSdr, dPdr + + +def calculate_S(x, xtal, des_ref, f, verbose=False): + """ + Optimization function used for structure generation + + Args: + x (float): input variable array to describe a crystal structure + xtal (instance): pyxtal object + des_ref: reference environment + f (callable): function to compute the enivronment + verbose (boolean): output more information + """ + if xtal.lattice is None: + des = np.zeros(des_ref.shape) + weights = 1 + else: + xtal.update_from_1d_rep(x) #; print(xtal) + if xtal.lattice is not None: + ids = [0] + weights = [] + for site in xtal.atom_sites: + ids.append(site.wp.multiplicity + ids[-1]) + weights.append(site.wp.multiplicity) + ids = ids[:-1] + weights = np.array(weights, dtype=float) + #try: + if True: + atoms = xtal.to_ase(resort=False, add_vaccum=False) + des = f.calculate(atoms, ids)['x'][ids] + #except: + # #print("bug in xtal2ase or bad structure", xtal.lattice) + # des = np.zeros(des_ref.shape) + # #raise ValueError('Skip the random error') + else: + des = np.zeros(des_ref.shape) + weights = 1 + + sim = np.sum((des-des_ref)**2, axis=1) #/des.shape[1] + obj = np.sum(sim*weights) + + #print(xtal); import sys; sys.exit() + return obj + + +def get_input_from_letters(spg, sites, dim=3): + """ + Prepare the input from letters + + Args: + spg (int): 1-230 + sites (list): [['4a'], ['4b']] or [[0], [1]] + dim (int): 0, 1, 2, 3 + + Returns: + space group instance + wp_list + dof + """ + g = Group(spg, dim=dim) + wp_combo = [] + dof = g.get_lattice_dof() + for site in sites: + wp_specie = [] + for s in site: + if type(s) in [int, np.int64]: + wp = g[s] + else: + wp = g.get_wyckoff_position(s) + wp_specie.append(wp) + dof += wp.get_dof() + wp_combo.append(wp_specie) + return g, wp_combo, dof + + +def create_trajectory(dumpfile, trjfile, modes=['Init', 'Iter'], dim=3): + """ + create the ase trajectory from the dump file. + This is used to analyze the performance of structure optimization + """ + keyword_start = 'Space Group' #\n' + keyword_end = 'END'#\n' + + with open(dumpfile, 'r') as f: + lines = f.readlines() + starts, ends = [], [] + for i in range(len(lines)): + l = lines[i].lstrip() + if l.startswith(keyword_start): + starts.append(i) + elif l.startswith(keyword_end): + ends.append(i) + #print(len(starts), len(ends)) + assert(len(starts)==len(ends)) + + atoms = [] + for i, j in zip(starts, ends): + # parse each run + xtal_strs = lines[i:j] + spg = int(xtal_strs[0].split(':')[1]) + l_type = Group(spg).lattice_type + elements = [] + numIons = [] + wps = [] + for count, tmp in enumerate(xtal_strs[1:]): + if tmp.startswith('Element'): + tmp1 = tmp.split(':') + tmp2 = tmp1[1].split() + elements.append(tmp2[0]) + numIons.append(int(tmp2[1])) + wps.append(tmp2[2:]) + else: + line_struc = count + 1 + break + + #print(spg, elements, numIons, wps) + g, wps, dof = get_input_from_letters(spg, wps, dim) + for line_number in range(line_struc, len(xtal_strs)): + tmp0 = xtal_strs[line_number].split(':') + tag = tmp0[0] + if tag in modes: + tmp = tmp0[1].split() + sim = float(tmp[0]) + x = [float(t) for t in tmp[1:]] + + # QZ: to fix + xtal = pyxtal() + xtal = from_1d_rep(x, wps, numIons, l_type, elements, dim) + + struc = xtal.to_ase() + struc.info = {'time': line_number-line_struc, 'fun': sim} + atoms.append(struc) + #print(xtal.lattice) + write(filename=trjfile, images=atoms, format='extxyz') + + + +class mof_builder(object): + """ + Class for generating MOF like structures + + Args: + elements (str): e.g. ['Si', 'O'] + composition (int): e.g. [1, 2] + dim (int): e.g., 0, 1, 2, 3 + db_file (str): default is 'mof.db' + log_file (str): default is 'mof.log' + + Examples + -------- + + To create a new structure instance + + >>> from mof_builder import mof_builder + >>> builder = mof_builder(['P', 'O', 'N'], [1, 1, 1], db_file='PON.db') + """ + + def __init__(self, elements, composition, dim=3, + db_file='mof.db', log_file='mof.log', + verbose=False): + + # Define the chemical system + self.dim = dim + self.elements = elements + self.composition = composition + + # Define the I/O + self.verbose = verbose + logging.getLogger().handlers.clear() + self.log_file = log_file + logging.basicConfig(format = "%(asctime)s| %(message)s", + filename = self.log_file, + level = logging.INFO) + + self.logging = logging + self.db_file = db_file + self.db = database_topology(db_file) + + # Initialize neccessary functions and attributes + self.calculator = None # will be a callable function + self.ref_environments = None # will be a numpy array + self.criteria = {} # will be a dictionary + + def __str__(self): + + s = "\n------MOF Builder------" + s += "\nSystem: " + for element, comp in zip(self.elements, self.composition): + s += "{:s}{:d} ".format(element, comp) + s += "\nDatabase: {}".format(self.db_file) + s += "\nLog_file: {}".format(self.log_file) + if self.calculator is not None: + s += "\nDescriptor: {}".format(self.calculator) + if self.ref_environments is not None: + (d1, d2) = self.ref_environments.shape + s += "Reference enviroments ({:}, {:})".format(d1, d2) + if len(self.criteria.keys()) > 0: + for key in self.criteria.keys(): + s += '\nCriterion_{:}: {:}'.format(key, self.criteria[key]) + return s + + def __repr__(self): + return str(self) + + def set_descriptor_calculator(self, dtype='SO3', mykwargs={}): + """ + Set up the calculator for descriptor computation. + Here we mostly use the pyxtal_ff module + + Arg: + dytpe (str): only SO3 is suppoted now + """ + if dtype == 'SO3': + kwargs = {'lmax': 4, + 'nmax': 2, + 'rcut': 2.2, + 'alpha': 1.5, + #'derivative': False, + 'weight_on': True, + 'neighborlist': 'ase', + } + kwargs.update(mykwargs) + + self.calculator = SO3(**kwargs) + #print(self.calculator)#; import sys; sys.exit() + + def set_reference_enviroments(self, cif_file, substitute=None): + """ + Get the reference enviroments + + Args: + cif_file (str): cif structure + substitute (dict): substitution directory + """ + + if self.calculator is None: + raise RuntimeError("Must call set_descriptor_calculator in advance") + + xtal = pyxtal() + xtal.from_seed(cif_file) + if substitute is not None: xtal.substitute(substitute)#; print(xtal) + xtal.resort_species(self.elements) + + ids = [0] * len(self.elements) + count = 0 + for site in xtal.atom_sites: + for i, element in enumerate(self.elements): + if element == site.specie: + ids[i] = count + break + count += site.multiplicity + if self.verbose: print("ids from Reference xtal", ids) + atoms = xtal.to_ase(resort=False) + self.ref_environments = self.calculator.calculate(atoms, ids)['x'][ids] + if self.verbose: print(self.ref_environments) + self.ref_xtal = xtal + + def set_criteria(self, CN=None, dimension=3, min_density=None, exclude_ii=False): + """ + define the criteria to check if a structure is good + + Args: + CN (int): coordination number + dimension (int): target dimensionality (e.g. we want the 3D structure) + min_density (float): minimum density + exclude_ii (bool): allow the bond between same element + """ + + if CN is not None: + self.criteria["CN"] = CN + if 'cutoff' in CN.keys(): + self.criteria['cutoff'] = CN['cutoff'] + else: + self.criteria['cutoff'] = None + if dimension is not None: + self.criteria["Dimension"] = dimension + if min_density is not None: + self.criteria["MIN_Density"] = min_density + self.criteria['exclude_ii'] = exclude_ii + + def get_input_from_letters(self, spg, wps): + """ + A short cut to get (spg, wps, dof) from the get_input functions + + Args: + spg (int): space group number 1-230 + wps (list): e.g. [['4a', '4b']] + """ + return get_input_from_letters(spg, wps, self.dim) + + def get_input_from_ref_xtal(self, xtal, substitute=None): + """ + Generate the input from a given pyxtal + + Args: + xtal: pyxtal object + substitute: a dictionary to describe chemical substitution + """ + + # Make sure it is a pyxtal object + if type(xtal) == str: + c = pyxtal() + c.from_seed(xtal) + xtal = c + if substitute is not None: xtal.substitute(substitute) + + g = xtal.group + sites = [ [] for _ in self.elements ] + dof = xtal.lattice.dof + for s in xtal.atom_sites: + for i, specie in enumerate(self.elements): + if s.specie == specie: + #wp_combo[i].append(s.wp) + sites[i].append(s.wp.get_label()) + break + dof += s.wp.get_dof() + + return g.number, sites, dof + + def get_similarity(self, xtal): + """ + Compute the similrity for a given xtal + QZ: call calculate_S() + + Args: + xtal: pyxtal object + + Returns: + the similarity value .w.r.t the reference environment + """ + x = xtal.get_1d_rep_x() + return calculate_S(x, xtal, self.ref_environments, + self.calculator) + + def optimize_xtals(self, xtals, ncpu=1, opt_type='local', + T=0.2, niter=20, early_quit=0.02, + add_db=True, symmetrize=False, + minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], + ): + xtals_opt = [] + xs = [] + count = 0 + if ncpu == 1: + for i, xtal in enumerate(xtals): + xtal, sim, _xs = self.optimize_xtal(xtal, i, opt_type, T, + niter, early_quit, + add_db, symmetrize, + minimizers=minimizers, + ) + if xtal is not None: + xtals_opt.append(xtal) + xs.append(_xs) + else: + N_cycle = int(np.ceil(len(xtals)/ncpu)) + print("\n# Parallel Calculation in optimize_xtals", ncpu, N_cycle) + + args_list = [] + for i in range(ncpu): + id1 = i * N_cycle + id2 = min([id1 + N_cycle, len(xtals)]) + wp_libs = [] + #print('cpu', i, id1, id2) + for id in range(id1, id2): + xtal = xtals[id] + x = xtal.get_1d_rep_x() + spg, wps, _ = self.get_input_from_ref_xtal(xtal) + wp_libs.append((x, xtal.group.number, wps)) + + args_list.append((self.dim, + wp_libs, + self.elements, + self.calculator, + self.ref_environments, + opt_type, + T, + niter, + early_quit, + minimizers)) + with ProcessPoolExecutor(max_workers=ncpu) as executor: + results = [executor.submit(minimize_from_x_par, *p) for p in args_list] + for result in results: + _xtals, _xs = result.result() + xtals_opt.extend(_xtals) + xs.extend(_xs) + + # Now process each of the results + valid_xtals = [] + for xtal, _xs in zip(xtals_opt, xs): + status = xtal.check_validity(self.criteria, verbose=self.verbose) + if status: + valid_xtals.append(xtal) + sim1 = self.get_similarity(xtal) + if symmetrize: + pmg = xtal.to_pymatgen() + xtal = pyxtal() + xtal.from_seed(pmg) + if add_db: + self.process_xtal(xtal, [0, sim1], count, xs=_xs) + count += 1 + else: + dicts = {'sim': "{:6.3f}".format(sim1)} + print(xtal.get_xtal_string(dicts)) + + return valid_xtals + + def optimize_xtal(self, xtal, count=0, opt_type='local', + T=0.2, niter=20, early_quit=0.02, + add_db=True, symmetrize=False, + minimizers=[('Nelder-Mead', 100), ('L-BFGS-B', 100)], + filename=None): + """ + Further optimize the input xtal w.r.t reference environment + + Args: + xtal (instance): pyxtal + """ + # Change the angle to a better rep + if xtal.dim==3 and xtal.lattice.ltype in ['triclinic', 'monoclinic']: + xtal.optimize_lattice(standard=True) + x = xtal.get_1d_rep_x() + _, wps, _ = self.get_input_from_ref_xtal(xtal) + + sim0 = self.get_similarity(xtal) + if xtal.lattice is not None: + result = minimize_from_x(x, xtal.dim, xtal.group.number, wps, + self.elements, self.calculator, + self.ref_environments, + opt_type = opt_type, + T = T, + niter = niter, + early_quit = early_quit, + minimizers = minimizers, + filename=filename) + xtal, xs = result + status = xtal.check_validity(self.criteria, verbose=self.verbose) + sim1 = self.get_similarity(xtal) + else: + xtal = None + xs = None + status = False + sim1 = None + xs = None + + if status: + if symmetrize: + pmg = xtal.to_pymatgen() + xtal = pyxtal() + xtal.from_seed(pmg) + if add_db: + self.process_xtal(xtal, [sim0, sim1], count, xs) + else: + dicts = {'sim': "{:6.3f} => {:6.3f}".format(sim0, sim1)} + print(xtal.get_xtal_string(dicts)) + + return xtal, sim1, xs + + def generate_xtal(self, spg, wps, niter, T=0.2, N_max=5, + early_quit=0.03, dump=False, verbose=None, + add_db=True, random_state=None): + """ + Generate a crystal with the desired local environment + + Args: + spg (int): group number + wps (list): list of wps for the disired crystal (e.g., [spg, wp1, wp2]) + n_iter (int): number of iterations for basin hopping + T (float): for basinhopping + N_max (int): number of maximum + early_quit (float): threshhold for early termination + dump (bool): whether or not dump the trajectory + """ + if verbose is None: verbose=self.verbose + xtal, sim = generate_xtal(self.dim, spg, wps, niter, + self.elements, + self.calculator, + self.ref_environments, + self.criteria, + T, N_max, early_quit, dump, + random_state, + verbose=verbose) + + if xtal is not None and xtal.check_validity(self.criteria): + if add_db: + self.process_xtal(xtal, [0, sim], 0) + + return xtal, sim + + def generate_xtals_from_wp_libs(self, wp_libs, N_max=5, ncpu=1, + T=0.2, factor=5, early_quit=0.02): + """ + Run multiple crystal generation from the given wp_libs + This is the core part for structure generation. + + Args: + wp_libs (tuple): (number, spg, wps, dof) + N_max (int): Number of maximum runs + ncpu (int): Num of parallel processes + T (float): basinhopping temperature + factor (int): the number of Basinhopping iterations = factor * dof + early_quit (float): early termination for basinhopping + """ + + # Generate xtals + _args = (self.dim, + self.elements, + self.calculator, + self.ref_environments, + self.criteria, + T, N_max, early_quit) + count = 0 + xtals, sims = [], [] + if ncpu == 1: + for wp_lib in wp_libs: + (number, spg, wps, dof) = wp_lib + xtal, sim = generate_xtal(spg, wps, factor*dof, *_args) + if xtal is not None: + xtals.append(xtal) + sims.append(sim) + + else: + N_cycle = int(np.ceil(len(wp_libs)/ncpu)) + print("\n# Parallel Calculation in generate_xtals_from_wp_libs", ncpu, N_cycle) + + args_list = [] + for i in range(ncpu): + id1 = i * N_cycle + id2 = min([id1 + N_cycle, len(wp_libs)]) + args_list.append((wp_libs[id1:id2], factor) + _args) + + with ProcessPoolExecutor(max_workers=ncpu) as executor: + results = [executor.submit(generate_xtal_par, *p) for p in args_list] + for result in results: + (_xtals, _sims) = result.result() + xtals.extend(_xtals) + sims.extend(_sims) + + for i, xtal in enumerate(xtals): + self.process_xtal(xtal, [0, sims[i]], i) + + return xtals + + def get_wp_libs_from_spglist(self, spg_list, + num_wp=(None, None), + num_fu=(None, None), + num_dof=(1, 10), + per_spg=30, + ncpu=1): + """ + Generate wp choices from the list of space groups + + Args: + spglist (list): list of space group numbers + num_wp (int): a tuple of (min_wp, max_wp) + num_fu (int): a tuple of (min_fu, max_fu) + num_dof (int): a tuple of (min_dof, max_dof) + per_spg (int): maximum number of wp combinations + ncpu (int): number of processors + """ + + print('\nGet wp_libs from the given spglist') + composition = self.composition + (min_wp, max_wp) = num_wp + if min_wp is None: min_wp = len(composition) + if max_wp is None: max_wp = max([min_wp, len(composition)]) + num_wp = (min_wp, max_wp) + + def process_wp_lib(spg, wp_lib): + strs = "{:d} wp combos in space group {:d}".format(len(wp_lib), spg) + print(strs); self.logging.info(strs) + + if len(wp_lib) > per_spg: + ids = np.random.choice(range(len(wp_lib)), per_spg, replace=False) + strs = "Randomly choose {:} wp combinations".format(per_spg) + wp_lib = [wp_lib[x] for x in ids] + print(strs); self.logging.info(strs) + return wp_lib + + # Get wp_libs + wp_libs_total = [] + if ncpu == 1: + for spg in spg_list: + wp_lib = generate_wp_lib([spg], composition, num_wp, num_fu, num_dof) + wp_lib = process_wp_lib(spg, wp_lib) + if len(wp_lib) > 0: + wp_libs_total.extend(wp_lib) + + else: + N_cycle = int(np.ceil(len(spg_list)/ncpu)) + args_list = [] + print("\n# Parallel Calculation in get_wp_libs_from_spglist", ncpu, N_cycle) + for i in range(ncpu): + id1 = i*N_cycle + id2 = min([id1+N_cycle, len(spg_list)]) + args_list.append((spg_list[id1:id2], + composition, + num_wp, + num_fu, + num_dof)) + + #collect the results + with ProcessPoolExecutor(max_workers=ncpu) as executor: + results = [executor.submit(generate_wp_lib_par, *p) for p in args_list] + for result in results: + (spgs, wp_libs) = result.result() + for spg, wp_lib in zip(spgs, wp_libs[0]): + wp_lib = process_wp_lib(spg, wp_lib) + wp_libs_total.extend(wp_lib) + + return sorted(wp_libs_total) + + def get_wp_libs_from_xtals(self, db_file=None, + num_wp=(None, None), + num_dof=(1, 20), + num_atoms=[1, 500]): + """ + For each struc in the database, get the (spg, wp) info + + Args: + df_file (str): database file + num_wp (int): tuple of (min_wp, max_wp) + num_dof (int): tuple of (min_dof, max_dof) + num_atoms (int): tuple of (min_at, max_at) + """ + + (min_dof, max_dof) = num_dof + (min_wp, max_wp) = num_wp + (min_at, max_at) = num_atoms + if min_wp is None: min_wp = len(self.composition) + if max_wp is None: max_wp = max([min_wp, len(self.composition)]) + + wp_libs_total = [] + if db_file is not None: + strs = "====Loading the structures from {:}".format(db_file) + print("\n", strs) + self.logging.info(strs) + + with connect(db_file) as db: + for row in db.select(): + atoms = db.get_atoms(row.id) + if min_at <= len(atoms) <= max_at: + xtal = pyxtal() + xtal.from_seed(atoms, tol=0.1) + if min_wp <= len(xtal.atom_sites) <= max_wp and \ + min_dof <= xtal.get_dof() <= max_dof and \ + xtal.check_validity(self.criteria): + t0 = time() + sim0 = self.get_similarity(xtal) + dicts = {'sim': "0 => {:6.3f}".format(sim0)} + strs = xtal.get_xtal_string(dicts) + print(strs, "{:.6f}".format(time()-t0)) + spg, wps, dof = self.get_input_from_ref_xtal(xtal) + wp_libs_total.append((sum(xtal.numIons), spg, wps, dof)) + + return sorted(wp_libs_total) + + + def import_structures(self, db_file, ids=(None, None), + check=True, same_group=True, + spglist=range(1, 231), bounds=[1, 500], + relax=True, + ): + """ + Import the structures from the external ase database + + Args: + db_file (str): ase database + check (boolean): whether or not + spglist (int): list of spg numbers + bounds: number of atoms + energy (boolean): add energy or not + """ + [lb, ub] = bounds + count = 0 + strs = "====Loading the structures from {:}".format(db_file) + print("\n", strs) + self.logging.info(strs) + with connect(db_file) as db: + (min_id, max_id) = ids + if min_id is None: min_id = 1 + if max_id is None: max_id = db.count() + 100000 + for row in db.select(): + if min_id <= row.id <= max_id: + atoms = row.toatoms() + xtal = pyxtal() + try: + xtal.from_seed(atoms, tol=0.1) + if lb <= len(atoms) <= ub and xtal.group.number in spglist: + status = True + else: + status = False + except: + print("Error in loading structure") + status = False + + # Check structures + if status: + # Relax the structure + if relax: + xtal, sim, _ = self.optimize_xtal(xtal, add_db=False) + + if xtal is not None: + energy = row.ff_energy if hasattr(row, 'ff_energy') else None + topology = row.topology if hasattr(row, 'topology') else None + self.process_xtal(xtal, [0, sim], count, + energy=energy, + topology=topology, + same_group=same_group) + count += 1 + + def process_xtal(self, xtal, sim, count=0, xs=None, energy=None, + topology=None, same_group=True, db=None, check=False): + """ + Check, print and add xtal to the database + + Args: + xtal: pyxtal object + sim: list of two similarity numbers + count (int): id of the structure in the database + xs (tuple): list of reps before and after optimization + energy (float): the energy value + same_group (bool): keep the same group or not + db (str): db path + check (bool): whether or not check if the structure is a duplicate + """ + if db is None: db = self.db + if check: + status = db.check_new_structure(xtal, same_group) + else: + status = True + header = "{:4d}".format(count) + dicts = {'energy': energy, + 'status': status, + 'sim': "{:12.3f} => {:6.3f}".format(sim[0], sim[1]) + } + strs = xtal.get_xtal_string(dicts, header) + print(strs) + self.logging.info(strs) + kvp = { + 'similarity0': sim[0], + 'similarity': sim[1], + } + if xs is not None: + kvp['x_init'] = np.array2string(xs[0]) + kvp['x_opt'] = np.array2string(xs[1]) + if energy is not None: kvp['ff_energy'] = energy + if topology is not None: kvp['topology'] = topology + if status: db.add_xtal(xtal, kvp) + + +""" +Custom step-function for basin hopping optimization +""" +class RandomDisplacementBounds(object): + """ + random displacement with bounds: + see: https://stackoverflow.com/a/21967888/2320035 + """ + def __init__(self, xmin, xmax, id1, id2, stepsize=0.1, dumpfile=None): + self.xmin = xmin # bound_min + self.xmax = xmax # bound_max + self.id1 = id1 + self.id2 = id2 + self.stepsize = stepsize + self.dumpfile = dumpfile + + def __call__(self, x): + """ + Move the step proportionally within the bound + """ + # To strongly rebounce the values hitting the wall + for i in range(self.id2): + if abs(x[i] - self.xmax[i]) < 0.1: + #print("To strongly rebounce max values", x[i], self.xmax[i]) + x[i] *= 0.5 + elif abs(x[i] - self.xmin[i]) < 0.1: + #print("To strongly rebounce min values", x[i], self.xmin[i]) + x[i] *= 2.0 + + random_step = np.random.uniform(low=-self.stepsize, + high=self.stepsize, + size=x.shape) + xnew = x + random_step + # Cell + coefs = 1+0.2*(np.random.sample(self.id2)-0.5) + xnew[:self.id2] *= coefs + + # xyz + xnew[self.id1:] -= np.floor(xnew[self.id1:]) + #xnew = np.maximum(self.xmin, xnew) + #xnew = np.minimum(self.xmax, xnew) + + # Randomly introduce compression to prevent non-pd xtal + if np.random.random() < 0.5: + id = np.argmax(xnew[:self.id2]) + xnew[id] *= 0.8 + xnew = np.maximum(self.xmin, xnew) + xnew = np.minimum(self.xmax, xnew) + + #min_step = np.maximum(self.xmin - x, -self.stepsize) + #max_step = np.minimum(self.xmax - x, self.stepsize) + #random_step = np.random.uniform(low=min_step, high=max_step, size=x.shape) + #xnew = x + random_step + if self.dumpfile is not None: + with open(self.dumpfile, 'a+') as f0: + # Initial value + strs = 'Init: {:9.3f} '.format(10.0) + for x0 in xnew: strs += '{:8.4f} '.format(x0) + strs += '\n' + f0.write(strs) + + return xnew + +if __name__ == "__main__": + + xtal = pyxtal() + xtal.from_spg_wps_rep(194, ['2c', '2b'], [2.46, 6.70]) + cif_file = xtal.to_pymatgen() + builder = mof_builder(['C'], [1], db_file='reaxff.db', verbose=False)#True) + builder.set_descriptor_calculator(mykwargs={'rcut': 1.9}) + builder.set_reference_enviroments(cif_file) + builder.set_criteria(CN={'C': [3]}) + print(builder) + print(builder.ref_xtal) + if False: + wp_libs = builder.get_wp_libs_from_spglist([191, 179], ncpu=1) + for wp_lib in wp_libs[:4]: print(wp_lib) + builder.generate_xtals_from_wp_libs(wp_libs[4:8], ncpu=2, N_max=4, early_quit=0.05) + + if False: + spg, wps = 179, ['6a', '6a', '6a', '6a'] + xtals = [] + for x in [ + #[ 9.6244, 2.5459, 0.1749, 0.7701, 0.4501, 0.6114], + #[15.0223, 1.5013, 0.8951, 0.6298, 0.4530, 0.1876], + #[10.0129, 2.6424, 0.3331, 0.7246, 0.4719, 0.8628], + #[10.2520, 3.1457, 0.2367, 0.6994, 0.2522, 0.6533], + #[ 9.3994, 2.5525, 0.3072, 0.8414, 0.3480, 0.9638], + [11.1120, 2.6428, 0.2973, 0.7513, 0.4236, 0.8777], + [7.9522, 2.6057, 0.5922, 0.9268, 0.6081, 0.3077], + ]: + xtal = pyxtal() + xtal.from_spg_wps_rep(spg, wps, x, ['C']*len(wps)) + xtals.append(xtal) + builder.optimize_xtal(xtal) + builder.optimize_xtals(xtals) diff --git a/requirements.txt b/requirements.txt index ea6fa64c..5d1d314e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,6 +18,10 @@ importlib-metadata==8.0.0 # via pyxtal (setup.py) joblib==1.4.2 # via pymatgen +juliacall==0.9.22 + # via pyxtal (setup.py) +juliapkg==0.1.13 + # via juliacall kiwisolver==1.4.5 # via matplotlib latexcodec==3.0.0 @@ -83,6 +87,8 @@ scipy==1.14.0 # pyxtal (setup.py) # ase # pymatgen +semver==3.0.2 + # via juliapkg six==1.16.0 # via # pybtex diff --git a/setup.py b/setup.py index 10287450..f1b4f9e0 100755 --- a/setup.py +++ b/setup.py @@ -11,7 +11,7 @@ setup( name="pyxtal", - version="1.0.1", + version="1.0.2", author="Scott Fredericks, Kevin Parrish, Qiang Zhu", author_email="alecfans@gmail.com", description="Python code for generation of crystal structures based on symmetry constraints.", @@ -22,6 +22,7 @@ "pyxtal", "pyxtal.database", "pyxtal.interface", + "pyxtal.lego", "pyxtal.optimize", "pyxtal.potentials", "pyxtal.database.cifs", @@ -41,6 +42,7 @@ "Operating System :: OS Independent", ], install_requires=[ + "juliacall", "spglib>=1.10.4", "pymatgen>=2024.3.1", "pandas>=2.0.2", @@ -48,6 +50,7 @@ "ase>=3.23.0", "scipy>=1.7.3", "numpy>=1.26,<2", # prevent the use of numpy2 + "vasprun-xml>=1.0.4", # prevent the use of numpy2 "importlib_metadata>=1.4", "typing-extensions>=4.12", ],