diff --git a/apex/core/common_prop.py b/apex/core/common_prop.py index dc179dd..d54f981 100644 --- a/apex/core/common_prop.py +++ b/apex/core/common_prop.py @@ -11,6 +11,7 @@ from apex.core.property.Surface import Surface from apex.core.property.Vacancy import Vacancy from apex.core.property.Phonon import Phonon +from apex.core.property.DecohesionEnergy import DecohesionEnergy from apex.core.lib.utils import create_path from apex.core.lib.util import collect_task from apex.core.lib.dispatcher import make_submission @@ -39,6 +40,8 @@ def make_property_instance(parameters, inter_param): return Gamma(parameters, inter_param) elif prop_type == "phonon": return Phonon(parameters, inter_param) + elif prop_type == "DecohesionEnergy": + return DecohesionEnergy(parameters, inter_param) else: raise RuntimeError(f"unknown APEX type {prop_type}") diff --git a/apex/core/property/DecohesionEnergy.py b/apex/core/property/DecohesionEnergy.py new file mode 100644 index 0000000..7950763 --- /dev/null +++ b/apex/core/property/DecohesionEnergy.py @@ -0,0 +1,300 @@ +import glob +import json +import logging +import os +import re +import dpdata +import numpy as np +from monty.serialization import dumpfn, loadfn +from pymatgen.core.structure import Structure +from pymatgen.core.surface import SlabGenerator + +from apex.core.calculator.lib import abacus_utils +from apex.core.calculator.lib import vasp_utils +from apex.core.property.Property import Property +from apex.core.refine import make_refine +from apex.core.reproduce import make_repro, post_repro +from dflow.python import upload_packages +upload_packages.append(__file__) + +class DecohesionEnergy(Property): + def __init__(self, parameter, inter_param=None): + ''' + core parameter for make_confs[POSCAR]： + min_slab_size max_vacuum_size vacuum_size_step + miller_index + ''' + parameter["reproduce"] = parameter.get("reproduce", False) + self.reprod = parameter["reproduce"] + if not self.reprod: + if not ("init_from_suffix" in parameter and "output_suffix" in parameter): + self.min_slab_size = parameter["min_slab_size"] + parameter["pert_xz"] = parameter.get("pert_xz", 0.01) + self.pert_xz = parameter["pert_xz"] + parameter["max_vacuum_size"] = parameter.get("max_vacuum_size", 15) + self.max_vacuum_size = parameter["max_vacuum_size"] + parameter["vacuum_size_step"]=parameter.get("vacuum_size_step", 1) + self.vacuum_size_step = parameter["vacuum_size_step"] + self.miller_index = tuple(parameter["miller_index"]) + parameter["cal_type"] = parameter.get("cal_type", "relaxation") + default_cal_setting = { + "relax_pos": False, + "relax_shape": False, + "relax_vol": False, + } + else: + parameter["cal_type"] = "static" + self.cal_type = parameter["cal_type"] + default_cal_setting = { + "relax_pos": False, + "relax_shape": False, + "relax_vol": False, + } + parameter["init_from_suffix"] = parameter.get("init_from_suffix", "00") + self.init_from_suffix = parameter["init_from_suffix"] + self.cal_type = parameter["cal_type"] + parameter["cal_setting"] = parameter.get("cal_setting", default_cal_setting) + for key in default_cal_setting: + parameter["cal_setting"].setdefault(key, default_cal_setting[key]) + self.cal_setting = parameter["cal_setting"] + self.parameter = parameter + self.inter_param = inter_param if inter_param != None else {"type": "vasp"} + + def make_confs(self, path_to_work, path_to_equi, refine=False): + path_to_work = os.path.abspath(path_to_work) + if os.path.exists(path_to_work): + logging.warning("%s already exists" % path_to_work) + else: + os.makedirs(path_to_work) + path_to_equi = os.path.abspath(path_to_equi) + + task_list = [] + cwd = os.getcwd() + + if self.reprod: + print("surface reproduce starts") + if "init_data_path" not in self.parameter: + raise RuntimeError("please provide the initial data path to reproduce") + init_data_path = os.path.abspath(self.parameter["init_data_path"]) + task_list = make_repro( + self.inter_param, + init_data_path, + self.init_from_suffix, + path_to_work, + self.parameter.get("reprod_last_frame", True), + ) + + else: + if refine: + logging.info("surface refine starts") + task_list = make_refine( + self.parameter["init_from_suffix"], + self.parameter["output_suffix"], + path_to_work, + ) + # record miller + init_from_path = re.sub( + self.parameter["output_suffix"][::-1], + self.parameter["init_from_suffix"][::-1], + path_to_work[::-1], + count=1, + )[::-1] + task_list_basename = list(map(os.path.basename, task_list)) + + for ii in task_list_basename: + init_from_task = os.path.join(init_from_path, ii) + output_task = os.path.join(path_to_work, ii) + os.chdir(output_task) + if os.path.isfile("decohesion_energy.json"): + os.remove("decohesion_energy.json") + if os.path.islink("decohesion_energy.json"): + os.remove("decohesion_energy.json") + os.symlink( + os.path.relpath(os.path.join(init_from_task, "decohesion_energy.json")), + "decohesion_energy.json",) + else: + if self.inter_param["type"] == "abacus": + CONTCAR = abacus_utils.final_stru(path_to_equi) + POSCAR = "STRU" + else: + # refine = false && reproduce = false && self.inter_param["type"]== "vasp" + CONTCAR = "CONTCAR" + POSCAR = "POSCAR" + equi_contcar = os.path.join(path_to_equi, CONTCAR) + if not os.path.exists(equi_contcar): + raise RuntimeError("please do relaxation first") + + if self.inter_param["type"] == "abacus": + stru = dpdata.System(equi_contcar, fmt="stru") + stru.to("contcar", "CONTCAR.tmp") + ptypes = vasp_utils.get_poscar_types("CONTCAR.tmp") + ss = Structure.from_file("CONTCAR.tmp") + os.remove("CONTCAR.tmp") + else: + ptypes = vasp_utils.get_poscar_types(equi_contcar) + # element type 读取 vasp 第五行 + ss = Structure.from_file(equi_contcar) + + # gen POSCAR of Slab + all_slabs = [] + vacuum = [] + num = 0 + while self.vacuum_size_step * num <= self.max_vacuum_size: + vacuum_size = self.vacuum_size_step * num + slabs = self.__gen_slab_pmg(ss, self.miller_index, self.min_slab_size, vacuum_size) + num = num + 1 + all_slabs.append(slabs) + vacuum.append(vacuum_size) + + os.chdir(path_to_work) + if os.path.exists(POSCAR): + os.remove( + POSCAR) + os.symlink(os.path.relpath(equi_contcar), POSCAR) + for ii in range(len(all_slabs)): + output_task = os.path.join(path_to_work, "task.%06d" % ii) + os.makedirs(output_task, exist_ok=True) + os.chdir(output_task) + for jj in [ + "INCAR", + "POTCAR", + "POSCAR", + "conf.lmp", + "in.lammps", + "STRU", + ]: + if os.path.exists(jj): + os.remove(jj) + task_list.append(output_task) + + logging.info( + "# %03d generate " % ii, + output_task, + " \t %d atoms" % len(all_slabs[ii].sites), + ) + + all_slabs[ii].to("POSCAR.tmp", "POSCAR") + vasp_utils.regulate_poscar("POSCAR.tmp", "POSCAR") + vasp_utils.sort_poscar("POSCAR", "POSCAR", ptypes) + vasp_utils.perturb_xz("POSCAR", "POSCAR", self.pert_xz) + if self.inter_param["type"] == "abacus": + abacus_utils.poscar2stru("POSCAR", self.inter_param, "STRU") + os.remove("POSCAR") + decohesion_energy = {"miller_index": self.miller_index, "vacuum_size":vacuum[ii]} + dumpfn(decohesion_energy, "decohesion_energy.json", indent=4) + os.chdir(cwd) + return task_list + + def post_process(self, task_list): + pass + + def task_type(self): + return self.parameter["type"] + + def task_param(self): + return self.parameter + + def _compute_lower(self, output_file, all_tasks, all_res) -> [dict, str]: + output_file = os.path.abspath(output_file) + res_data = {} + ptr_data = os.path.dirname(output_file) + "\n" + if not self.reprod: + ''' + equi_path = os.path.abspath( + os.path.join( + os.path.dirname(output_file), "../relaxation/relax_task" + ) + ) + equi_result = loadfn(os.path.join(equi_path, "result.json")) + equi_epa = equi_result["energies"][-1] / np.sum( + equi_result["atom_numbs"] + ) + ''' + vacuum_size_step = loadfn(os.path.join(os.path.dirname(output_file), "param.json"))["vacuum_size_step"] + ptr_data += ("Miller Index: " + str(loadfn(os.path.join(os.path.dirname(output_file), "param.json"))["miller_index"]) + "\n") + ptr_data += "Vacuum_size(e-10 m):\tDecohesion_E(J/m^2) Decohesion_S(Pa)\n" + pre_task_result = loadfn(os.path.join(all_tasks[0],"result_task.json")) + pre_evac = 0 + equi_evac = pre_task_result["energies"][-1] + for ii in all_tasks: + task_result = loadfn(os.path.join(ii, "result_task.json")) + AA = np.linalg.norm( + np.cross(task_result["cells"][0][0], task_result["cells"][0][1]) + ) + structure_dir = os.path.basename(ii) + Cf = 1.60217657e-16 / 1e-20 * 0.001 + evac = (task_result["energies"][-1] - equi_evac) / AA * Cf + vacuum_size = loadfn(os.path.join(ii, "decohesion_energy.json"))["vacuum_size"] + stress = (evac - pre_evac) / vacuum_size_step * 1e10 + + ptr_data += "%-30s % 7.3f %10.3e \n" % ( + str(vacuum_size) + "-" + structure_dir + ":", + evac, + stress, + ) + res_data[str(vacuum_size) + "_" + structure_dir] = [ + evac, + stress, + vacuum_size, + ] + pre_evac = evac + + else: + if "init_data_path" not in self.parameter: + raise RuntimeError("please provide the initial data path to reproduce") + init_data_path = os.path.abspath(self.parameter["init_data_path"]) + res_data, ptr_data = post_repro( + init_data_path, + self.parameter["init_from_suffix"], + all_tasks, + ptr_data, + self.parameter.get("reprod_last_frame", True), + ) + + with open(output_file, "w") as fp: + json.dump(res_data, fp, indent=4) + + return res_data, ptr_data + + def __gen_slab_pmg(self, structure: Structure, + plane_miller, slab_size, vacuum_size) -> Structure: + + # Generate slab via Pymatgen + slabGen = SlabGenerator(structure, miller_index=plane_miller, + min_slab_size=slab_size, min_vacuum_size=0, + center_slab=True, in_unit_planes=False, + lll_reduce=True, reorient_lattice=False, + primitive=False) + slabs_pmg = slabGen.get_slabs(ftol=0.001) + slab = [s for s in slabs_pmg if s.miller_index == plane_miller][0] + # If a transform matrix is passed, reorient the slab + order = zip(slab.frac_coords, slab.species) + c_order = sorted(order, key=lambda x: x[0][2]) + sorted_frac_coords = [] + sorted_species = [] + for (frac_coord, species) in c_order: + sorted_frac_coords.append(frac_coord) + sorted_species.append(species) + # add vacuum layer to the slab with height unit of angstrom + a, b, c = slab.lattice.matrix + slab_height = slab.lattice.matrix[2][2] + if slab_height >= 0: + self.is_flip = False + elong_scale = 1 + (vacuum_size / slab_height) + else: + self.is_flip = True + elong_scale = 1 + (-vacuum_size / slab_height) + new_lattice = [a, b, elong_scale * c] + new_frac_coords = [] + for ii in range(len(sorted_frac_coords)): + coord = sorted_frac_coords[ii].copy() + coord[2] = coord[2] / elong_scale + new_frac_coords.append(coord) + # 建立 slab 结构 + + slab_new = Structure(lattice=np.matrix(new_lattice), + coords=new_frac_coords, species=sorted_species) + + return slab_new + + diff --git a/apex/reporter/DashReportApp.py b/apex/reporter/DashReportApp.py index 76fa746..d77e593 100644 --- a/apex/reporter/DashReportApp.py +++ b/apex/reporter/DashReportApp.py @@ -33,6 +33,8 @@ def return_prop_class(prop_type: str): return GammaReport elif prop_type == 'phonon': return PhononReport + elif prop_type == 'DecohesionEnergy': + return DecohesionEnergyReport def return_prop_type(prop: str): @@ -140,6 +142,9 @@ def generate_layout(self): if default_dataset: break + radio_inline = False + if len(self.all_dimensions) > 10: + radio_inline = True layout = html.Div( [ html.H1("APEX Results Visualization Report", style={'textAlign': 'center'}), @@ -147,7 +152,7 @@ def generate_layout(self): dcc.RadioItems( id='confs-radio', options=[{'label': name, 'value': name} for name in self.all_dimensions], - value=default_dimension, + value=default_dimension, inline=radio_inline, style={"fontSize": UI_FRONTSIZE} ), html.Br(), diff --git a/apex/reporter/property_report.py b/apex/reporter/property_report.py index a64a066..9c4d228 100644 --- a/apex/reporter/property_report.py +++ b/apex/reporter/property_report.py @@ -117,6 +117,87 @@ def dash_table(res_data: dict, decimal: int = 3, **kwargs) -> dash_table.DataTab return table, df +class DecohesionEnergyReport(PropertyReport): + @staticmethod + def plotly_graph(res_data: dict, name: str, **kwargs): + decohesion_e = [values[0] for values in res_data.values()] + stress = [values[1] for values in res_data.values()] + vacuum_size = [values[2] for values in res_data.values()] + vacuum_size = [str(item) for item in vacuum_size] + df = pd.DataFrame({ + "separation distance (A)": vacuum_size, + "Decohesion energy (J/m^2)": decohesion_e, + "Decohesion stress (GPa)": [s / 1e9 for s in stress], + }) + trace_E = go.Scatter( + name=f"{name} Decohesion Energy", + x=df['separation distance (A)'], + y=df['Decohesion energy (J/m^2)'], + mode='lines+markers', + yaxis='y1' + ) + + trace_S = go.Scatter( + name=f"{name} Decohesion Stress", + x=df['separation distance (A)'], + y=df['Decohesion stress (GPa)'], + mode='lines+markers', + yaxis='y2' + ) + layout = go.Layout( + title=dict( + text='Decohesion Energy and Stress', + x=0.5, # 标题居中 + xanchor='center' + ), + xaxis=dict( + title_text="separation distance (A)", + title_font=dict( + size=18, + color="#7f7f7f" + ) + ), + yaxis=dict( + title="Decohesion energy (J/m^2)", + titlefont=dict( + size=18, + color="#7f7f7f" + ) + ), + yaxis2=dict( + title="Decohesion stress (GPa)", + titlefont=dict( + size=18, + color="#7f7f7f" + ), + overlaying='y', + side='right' + ) + ) + trace = [trace_E, trace_S] + return trace, layout + + @staticmethod + def dash_table(res_data: dict, decimal: int = 3, **kwargs) -> dash_table.DataTable: + decohesion_e = [values[0] for values in res_data.values()] + stress = [values[1] for values in res_data.values()] + vacuum_size = [values[2] for values in res_data.values()] + vacuum_size = [str(item) for item in vacuum_size] + df = pd.DataFrame({ + "separation distance (A)": vacuum_size, + "Decohesion energy (J/m^2)": round_format(decohesion_e, decimal), + "Decohesion stress (GPa)": round_format([s / 1e9 for s in stress], decimal), + }) + table = dash_table.DataTable( + data=df.to_dict('records'), + columns=[{'name': i, 'id': i} for i in df.columns], + style_table={'width': TABLE_WIDTH, + 'minWidth': TABLE_MIN_WIDTH, + 'overflowX': 'auto'}, + style_cell={'textAlign': 'left'} + ) + return table, df + class ElasticReport(PropertyReport): @staticmethod diff --git a/setup.py b/setup.py index 9cb7cad..0020dd5 100644 --- a/setup.py +++ b/setup.py @@ -14,6 +14,7 @@ url="https://github.com/deepmodeling/APEX.git", packages=setuptools.find_packages(), install_requires=[ + "numpy==1.26.4", "pydflow>=1.7.83", "pymatgen>=2023.8.10", 'pymatgen-analysis-defects>=2023.8.22',