Merge pull request #89 from deepmodeling/devel-1.2.0

We added the option to specify the modulus approximation method (“voigt”, “reuss” or “vrh”) in the elasticity calculation workflow in 1.2.12.

Author: kevinwenminion

.github/workflows/main.yml

@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.10"]
 
     steps:
     - uses: actions/checkout@master

README.md

@@ -413,6 +413,7 @@ Below are three examples (for detailed explanations of each parameter, please re
   | shear_deform | Float          | 0.01    | The deformation in other directions, default = 1e-2                                                                               |
   | conventional | Bool           | False   | Whether adopt conventional cell for deformation                                                                                   |
   | ieee         | Bool           | False   | Whether rotate relaxed structure into IEEE-standard format before deformation ([ref](https://ieeexplore.ieee.org/document/26560)) |
+  | modulus_type | String         | "voigt"   | Bulk and shear modulus average type (default: "voigt"). Choose from "voigt", "reuss" and "vrh"                                    |
 
 ##### 3.1.2.3. Surface
   | Key words | Data structure | Example | Description |

apex/__init__.py

@@ -1,5 +1,5 @@
 import os
-__version__ = '1.2.10'
+__version__ = '1.2.12'
 LOCAL_PATH = os.getcwd()
 
 

apex/core/calculator/Lammps.py

@@ -14,7 +14,8 @@
     inter_snap,
     inter_gap,
     inter_rann,
-    inter_mace
+    inter_mace,
+    inter_nep
 )
 from .Task import Task
 from dflow.python import upload_packages
@@ -55,6 +56,8 @@ def set_inter_type_func(self):
             self.inter_func = inter_rann
         elif self.inter_type == "mace":
             self.inter_func = inter_mace
+        elif self.inter_type == "nep":
+            self.inter_func = inter_nep
         else:
             self.inter_func = inter_eam_alloy
 

apex/core/calculator/__init__.py

@@ -8,4 +8,5 @@
     'gap',
     'rann',
     'mace',
+    'nep'
 ]

apex/core/calculator/lib/abacus_utils.py

@@ -4,7 +4,7 @@
 
 import dpdata
 import numpy as np
-from pymatgen.core.structure import Structure
+from pymatgen.core import Structure
 
 from . import abacus_scf
 from dflow.python import upload_packages

apex/core/calculator/lib/lammps_utils.py

@@ -138,6 +138,16 @@ def inter_mace(param):
     ret += line
     return ret
 
+def inter_nep(param):
+    ret = ""
+    line = "pair_style      nep \n"
+    line += "pair_coeff      * * %s " % param["model_name"][0]
+    for ii in param["param_type"]:
+        line += ii + " "
+    line += "\n"
+    ret += line
+    return ret
+
 
 def inter_snap(param):
     ret = ""

apex/core/common_equi.py

@@ -5,6 +5,7 @@
 import dpdata
 from monty.serialization import dumpfn
 from pymatgen.core.structure import Structure
+from pymatgen.analysis.structure_matcher import StructureMatcher
 from apex.core.calculator.lib import abacus_utils
 from apex.core.lib import crys
 from apex.core.calculator.calculator import make_calculator
@@ -14,7 +15,7 @@
 from apex.core.structure import StructureInfo
 from dflow.python import upload_packages
 upload_packages.append(__file__)
-lammps_task_type = ['deepmd', 'eam_alloy', 'meam', 'eam_fs', 'meam_spline', 'snap', 'gap', 'rann', 'mace']
+lammps_task_type = ['deepmd', 'eam_alloy', 'meam', 'eam_fs', 'meam_spline', 'snap', 'gap', 'rann', 'mace', 'nep']
 
 
 def make_equi(confs, inter_param, relax_param):
@@ -201,13 +202,23 @@ def post_equi(confs, inter_param):
         except FileNotFoundError:
             logging.warning(f"No CONTCAR found in {ii}, skip")
             continue
+        try:
+            init_ss = Structure.from_file(poscar)
+        except FileNotFoundError:
+            logging.warning(f"No POSCAR found in {ii}, skip")
+            continue
         st = StructureInfo(ss)
+        matcher = StructureMatcher()
+        is_match = matcher.fit(init_ss, ss)
+        if not is_match:
+            logging.warning(f"Structure mismatch after relaxation in {ii}")
         struct_info_dict = {
             "space_group_symbol": st.space_group_symbol,
             "space_group_number": st.space_group_number,
             "point_group_symbol": st.point_group_symbol,
             "crystal_system": st.crystal_system,
             "lattice_type": st.lattice_type,
+            "mismatch": not is_match,
         }
 
         dumpfn(struct_info_dict, os.path.join(ii, "structure.json"), indent=4)

apex/core/common_prop.py

@@ -1,7 +1,7 @@
 import glob
 import os
 from multiprocessing import Pool
-from monty.serialization import dumpfn
+from monty.serialization import dumpfn, loadfn
 
 from apex.core.calculator.calculator import make_calculator
 from apex.core.property.Elastic import Elastic
@@ -54,6 +54,12 @@ def make_property(confs, inter_param, property_list):
     conf_dirs.sort()
     for ii in conf_dirs:
         sepline(ch=ii, screen=True)
+        path_to_equi = os.path.join(ii, "relaxation", "relax_task")
+        try:
+            structure_dict = loadfn(os.path.join(path_to_equi, "structure.json"))
+        except FileNotFoundError:
+            structure_dict = {}
+        mismatch = structure_dict.get("mismatch", False)
         for jj in property_list:
             do_refine, suffix = handle_prop_suffix(jj)
             if not suffix:
@@ -63,8 +69,11 @@ def make_property(confs, inter_param, property_list):
             # determine the suffix: from scratch or refine
 
             property_type = jj["type"]
-            path_to_equi = os.path.join(ii, "relaxation", "relax_task")
             path_to_work = os.path.join(ii, property_type + "_" + suffix)
+            skip_mismatch = jj.get("skip_mismatch", False)
+            if mismatch and skip_mismatch:
+                print("Skip mismatched structure")
+                continue
 
             create_path(path_to_work)
 

apex/core/property/Elastic.py

@@ -35,6 +35,8 @@ def __init__(self, parameter, inter_param=None):
             self.conventional = parameter["conventional"]
             parameter.setdefault("ieee", False)
             self.ieee = parameter["ieee"]
+            parameter.setdefault("modulus_type", "voigt")
+            self.modulus_type = parameter["modulus_type"]
         parameter.setdefault("cal_type", "relaxation")
         self.cal_type = parameter["cal_type"]
         default_cal_setting = {
@@ -272,19 +274,27 @@ def _compute_lower(self, output_file, all_tasks, all_res):
             res_data["elastic_tensor"].append(c_ii)
             ptr_data += "\n"
 
-        BV = et.k_voigt / 1e4
-        GV = et.g_voigt / 1e4
+        if self.modulus_type == "voigt":
+            BV = et.k_voigt / 1e4
+            GV = et.g_voigt / 1e4
+        elif self.modulus_type == "reuss":
+            BV = et.k_reuss / 1e4
+            GV = et.g_reuss / 1e4
+        elif self.modulus_type == "vrh":
+            BV = et.k_vrh / 1e4
+            GV = et.g_vrh / 1e4
+
         EV = 9 * BV * GV / (3 * BV + GV)
         uV = 0.5 * (3 * BV - 2 * GV) / (3 * BV + GV)
 
-        res_data["BV"] = BV
-        res_data["GV"] = GV
-        res_data["EV"] = EV
-        res_data["uV"] = uV
-        ptr_data += "# Bulk   Modulus BV = %.2f GPa\n" % BV
-        ptr_data += "# Shear  Modulus GV = %.2f GPa\n" % GV
-        ptr_data += "# Youngs Modulus EV = %.2f GPa\n" % EV
-        ptr_data += "# Poission Ratio uV = %.2f\n " % uV
+        res_data["B"] = BV
+        res_data["G"] = GV
+        res_data["E"] = EV
+        res_data["u"] = uV
+        ptr_data += "# Bulk   Modulus B = %.2f GPa\n" % BV
+        ptr_data += "# Shear  Modulus G = %.2f GPa\n" % GV
+        ptr_data += "# Youngs Modulus E = %.2f GPa\n" % EV
+        ptr_data += "# Poission Ratio u = %.2f\n " % uV
 
         dumpfn(res_data, output_file, indent=4)