Merge pull request #47 from ZLI-afk/v1.0

v1.1.1

README.md

@@ -3,6 +3,7 @@
 [APEX](https://github.com/deepmodeling/APEX): Alloy Property EXplorer using simulations, is a component of the [AI Square](https://aissquare.com/) project that involves the restructuring of the [DP-Gen](https://github.com/deepmodeling/dpgen) `auto_test` module to develop a versatile and extensible Python package for general alloy property testing. This package enables users to conveniently establish a wide range of property-test workflows by utilizing various computational approaches, including support for LAMMPS, VASP, and ABACUS.
 
 ## New Features Update (v1.0)
+* Add calculation function of `phonon` spectrum (v1.1.0)
 * Decouple property calculations into individual sub-workflow to facilitate the customization of complex property functions
 * Support one-click parallel submission of multiple workflows
 * Support `run` step in the single step test mode (Interaction method similar to `auto_test`)
@@ -28,6 +29,7 @@
         - [3.1.2.4. Vacancy](#3124-vacancy)
         - [3.1.2.5. Interstitial](#3125-interstitial)
         - [3.1.2.6. Gamma Line](#3126-gamma-line)
+        - [3.1.2.7. Phonon Spectrum](#3127-phonon-spectrum)
     - [3.2. Command](#32-command)
       - [3.2.1. Workflow Submission](#321-workflow-submission)
       - [3.2.2. Single-Step Test](#322-single-step-test)
@@ -43,7 +45,7 @@ APEX adopts the functionality of the second-generation alloy properties calculat
 The comprehensive architecture of APEX is demonstrated below:
 
 <div>
-    <img src="./docs/images/apex_demo_v1.0.png" alt="Fig1" style="zoom: 35%;">
+    <img src="./docs/images/apex_demo.png" alt="Fig1" style="zoom: 35%;">
     <p style='font-size:1.0rem; font-weight:none'>Figure 1. APEX schematic diagram</p>
 </div>
 
@@ -344,7 +346,31 @@ The parameters related to Gamma line calculation are listed below:
       "n_steps":         10
 	}
   ```
-  **It should be noted that for various crystal structures, users can further define slip parameters within the respective nested dictionaries, which will be prioritized for adoption. In above example, the slip system configuration within the "hcp" dictionary will be utilized.**
+  It should be noted that for various crystal structures, **users can further define slip parameters within the respective nested dictionaries, which will be prioritized for adoption**. In above example, the slip system configuration within the "hcp" dictionary will be utilized.
+
+##### 3.1.2.7. Phonon Spectrum
+This function incorporates part of [dflow-phonon](https://github.com/Chengqian-Zhang/dflow-phonon) codes into APEX to make it more complete. This workflow is realized via [Phonopy](https://github.com/phonopy/phonopy), and plus [phonoLAMMPS](https://github.com/abelcarreras/phonolammps) for LAMMPS calculation. 
+
+*IMPORTANT!!*: it should be noticed that one will need the **phonoLAMMPS** package pre-installed within one's `run_image` for proper `LAMMPS` calculation of phonon spectrum.
+
+The parameters related to `Phonon` calculation are listed below:
+  | Key words | Data structure | Default | Description |
+  | :------------ | ----- | ----- | ------------------- |
+  | primitive | Bool | False | Whether to find primitive lattice structure for phonon calculation |
+  | approach | String | "linear" | Specify phonon calculation method when using `VASP`; Two options: 1. "linear" for the *Linear Response Method*, and 2. "displacement" for the *Finite Displacement Method* |
+  | supercell_size | Sequence[Int] | [2, 2, 2] | Size of supercell created for calculation |
+  | MESH | Sequence[Int] | None | Specify the dimensions of the grid in reciprocal space for which the phonon frequencies and eigenvectors are to be calculated. For example: [8, 8, 8]; Refer to [Phonopy MESH](http://phonopy.github.io/phonopy/setting-tags.html#mesh-sampling-tags) |
+  | PRIMITIVE_AXES | String | None | To define the basis vectors of a primitive cell in terms of the basis vectors of a conventional cell for input cell transformation. For example: "0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0"; Refer to [Phonopy PRIMITIVE_AXES](http://phonopy.github.io/phonopy/setting-tags.html#primitive-axes-or-primitive-axis) |
+  | BAND | String | None | Indicate band path in reciprocal space as format of [Phonopy BAND](http://phonopy.github.io/phonopy/setting-tags.html#band-and-band-points); For example: "0 0 0 1/2 0 1/2, 1/2 1/2 1 0 0 0 1/2 1/2 1/2" |
+  | BAND_POINTS | Int | 51 | Number of sampling points including the path ends |
+  | BAND_CONNECTION | Bool | True | With this option, band connections are estimated from eigenvectors and band structure is drawn considering band crossings. In sensitive cases, to obtain better band connections, it requires to increase number of points calculated in band segments by the `BAND_POINTS` tag. |
+
+When utilize the `VASP`, you have **two** primary calculation methods at your disposal: the **Linear Response Method** and the **Finite Displacement Method**.
+
+The **Linear Response Method** has an edge over the Finite Displacement Method in that it eliminates the need for creating super-cells, thereby offering computational efficiency in certain cases. Additionally, this method is particularly well-suited for systems with anomalous phonon dispersion (like systems with Kohn anomalies), as it can precisely calculate the phonons at the specified points.
+
+On the other hand, the **Finite Displacement Method**'s advantage lies in its versatility; it functions as an add-on compatible with any code, including those beyond the scope of density functional theory. The only requirement is that the external code can compute forces. For instance, ABACUS may lack an implementation of the Linear Response Method, but can effectively utilize the Finite Displacement Method implemented in phonon calculation.
+
 
 ### 3.2. Command
 APEX currently supports two seperate run modes: **workflow submission** (running via dflow) and **single-step test** (running without dflow).

apex/__init__.py

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

apex/__main__.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-from .main import main
+from apex.main import main
 
 if __name__ == '__main__':
     main()

apex/core/calculator/ABACUS.py

@@ -4,8 +4,7 @@
 from dpdata import LabeledSystem
 from monty.serialization import dumpfn
 
-import apex.core.calculator.lib.abacus as abacus
-import apex.core.calculator.lib.abacus_scf as abacus_scf
+from apex.core.calculator.lib import abacus_utils, abacus_scf
 #from dpgen import dlog
 from apex.core.calculator.Task import Task
 from apex.utils import sepline
@@ -98,7 +97,7 @@ def make_input_file(self, output_dir, task_type, task_param):
         incar_relax = abacus_scf.get_abacus_input_parameters(relax_incar_path)
 
         # deal with relaxation
-
+        prop_type = task_param.get("type", "relaxation")
         cal_type = task_param["cal_type"]
         cal_setting = task_param["cal_setting"]
 
@@ -158,7 +157,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                 raise RuntimeError("not supported calculation type for ABACUS")
 
             # modify STRU file base on the value of fix_atom
-            abacus.stru_fix_atom(os.path.join(output_dir, "STRU"), fix_atom)
+            abacus_utils.stru_fix_atom(os.path.join(output_dir, "STRU"), fix_atom)
 
         if "basis_type" not in incar:
             logging.info("'basis_type' is not defined, set to be 'pw'!")
@@ -169,7 +168,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                 % incar["basis_type"]
             )
             raise RuntimeError(mess)
-        abacus.write_input(os.path.join(output_dir, "../INPUT"), incar)
+        abacus_utils.write_input(os.path.join(output_dir, "../INPUT"), incar)
         cwd = os.getcwd()
         os.chdir(output_dir)
         if not os.path.islink("INPUT"):
@@ -182,7 +181,7 @@ def make_input_file(self, output_dir, task_type, task_param):
         if "kspacing" in incar:
             kspacing = float(incar["kspacing"])
             if os.path.isfile(os.path.join(output_dir, "STRU")):
-                kpt = abacus.make_kspacing_kpt(
+                kpt = abacus_utils.make_kspacing_kpt(
                     os.path.join(output_dir, "STRU"), kspacing
                 )
                 kpt += [0, 0, 0]
@@ -195,7 +194,7 @@ def make_input_file(self, output_dir, task_type, task_param):
             mess += "You can set key word 'kspacing' (unit in 1/bohr) as a float value in INPUT\n"
             mess += "or set key word 'K_POINTS' as a list in 'cal_setting', e.g. [1,2,3,0,0,0]\n"
             raise RuntimeError(mess)
-        abacus.write_kpt(os.path.join(output_dir, "KPT"), kpt)
+        abacus_utils.write_kpt(os.path.join(output_dir, "KPT"), kpt)
 
     def compute(self, output_dir):
         if not os.path.isfile(os.path.join(output_dir, "INPUT")):

apex/core/calculator/Lammps.py

@@ -4,15 +4,15 @@
 
 from monty.serialization import dumpfn, loadfn
 
-import apex.core.calculator.lib.lammps as lammps
+from apex.core.calculator.lib import lammps_utils
 #from dpgen import dlog
-from apex.core.calculator.lib.lammps import (
+from apex.core.calculator.lib.lammps_utils import (
     inter_deepmd,
     inter_eam_alloy,
     inter_eam_fs,
     inter_meam,
 )
-from apex.core.calculator.Task import Task
+from .Task import Task
 from dflow.python import upload_packages
 upload_packages.append(__file__)
 
@@ -129,10 +129,10 @@ def make_potential_files(self, output_dir):
         dumpfn(self.inter, os.path.join(output_dir, "inter.json"), indent=4)
 
     def make_input_file(self, output_dir, task_type, task_param):
-        lammps.cvt_lammps_conf(
+        lammps_utils.cvt_lammps_conf(
             os.path.join(output_dir, "POSCAR"),
             os.path.join(output_dir, "conf.lmp"),
-            lammps.element_list(self.type_map),
+            lammps_utils.element_list(self.type_map),
         )
 
         # dumpfn(task_param, os.path.join(output_dir, 'task.json'), indent=4)
@@ -195,7 +195,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                 relax_vol = cal_setting["relax_vol"]
 
                 if [relax_pos, relax_shape, relax_vol] == [True, False, False]:
-                    fc = lammps.make_lammps_equi(
+                    fc = lammps_utils.make_lammps_equi(
                         "conf.lmp",
                         self.type_map,
                         self.inter_func,
@@ -207,7 +207,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                         False,
                     )
                 elif [relax_pos, relax_shape, relax_vol] == [True, True, True]:
-                    fc = lammps.make_lammps_equi(
+                    fc = lammps_utils.make_lammps_equi(
                         "conf.lmp",
                         self.type_map,
                         self.inter_func,
@@ -225,7 +225,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                 ] and not task_type == "eos":
                     if "scale2equi" in task_param:
                         scale2equi = task_param["scale2equi"]
-                        fc = lammps.make_lammps_press_relax(
+                        fc = lammps_utils.make_lammps_press_relax(
                             "conf.lmp",
                             self.type_map,
                             scale2equi[int(output_dir[-6:])],
@@ -239,7 +239,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                             maxeval,
                         )
                     else:
-                        fc = lammps.make_lammps_equi(
+                        fc = lammps_utils.make_lammps_equi(
                             "conf.lmp",
                             self.type_map,
                             self.inter_func,
@@ -256,7 +256,7 @@ def make_input_file(self, output_dir, task_type, task_param):
                     False,
                 ] and task_type == "eos":
                     task_param["cal_setting"]["relax_shape"] = False
-                    fc = lammps.make_lammps_equi(
+                    fc = lammps_utils.make_lammps_equi(
                         "conf.lmp",
                         self.type_map,
                         self.inter_func,
@@ -268,15 +268,15 @@ def make_input_file(self, output_dir, task_type, task_param):
                         False,
                     )
                 elif [relax_pos, relax_shape, relax_vol] == [False, False, False]:
-                    fc = lammps.make_lammps_eval(
+                    fc = lammps_utils.make_lammps_eval(
                         "conf.lmp", self.type_map, self.inter_func, self.model_param
                     )
 
                 else:
                     raise RuntimeError("not supported calculation setting for LAMMPS")
 
             elif cal_type == "static":
-                fc = lammps.make_lammps_eval(
+                fc = lammps_utils.make_lammps_eval(
                     "conf.lmp", self.type_map, self.inter_func, self.model_param
                 )
 
@@ -442,7 +442,7 @@ def compute(self, output_dir):
 
             _tmp = self.type_map
             #dlog.debug(_tmp)
-            type_map_list = lammps.element_list(self.type_map)
+            type_map_list = lammps_utils.element_list(self.type_map)
 
             type_map_idx = list(range(len(type_map_list)))
             atom_numbs = []

apex/core/calculator/VASP.py

@@ -6,10 +6,10 @@
 from pymatgen.core.structure import Structure
 from pymatgen.io.vasp import Incar, Kpoints
 
-import apex.core.calculator.lib.vasp as vasp
 #from dpgen import dlog
 from apex.core.calculator.Task import Task
-from apex.core.calculator.lib.vasp import incar_upper
+from apex.core.calculator.lib import vasp_utils
+from apex.core.calculator.lib.vasp_utils import incar_upper
 from apex.utils import sepline
 from dflow.python import upload_packages
 upload_packages.append(__file__)
@@ -81,18 +81,36 @@ def make_input_file(self, output_dir, task_type, task_param):
         incar_relax = incar_upper(Incar.from_file(relax_incar_path))
 
         # deal with relaxation
-
+        prop_type = task_param.get("type", "relaxation")
         cal_type = task_param["cal_type"]
         cal_setting = task_param["cal_setting"]
 
         # user input INCAR for apex calculation
         if "input_prop" in cal_setting and os.path.isfile(cal_setting["input_prop"]):
             incar_prop = os.path.abspath(cal_setting["input_prop"])
             incar = incar_upper(Incar.from_file(incar_prop))
+            logging.info(f"Will use user specified INCAR (path: {incar_prop}) for {prop_type} calculation")
 
         # revise INCAR based on the INCAR provided in the "interaction"
         else:
-            incar = incar_relax
+            if prop_type == "phonon":
+                approach = task_param.get("approach", "linear")
+                logging.info(f"No specification of INCAR for {prop_type} calculation, will auto generate")
+                if approach == "linear":
+                    incar = incar_upper(Incar.from_str(
+                        vasp_utils.make_vasp_phonon_dfpt_incar(
+                            ecut=650, ediff=0.0000001, npar=None, kpar=None
+                        )))
+                elif approach == "displacement":
+                    incar = incar_upper(Incar.from_str(
+                        vasp_utils.make_vasp_static_incar(
+                            ecut=650, ediff=0.0000001, npar=8, kpar=1
+                        )))
+            else:
+                if not prop_type == "relaxation":
+                    logging.info(f"No specification of INCAR for {prop_type} calculation, will use INCAR for relaxation")
+                incar = incar_relax
+
             if cal_type == "relaxation":
                 relax_pos = cal_setting["relax_pos"]
                 relax_shape = cal_setting["relax_shape"]
@@ -192,7 +210,7 @@ def make_input_file(self, output_dir, task_type, task_param):
             os.remove("INCAR")
             os.symlink("../INCAR", "INCAR")
         os.chdir(cwd)
-        ret = vasp.make_kspacing_kpoints(self.path_to_poscar, kspacing, kgamma)
+        ret = vasp_utils.make_kspacing_kpoints(self.path_to_poscar, kspacing, kgamma)
         kp = Kpoints.from_str(ret)
         kp.write_file(os.path.join(output_dir, "KPOINTS"))
 

apex/core/calculator/calculator.py

@@ -4,6 +4,7 @@
 from dflow.python import upload_packages
 upload_packages.append(__file__)
 
+LAMMPS_TYPE = ["deepmd", "meam", "eam_fs", "eam_alloy"]
 
 def make_calculator(inter_parameter, path_to_poscar):
     """
@@ -14,7 +15,7 @@ def make_calculator(inter_parameter, path_to_poscar):
         return VASP(inter_parameter, path_to_poscar)
     elif inter_type == "abacus":
         return ABACUS(inter_parameter, path_to_poscar)
-    elif inter_type in ["deepmd", "meam", "eam_fs", "eam_alloy"]:
+    elif inter_type in LAMMPS_TYPE:
         return Lammps(inter_parameter, path_to_poscar)
     #    if inter_type == 'siesta':
     #        return Siesta(inter_parameter, path_to_poscar)

apex/core/calculator/lib/abacus_scf.py

@@ -3,7 +3,7 @@
 import numpy as np
 from dpdata.abacus.scf import get_cell, get_coords, get_nele_from_stru
 
-from apex.core.calculator.lib import vasp
+from . import vasp_utils
 from dflow.python import upload_packages
 from dargs.dargs import Argument
 upload_packages.append(__file__)
@@ -423,7 +423,7 @@ def make_kspacing_kpoints_stru(stru, kspacing):
     if type(kspacing) is not list:
         kspacing = [kspacing, kspacing, kspacing]
     box = stru["cells"]
-    rbox = vasp.reciprocal_box(box)
+    rbox = vasp_utils.reciprocal_box(box)
     kpoints = [
         max(1, (np.ceil(2 * np.pi * np.linalg.norm(ii) / ks).astype(int)))
         for ii, ks in zip(rbox, kspacing)

apex/core/calculator/lib/abacus.py → apex/core/calculator/lib/abacus_utils.py

@@ -1,11 +1,12 @@
 #!/usr/bin/python3
-import os
+import logging
+import os, re
 
 import dpdata
 import numpy as np
 from pymatgen.core.structure import Structure
 
-import apex.core.calculator.lib.abacus_scf as abacus_scf
+from . import abacus_scf
 from dflow.python import upload_packages
 upload_packages.append(__file__)
 
@@ -499,3 +500,22 @@ def modify_stru_path(strucf, tpath):
     with open(strucf, "w") as f1:
         f1.writelines(lines)
 
+
+def append_orb_file_to_stru(stru, orb: dict = None, prefix: str = ""):
+    with open(stru, 'r') as fin:
+        all_string = fin.read()
+
+    if not re.search("NUMERICAL_ORBITAL", all_string) and orb:
+        logging.info(msg="Append user specific orbital files to STRU")
+        orbital_string = "NUMERICAL_ORBITAL\n"
+        for key, item in orb.items():
+            orbital_string += str(os.path.join(prefix, item))
+            orbital_string += '\n'
+        all_string_new = orbital_string + "\n" + all_string
+
+        if os.path.isfile(stru) or os.path.islink(stru):
+            os.remove(stru)
+
+        with open(stru, 'w') as fout:
+            fout.write(all_string_new)
+

apex/core/calculator/lib/lammps.py → apex/core/calculator/lib/lammps_utils.py

@@ -9,7 +9,7 @@
 from dpdata.periodic_table import Element
 from packaging.version import Version
 
-import apex.core.lib.util as util
+from apex.core.lib import util
 from dflow.python import upload_packages
 upload_packages.append(__file__)