Repository: SPARC-X/SPARC-X-API (language: python)

Issue:
Add contributor template
- [] template for contribution
- [] issue template
- [] Reame section for contributors
various fixes
fixes to #24 #19 #22 and minor bugs


Discussion / hints:
issue and PR templates require maintainer privilege. Will add once v0.2 branch is merged


Fix patch:
diff --git a/.gitignore b/.gitignore
index 1ad89ca5..2a602edb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -771,3 +771,5 @@ ex0-*/
 al-eos-sparc.traj
 */ex1-sparc/
 examples/ex1-ase/
+/SPARC-master/
+/master.zip
diff --git a/README.md b/README.md
index dba1ac1f..2f8acb06 100644
--- a/README.md
+++ b/README.md
@@ -1,9 +1,9 @@
-# `sparc-dft-api`: A Python API for the SPARC DFT Code
+# `SPARC-X-API`: A Python API for the SPARC-X DFT Code
 [![Package](https://raw.githubusercontent.com/alchem0x2A/sparc-dft-api/badges/badges/package.svg)](https://raw.githubusercontent.com/alchem0x2A/sparc-dft-api/badges/badges/package.svg)
 [![Coverage](https://raw.githubusercontent.com/alchem0x2A/sparc-dft-api/badges/badges/coverage.svg)](https://raw.githubusercontent.com/alchem0x2A/sparc-dft-api/badges/badges/coverage.svg)
 [![Unit tests](https://github.com/alchem0x2A/sparc-dft-api/actions/workflows/installation_test.yml/badge.svg)](https://github.com/alchem0x2A/sparc-dft-api/actions/workflows/installation_test.yml)
 
-`sparc-dft-api` is an [ASE](https://wiki.fysik.dtu.dk/ase/)-compatible Python API for the density functional theory (DFT) code [SPARC](https://github.com/SPARC-X/SPARC). It offers:
+`SPARC-X-API` is an [ASE](https://wiki.fysik.dtu.dk/ase/)-compatible Python API for the density functional theory (DFT) code [SPARC](https://github.com/SPARC-X/SPARC). It offers:
 
 1. ASE-compatible I/O format for SPARC files
 2. A JSON API interfacing with SPARC's C-code for parameter validation and conversion
@@ -23,7 +23,7 @@ which includes the pseudopotential files:
 # Change 'sparc-env' to your desired name if needed
 conda create -n sparc-env
 conda activate sparc-env
-conda install -c alchem0x2a sparc-dft-api
+conda install -c alchem0x2a sparc-x-api
 ```
 
 On Linux platforms (x86_64, aarch64), you can also install the
@@ -34,13 +34,13 @@ conda install -c alchem0x2a sparc
 conda activate sparc-env   # Re-activate to have the env variables effective
 ```
 
-*Note: Packaging of sparc-dft-api on conda-forge is in progress.*
+*Note: Packaging of SPARC-X-API on conda-forge is in progress.*
 
 ### 2. Manual installation from source with `pip`
 
 
 ```bash
-python -m pip install git+https://github.com/SPARC-X/sparc-dft-api
+python -m pip install git+https://github.com/SPARC-X/SPARC-X-API
 ```
 
 Optionally, you can download the latest SPMS pseudopotentials and unpacks the pseudopotential files into `<python-lib-root>/site-packages/sparc/psp`:
@@ -64,7 +64,7 @@ python -m sparc.quicktest
 
 A proper setup will display the following sections at the output's conclusion:
 
-<img width="500" alt="image" src="https://github.com/alchem0x2A/sparc-dft-api/assets/6829706/95cb712e-4c77-4b14-8130-4961e3c50278">
+<img width="500" alt="image" src="https://github.com/alchem0x2A/SPARC-X-API/assets/6829706/95cb712e-4c77-4b14-8130-4961e3c50278">
 
 For using the API to parse SPARC input and output files, it's
 essential that the "Import" and "JSON API" tests are successful. For
@@ -76,7 +76,7 @@ configuring the environment variables. If you run into further problems, consult
 [Trouble Shooting](doc/troubleshooting.md).
 
 ## Setting up the environment
-`sparc-dft-api` is designed to automate the discovery of
+`SPARC-X-API` is designed to automate the discovery of
 pseudopotential files, the JSON API, and the SPARC binary. However,
 you can exert fine-grained control over their setup:
 
@@ -87,7 +87,7 @@ order:
 1) Via the `psp_dir` argument passed to the `sparc.SPARC` calculator.
 2) Through the environment variables `$SPARC_PSP_PATH` or `$SPARC_PP_PATH` (this is the
  method employed by [`conda` installation](#1-via-anaconda-or-miniconda-recommended)).
-3) By using `psp8` files bundled with the sparc-dft-api installation (see the
+3) By using `psp8` files bundled with the SPARC-X-API installation (see the
 [manual installation](#2-manual-installation-from-source-with-pip)).
 
 To specify a custom path for your psp8 files, set the `$SPARC_PSP_PATH` or `$SPARC_PP_PATH` variable as follows:
@@ -101,7 +101,7 @@ python -c "from sparc.common import psp_dir; print(psp_dir)"
 ```
 
 ### B) JSON schema
-`sparc-dft-api` is engineered for compatibility with the SPARC
+`SPARC-X-API` is engineered for compatibility with the SPARC
 C-code. It achieves this by loading a JSON schema for
 parameter validation and unit conversion. You can review the default
 schema used by the API at sparc.sparc_json_api.default_json_api
@@ -121,7 +121,7 @@ schema used by the API at sparc.sparc_json_api.default_json_api
 ```
 
 The schema file is generated from SPARC's LaTeX documentation.  In
-upcoming releases of `sparc-dft-api`, we're aiming to provide users
+upcoming releases of `SPARC-X-API`, we're aiming to provide users
 the flexibility to use their own custom schema files. This would be
 particularly useful for those who might be testing a development
 branch of SPARC.
@@ -132,7 +132,7 @@ The command to execute SPARC calculations is determined based on the following p
 
 1) The command argument provided directly to the `sparc.SPARC` calculator.
 2) The environment variable `$ASE_SPARC_COMMAND`
-3) If neither of the above is defined, `sparc-dft-api` looks for the SPARC binary under current `$PATH` and combine with the suitable `mpi` command prefix.
+3) If neither of the above is defined, `SPARC-X-API` looks for the SPARC binary under current `$PATH` and combine with the suitable `mpi` command prefix.
 
 Example:
 ```bash
@@ -146,7 +146,7 @@ export ASE_SPARC_COMMAND="mpirun -n 8 /path/to/sparc -name PREFIX"
 ### 1. Read / write SPARC files
 
 In contrast to many other DFT codes, where the ASE I/O formats refer
-to a single file, `sparc-dft-api` operates on the whole calculation
+to a single file, `SPARC-X-API` operates on the whole calculation
 directory, also known as a "SPARC bundle". This API integrates
 seamlessly with ASE, allowing for the automatic detection of the SPARC
 file format:
@@ -177,9 +177,9 @@ For a deeper dive into the bundle I/O format, see [Advanced Topics](doc/advanced
 A recurring challenge of Python interfaces to DFT codes it the
 inconsistencies between low-level codes (Fortran/C/C++) and outdated
 upper-level APIs regarding parameter sets and default values. To
-address this issue, `sparc-dft-api` handles DFT parameters through a
+address this issue, `SPARC-X-API` handles DFT parameters through a
 JSON schema translated from SPARC's LaTeX documentation.  Each release
-of `sparc-dft-api` is linked with a specific version of the SPARC
+of `SPARC-X-API` is linked with a specific version of the SPARC
 source code, ensuring compatibility and consistency with the default
 parameter set. The main driver of this feature is the
 `sparc.api.SparcAPI` class.
@@ -198,7 +198,7 @@ Topics](doc/advanced_topics.md).
 
 ### 3. Calculator interface
 
-`sparc-dft-api` offers a calculator interface that aligns with many
+`SPARC-X-API` offers a calculator interface that aligns with many
 other ASE calculators.  If you've worked with ASE modules like `Vasp`,
 `QuantumEspresso`, or `GPAW`, you'll find this package intuitive,
 as shown in the following examples:
@@ -252,7 +252,7 @@ If you want to extract more information about the MD simulation steps, take a lo
 
 4. Geometric optimization using ASE's optimizers
 
-The power of `sparc-dft-api` is to combine single point `SPARC` calculations with advanced ASE optimizers, such as BFGS, FIRE or GPMin. Example 2 can be re-written as:
+The power of `SPARC-X-API` is to combine single point `SPARC` calculations with advanced ASE optimizers, such as BFGS, FIRE or GPMin. Example 2 can be re-written as:
 
 ```python
 from sparc.calculator import SPARC
@@ -267,7 +267,7 @@ opt.run(fmax=0.02)
 
 ### 4. Command-line tools
 
-`sparc-dft-api` provides a simple command wrapper `sparc-ase` to add
+`SPARC-X-API` provides a simple command wrapper `sparc-ase` to add
 support of SPARC file formats to the `ase` cli tools. Simple
 replace `ase [subcommand] [args]` with `sparc-ase [subcommand] [args]`
 to access your SPARC bundle files as you would use for other file
@@ -279,9 +279,9 @@ images.
 Below is a screenshot showing the usage of `sparc-ase gui` to visualize a
 short [MD trajectory](tests/outputs/NH3_sort_lbfgs_opt.sparc).
 
-<img width="1200" alt="image" src="https://github.com/alchem0x2A/sparc-dft-api/assets/6829706/e72329ff-7194-4819-94f8-486ef2218844">
+<img width="1200" alt="image" src="https://github.com/alchem0x2A/SPARC-X-API/assets/6829706/e72329ff-7194-4819-94f8-486ef2218844">
 
-### 5. Units used in `sparc-dft-api`
+### 5. Units used in `SPARC-X-API`
 
 In the SPARC DFT code, all input parameters conventionally employ atomic units, such as Hartree and Bohr. Conversely, ASE objects (like `Atoms.positions`, `Atoms.cell`, `Atoms.get_potential_energy()`) utilize eV/Angstrom units.
 
@@ -291,7 +291,63 @@ atoms.calc = SPARC(h=0.25, REFERENCE_CUTOFF=0.5, EXX_RANGE_PBE=0.16, **params)
 ```
 inputs following ASE's convention (e.g., `h`) adopt eV/Angstrom units (thus the same setting can be applied to other DFT calculators),
 On the other hand, all SPARC-specific parameters, which can often be recognized by their capitalized format (like `REFERENCE_CUTOFF`, `EXX_RANGE_PBE`), retain their original values consistent with their representation in the `.inpt` files.
-The reasoning and details about unit conversion can be found in the [Rules for Input Parameters](https://github.com/alchem0x2A/sparc-dft-api/blob/master/doc/advanced_topics.md#rules-for-input-parameters-in-sparcsparc-calculator)  in Advanced Topics.
+The reasoning and details about unit conversion can be found in the [Rules for Input Parameters](https://github.com/alchem0x2A/SPARC-X-API/blob/master/doc/advanced_topics.md#rules-for-input-parameters-in-sparcsparc-calculator)  in Advanced Topics.
+
+
+In order for `SPARC-X-API` to be compatible with other ASE-based DFT calculators,
+there is a list of special parameters consistent with the ASE convention and uses Å / eV / GPa / fs 
+unit system:
+
+| parameter name | meaning                         | example        | equivalent `SPARC` input |
+|----------------|---------------------------------|----------------|--------------------------|
+| `xc`           | Exchange-correlation functional | `xc=pbe` | `EXCHANGE_CORRELATION: GGA_PBE` |
+|                |                                 | `xc=lda` | `EXCHANGE_CORRELATION: LDA_PZ` |
+|                |                                 | `xc=rpbe` | `EXCHANGE_CORRELATION: GGA_RPBE` |
+|                |                                 | `xc=pbesol` | `EXCHANGE_CORRELATION: GGA_PBEsol` |
+|                |                                 | `xc=pbe0` | `EXCHANGE_CORRELATION: PBE0` |
+|                |                                 | `xc=hf` | `EXCHANGE_CORRELATION: HF` |
+|                |                                 | `xc=hse` or `xc=hse03` | `EXCHANGE_CORRELATION: HSE` |
+|                |                                 | `xc=vdwdf1` or `xc=vdw-df` | `EXCHANGE_CORRELATION: vdWDF1` |
+|                |                                 | `xc=vdwdf2` or `xc=vdw-df2` | `EXCHANGE_CORRELATION: vdWDF2` |
+|                |                                 | `xc=scan` | `EXCHANGE_CORRELATION: SCAN` |
+| `h`            | Real grid spacing    (Å)           | `h=0.2`        | `FD_GRID: Nx Ny Nz`  (nearest int value)         |
+| `gpts`         | Explicit grid points |   `gpts=[10, 10, 10]` | `FD_GRID: 10 10 10` |
+| `kpts`         | Kpoint mesh          |   `kpts=[3, 3, 3]`    | `KPOINT_GRID: 3 3 3` |
+| `convergence`  | Dict of convergence criteria (see below) |  | |
+|                | `energy`  eV/atom         | `convergence={"energy": 1e-4}` | `TOL_SCF: 3e-6` |
+|                | `relax` (forces)  eV/Å            | `convergence={"relax": 1e-2}` | `TOL_RELAX: 2e-4` |
+|                | `density` e/atom          | `convergence={`density`: 1e-6}`| `TOL_PSEUDOCHARGE: 1e-6` |
+
+Users from other DFT codes can easily port their ASE codes to `SPARC-X-API` using the special parameters with minimal modification:
+
+Example 1: VASP vs SPARC
+
+```python
+# Using VASP
+from ase.calculators.vasp import Vasp
+calc = Vasp(xc="rpbe", kpts=(9, 9, 9), directory="vasp-calc")
+```
+vs
+```python
+# Using SPARC
+from sparc.calculator import SPARC
+calc = SPARC(xc="rpbe", kpts=(9, 9, 9), directory="sparc-calc.sparc")
+```
+
+Example 2: GPAW (another real-space DFT code) vs SPARC
+```python
+# Using GPAW
+from gpaw import GPAW
+calc = GPAW(xc="PBE", kpts=(9, 9, 9), h=0.25, directory="vasp-calc", convergence={"energy": 1.e-4})
+```
+vs
+```python
+# Using SPARC
+from sparc.calculator import SPARC
+calc = SPARC(xc="PBE", kpts=(9, 9, 9), h=0.25, directory="sparc-calc.sparc", convergence={"energy": 1.e-4})
+```
+
+
 
 ## Troubleshooting
 Please refer to the [troubleshooting](doc/troubleshooting.md) guidelines
@@ -300,7 +356,7 @@ Please refer to the [troubleshooting](doc/troubleshooting.md) guidelines
 A detailed description about how the API works can be found [here](doc/advanced_topics.md)
 
 ## API changes
-The API changes compared to the older release ([v0.1](https://github.com/SPARC-X/sparc-dft-api/tree/eac557f214b402122a506f88f38c7a8767283503)) are summarized [here](doc/changes_v0.1.md)
+The API changes compared to the older release ([v0.1](https://github.com/SPARC-X/SPARC-X-API/tree/eac557f214b402122a506f88f38c7a8767283503)) are summarized [here](doc/changes_v0.1.md)
 
 ## How to contribute
 Please refer to our [guidelines for contributors](doc/contribution_guideline.md)
diff --git a/doc/advanced_topics.md b/doc/advanced_topics.md
index 2a3612a0..347bf04b 100644
--- a/doc/advanced_topics.md
+++ b/doc/advanced_topics.md
@@ -93,24 +93,6 @@ old_atoms.get_potential_energy()
 ```
 
 
-### Special inputs for `sparc.SPARC` calculator
-
-The following input parameters have special meaning in `sparc.SPARC` calculator,
-they are consistent with definitions in other ASE calculators and uses Å / eV / GPa / fs 
-unit system:
-
-| parameter name | meaning                         | example        | equivalent `SPARC` input |
-|----------------|---------------------------------|----------------|--------------------------|
-| `xc`           | Exchange-correlation functional | `xc=pbe` | `EXCHANGE_CORRELATION: GGA_PBE` |
-| `h`            | Real grid spacing    (Å)           | `h=0.2`        | `FD_GRID: Nx Ny Nz`  (calculated values)         |
-| `gpts`         | Explicit grid points |   `gpts=[10, 10, 10]` | `FD_GRID: 10 10 10` |
-| `kpts`         | Kpoint mesh          |   `kpts=[3, 3, 3]`    | `KPOINT_GRID: 3 3 3` |
-| `convergence`  | Dict of convergence criteria (see below) |  | |
-|                | `energy`  eV/atom         | `convergence={"energy": 1e-4}` | `SCF_ENERGY_ACC: 3e-6` |
-|                | `forces`  eV/Å            | `convergence={"forces": 1e-2}` | `TOL_RELAX: 2e-4` |
-|                | `density` e/atom          | `convergence={`density`: 1e-6}`| `TOL_PSEUDOCHARGE: 1e-6` |
-
-*WIP*: support more advanced settings like D3, HSE and DFT+U in `xc` settings
 
 ### Rules for input parameters in `sparc.SPARC` calculator
 
diff --git a/doc/changes_v0.1.md b/doc/changes_v0.1.md
index a63a88d0..07ca9741 100644
--- a/doc/changes_v0.1.md
+++ b/doc/changes_v0.1.md
@@ -1,12 +1,12 @@
 ## Major changes from `sparc-dft-api` [v0.1](https://github.com/SPARC-X/sparc-dft-api/tree/eac557f214b402122a506f88f38c7a8767283503)
 
-`sparc-dft-api` has been heavily refactored in v0.2. If you're using legacy Python codes
+`sparc-dft-api` has been heavily refactored in v1.0. If you're using legacy Python codes
 that are written under v0.1 API, there are a few major changes that require your attention:
 
 1. Support for single `.ion` file format is deprecated. Instead, `v0.2` API treats the whole SPARC directory as a bundle format. Please use `read_sparc` and `write_sparc` methods for basic file I/O instead.
 Nevertheless, reading calculation results generated by a v0.1 API code will not be affected.
 
-2. v0.2 API uses a different mapping scheme for the sorting of ASE atoms objects (similar to `Vasp`), add a comment section in `.ion` file similar to follows:
+2. v1.0 API uses a different mapping scheme for the sorting of ASE atoms objects (similar to `Vasp`), add a comment section in `.ion` file similar to follows:
 ```python
 # ASE-SORT:
 # 3 2 1 0
@@ -14,7 +14,7 @@ Nevertheless, reading calculation results generated by a v0.1 API code will not
 ```
 which maps atoms 3, 2, 1, 0 from the SPARC .ion file order to atoms 0, 1, 2, 3 in ASE order. This is useful for systems that are constructed by ASE's `add_adsorbate` method.
 
-3. v0.2 API accepts all SPARC internal parameters (i.e. **CAPITALIZED**) in *atomic units* for consistency reason. 
+3. v1.0 API accepts all SPARC internal parameters (i.e. **CAPITALIZED**) in *atomic units* for consistency reason. 
 However, we also keep a list of "special input params" that are conventionally used in other ASE calculators, that use Å / eV / GPa / fs unit system.
 
 4. Defining `LATVEC`, `LATVEC_SCALE`, or `CELL` via the calculator parameters is no longer encouraged. Instead, all structure changes should be made to the `Atoms` object.
@@ -24,7 +24,7 @@ For more discussion please see [Advanced Topic] section.
 Below are a list of v0.1 method of the `SPARC` calculator and their current status in v0.2 API. 
 `calc` is an instance of `sparc.SPARC`.
 
-|  old methods           | status in v0.2 API |     alternatives                   |
+|  old methods           | status in v1.0 API |     alternatives                   |
 |------------------------|--------------------|------------------------------------|
 | `interpret_grid_input` | <span style="color:red">deprecated</span>         | `calc.set(fd_grid=[20, 20, 20])`   |
 | `interpret_kpoint_input` | <span style="color:red">deprecated</span>         | `calc.set(kpts=[3, 3, 3])`   |
@@ -47,4 +47,4 @@ Below are a list of v0.1 method of the `SPARC` calculator and their current stat
 | `parse_input_args` | <span style="color:red">deprecated</span> | `calc.set(**kwargs)` |
 | `recover_index_order_from_ion_file` | <span style="color:red">deprecated</span> | Use `calc.sort` and `calc.resort` |
 | `atoms_dict` | <span style="color:red">deprecated</span> | Use third party library like `bson` |
-| `dict_atoms` | <span style="color:red">deprecated</span> | Use third party library like `bson` |
\ No newline at end of file
+| `dict_atoms` | <span style="color:red">deprecated</span> | Use third party library like `bson` |
diff --git a/examples/ex0-eos.py b/examples/ex0-eos.py
index 92c39897..4ecf0989 100644
--- a/examples/ex0-eos.py
+++ b/examples/ex0-eos.py
@@ -15,9 +15,7 @@ def main():
     calc = SPARC(h=0.25, kpts=(3, 3, 3), xc="pbe", directory="ex0-eos")
     vol = atoms.get_volume()
     atoms.calc = calc
-    eos = calculate_eos(
-        atoms, npoints=5, eps=0.05, trajectory="al-eos-sparc.traj"
-    )
+    eos = calculate_eos(atoms, npoints=5, eps=0.05, trajectory="al-eos-sparc.traj")
     print("Original volume: Ang^3", vol)
     v, e, B = eos.fit()
     print("Fitted volume (Ang^3), energy (eV), modulus (eV/Ang^3)")
diff --git a/examples/ex1-ase-optimize.py b/examples/ex1-ase-optimize.py
index 58e3b9eb..69289914 100644
--- a/examples/ex1-ase-optimize.py
+++ b/examples/ex1-ase-optimize.py
@@ -50,7 +50,7 @@ def optimize_ase_lbfgs():
     )
     atoms.calc = calc
     opt = BFGS(atoms)
-    #breakpoint()
+    # breakpoint()
     opt.run(fmax=0.02)
     e_fin = atoms.get_potential_energy()
     f_fin = atoms.get_forces()
diff --git a/setup.py b/setup.py
index 553e7b80..f7911f68 100644
--- a/setup.py
+++ b/setup.py
@@ -17,13 +17,13 @@
 ]
 
 setup(
-    name="sparc-dft-api",
+    name="sparc-x-api",
     version="1.0.0",
     python_requires=">=3.8",
     description="Python API for the SPARC DFT Code",
     author="Tian Tian, Ben Comer",
     author_email="alchem0x2a@gmail.com, ben.comer@gatech.edu",
-    url="https://github.com/SPARC-X/sparc-dft-api",
+    url="https://github.com/SPARC-X/SPARC-X-API",
     packages=find_packages(),
     install_requires=["ase>=3.22.0"],
     entry_points={
diff --git a/sparc/api.py b/sparc/api.py
index 5f3f9215..06e141ae 100644
--- a/sparc/api.py
+++ b/sparc/api.py
@@ -30,8 +30,7 @@ def get_parameter_dict(self, parameter):
         parameter = parameter.upper()
         if parameter not in self.parameters.keys():
             raise KeyError(
-                f"Parameter {parameter} is not known to "
-                f"SPARC {self.sparc_version}!"
+                f"Parameter {parameter} is not known to " f"SPARC {self.sparc_version}!"
             )
         return self.parameters[parameter]
 
@@ -90,8 +89,13 @@ def validate_input(self, parameter, input):
                 float(input)
                 return True
             except (TypeError, ValueError):
-                return False
+                try:
+                    float(input.split()[0])
+                    return True
+                except Exception:
+                    return False
         elif "array" in dtype:
+            # import pdb; pdb.set_trace()
             if is_input_string:
                 if ("." in input) and ("integer" in dtype):
                     warn(
@@ -103,7 +107,8 @@ def validate_input(self, parameter, input):
                         )
                     )
                 try:
-                    arr = np.genfromtxt(input.splitlines(), dtype=float)
+                    # import pdb; pdb.set_trace()
+                    arr = np.genfromtxt(input.splitlines(), dtype=float, ndmin=1)
                     # In valid input with nan
                     if np.isnan(arr).any():
                         arr = np.array(0.0)
@@ -111,7 +116,7 @@ def validate_input(self, parameter, input):
                     arr = np.array(0.0)
             else:
                 try:
-                    arr = np.asarray(input)
+                    arr = np.atleast_1d(np.asarray(input))
                     if (arr.dtype not in (int, bool)) and ("integer" in dtype):
                         warn(
                             (
@@ -124,6 +129,9 @@ def validate_input(self, parameter, input):
                 except Exception:
                     arr = np.array(0.0)
             return len(arr.shape) > 0
+        # elif dtype == "other":
+        #     # Any "other"-type inputs should be provided only using string
+        #     return is_input_string
         else:
             raise ValueError(f"Data type {dtype} is not supported!")
 
@@ -132,6 +140,8 @@ def convert_string_to_value(self, parameter, string):
 
         # Special case, the string may be a multiline string-array!
         if isinstance(string, list):
+            # Make sure there is a line break at the end, for cases like ["2."]
+            string.append("")
             string = [s.strip() for s in string]
             string = "\n".join(string)
 
@@ -153,15 +163,25 @@ def convert_string_to_value(self, parameter, string):
             if allow_bool_input:
                 value = bool(value)
         elif dtype == "double":
-            value = float(string)
+            # Some inputs, like TARGET_PRESSURE, may be accepted with a unit
+            # like 0.0 GPa. Only accept the first part
+            try:
+                value = float(string)
+            except ValueError as e:
+                try:
+                    value = float(string.split()[0])
+                except Exception:
+                    raise e
         elif dtype == "integer array":
-            value = np.genfromtxt(string.splitlines(), dtype=int)
+            value = np.genfromtxt(string.splitlines(), dtype=int, ndmin=1)
             if allow_bool_input:
                 value = value.astype(bool)
         elif dtype == "double array":
-            value = np.genfromtxt(string.splitlines(), dtype=float)
-        else:
+            value = np.genfromtxt(string.splitlines(), dtype=float, ndmin=1)
+        elif dtype == "other":
+            value = string
             # should not happen since validate_input has gatekeeping
+        else:
             raise ValueError(f"Unsupported type {dtype}")
 
         return value
@@ -190,6 +210,10 @@ def convert_value_to_string(self, parameter, value):
             string = "{:.14f}".format(float(value))
         elif dtype in ("integer array", "double array"):
             string = _array_to_string(value, dtype)
+        elif dtype == "other":
+            if not is_input_string:
+                raise ValueError("Only support string value when datatype is other")
+            string = value
         else:
             # should not happen since validate_input has gatekeeping
             raise ValueError(f"Unsupported type {dtype}")
@@ -206,9 +230,7 @@ def _array_to_string(arr, format):
         fmt = "%d"
     elif format in ("double array", "double"):
         fmt = "%.14f"
-    np.savetxt(
-        buf, arr, delimiter=" ", fmt=fmt, header="", footer="", newline="\n"
-    )
+    np.savetxt(buf, arr, delimiter=" ", fmt=fmt, header="", footer="", newline="\n")
     # Return the string output of the buffer with
     # whitespaces removed
     return buf.getvalue().strip()
diff --git a/sparc/calculator.py b/sparc/calculator.py
index 6c72f3fa..aee2a609 100644
--- a/sparc/calculator.py
+++ b/sparc/calculator.py
@@ -5,7 +5,7 @@
 import subprocess
 
 from .io import SparcBundle
-from .utils import _find_default_sparc, h2gpts
+from .utils import _find_default_sparc, h2gpts, deprecated
 from warnings import warn, warn_explicit
 from .api import SparcAPI
 import datetime
@@ -27,22 +27,6 @@
 defaultAPI = SparcAPI()
 
 
-def deprecated(message):
-    def decorator(func):
-        def new_func(*args, **kwargs):
-            warn(
-                "Function {} is deprecated sparc-dft-api >= v0.2! {}".format(
-                    func.__name__, message
-                ),
-                category=DeprecationWarning,
-            )
-            return func(*args, **kwargs)
-
-        return new_func
-
-    return decorator
-
-
 class SPARC(FileIOCalculator):
     # TODO: magmom should be a possible input
     implemented_properties = ["energy", "forces", "fermi", "stress"]
@@ -177,9 +161,7 @@ def check_state(self, atoms, tol=1e-9):
         # A few hard-written rules. Wrapping should only affect the position
         if "positions" in system_changes:
             atoms_copy.wrap()
-            new_system_changes = FileIOCalculator.check_state(
-                self, atoms_copy, tol=tol
-            )
+            new_system_changes = FileIOCalculator.check_state(self, atoms_copy, tol=tol)
             # TODO: make sure such check only happens for PBC
             # the position is wrapped, accept as the same structure
             if "positions" not in new_system_changes:
@@ -220,16 +202,16 @@ def _make_command(self, extras=""):
             self.command = command_env
         return f"{self.command} {extras}"
 
-    def calculate(
-        self, atoms=None, properties=["energy"], system_changes=all_changes
-    ):
+    def calculate(self, atoms=None, properties=["energy"], system_changes=all_changes):
         """Perform a calculation step"""
-        # For v1.0.0, we'll only allow pbc=True to make ourselves easier
-        # TODO: need to have more flexible support for pbc types and check_state
-        if not all(atoms.pbc):
-            raise NotImplementedError(
-                "Non-pbc atoms input has not been tested in the api. Please use pbc=True for now."
-            )
+        # Check if the user accidentally provides atoms unit cell without vacuum
+
+        if atoms and np.any(atoms.cell.cellpar()[:3] == 0):
+            # TODO: choose a better error name
+            msg = "Cannot setup SPARC calculation because at least one of the lattice dimension is zero!"
+            if any([bc_ is False for bc_ in atoms.pbc]):
+                msg += " Please add a vacuum in the non-periodic direction of your input structure."
+            raise ValueError(msg)
         Calculator.calculate(self, atoms, properties, system_changes)
         self.write_input(self.atoms, properties, system_changes)
         self.execute()
@@ -248,7 +230,15 @@ def calculate(
                     self.atoms.get_initial_magnetic_moments()
                 )
 
-            # atoms = self.atoms.copy()
+    def get_stress(self, atoms=None):
+        """Warn user the dimensionality change when using stress"""
+        if "stress_equiv" in self.results:
+            raise NotImplementedError(
+                "You're requesting stress in a low-dimensional system. Please use `calc.results['stress_equiv']` instead!"
+            )
+        return super().get_stress(atoms)
+
+        # atoms = self.atoms.copy()
 
     # def update_atoms(self, atoms):
     #     """Update atoms after calculation if the positions are changed
@@ -284,9 +274,7 @@ def _check_input_exclusion(self, input_parameters, atoms=None):
             )
 
         # Rule 2: LATVEC_SCALE, CELL
-        if ("LATVEC_SCALE" in input_parameters) and (
-            "CELL" in input_parameters
-        ):
+        if ("LATVEC_SCALE" in input_parameters) and ("CELL" in input_parameters):
             # TODO: change to ExclusionParameterError
             raise ValueError(
                 "LATVEC_SCALE and CELL cannot be specified simultaneously!"
@@ -296,12 +284,7 @@ def _check_input_exclusion(self, input_parameters, atoms=None):
         # LATVEC, LATVEC_SCALE or CELL
         # TODO: make sure the rule makes sense for molecules
         if atoms is not None:
-            if any(
-                [
-                    p in input_parameters
-                    for p in ["LATVEC", "LATVEC_SCALE", "CELL"]
-                ]
-            ):
+            if any([p in input_parameters for p in ["LATVEC", "LATVEC_SCALE", "CELL"]]):
                 raise ValueError(
                     "When passing an ase atoms object, LATVEC, LATVEC_SCALE or CELL cannot be set simultaneously!"
                 )
@@ -317,10 +300,7 @@ def _check_minimal_input(self, input_parameters):
                 raise ValueError(f"Parameter {param} is not provided.")
         # At least one from ECUT, MESH_SPACING and FD_GRID must be provided
         if not any(
-            [
-                param in input_parameters
-                for param in ("ECUT", "MESH_SPACING", "FD_GRID")
-            ]
+            [param in input_parameters for param in ("ECUT", "MESH_SPACING", "FD_GRID")]
         ):
             raise ValueError(
                 "You should provide at least one of ECUT, MESH_SPACING or FD_GRID."
@@ -395,10 +375,7 @@ def execute(self):
         errorcode = self.proc.returncode
 
         if errorcode > 0:
-            msg = (
-                f"SPARC failed with command {command}"
-                f"with error code {errorcode}"
-            )
+            msg = f"SPARC failed with command {command}" f"with error code {errorcode}"
             raise RuntimeError(msg)
 
         return
@@ -416,9 +393,7 @@ def read_results(self):
         """Parse from the SparcBundle"""
         # TODO: try use cache?
         # self.sparc_bundle.read_raw_results()
-        last = self.sparc_bundle.convert_to_ase(
-            indices=-1, include_all_files=False
-        )
+        last = self.sparc_bundle.convert_to_ase(indices=-1, include_all_files=False)
         self.atoms = last.copy()
         self.results.update(last.calc.results)
 
@@ -497,13 +472,37 @@ def _convert_special_params(self, atoms=None):
         params = self.special_params.copy()
 
         # xc --> EXCHANGE_CORRELATION
-        # TODO: more XC options
         if "xc" in params:
             xc = params.pop("xc")
             if xc.lower() == "pbe":
                 converted_sparc_params["EXCHANGE_CORRELATION"] = "GGA_PBE"
             elif xc.lower() == "lda":
-                converted_sparc_params["EXCHANGE_CORRELATION"] = "LDA_PW"
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "LDA_PZ"
+            elif xc.lower() == "rpbe":
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "GGA_RPBE"
+            elif xc.lower() == "pbesol":
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "GGA_PBEsol"
+            elif xc.lower() == "pbe0":
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "PBE0"
+            elif xc.lower() == "hf":
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "HF"
+            # backward compatibility for HSE03. Note HSE06 is not supported yet
+            elif (xc.lower() == "hse") or (xc.lower() == "hse03"):
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "HSE"
+            # backward compatibility for VASP-style XCs
+            elif (
+                (xc.lower() == "vdwdf1")
+                or (xc.lower() == "vdw-df")
+                or (xc.lower() == "vdw-df1")
+            ):
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "vdWDF1"
+            elif (xc.lower() == "vdwdf2") or (xc.lower() == "vdw-df2"):
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "vdWDF2"
+            elif xc.lower() == "scan":
+                converted_sparc_params["EXCHANGE_CORRELATION"] = "SCAN"
+            else:
+                # TODO: alternative exception
+                raise ValueError(f"xc keyword value {xc} is invalid!")
 
         # h --> gpts
         if "h" in params:
@@ -517,10 +516,7 @@ def _convert_special_params(self, atoms=None):
                     "Must have an active atoms object to convert h --> gpts!"
                 )
             if any(
-                [
-                    p in self.valid_params
-                    for p in ("FD_GRID", "ECUT", "MESH_SPACING")
-                ]
+                [p in self.valid_params for p in ("FD_GRID", "ECUT", "MESH_SPACING")]
             ):
                 warn(
                     "You have specified one of FD_GRID, ECUT or MESH_SPACING, "
@@ -538,9 +534,7 @@ def _convert_special_params(self, atoms=None):
                 converted_sparc_params["FD_GRID"] = gpts
             else:
                 # TODO: customize error
-                raise ValueError(
-                    f"Input parameter gpts has invalid value {gpts}"
-                )
+                raise ValueError(f"Input parameter gpts has invalid value {gpts}")
 
         # kpts
         if "kpts" in params:
@@ -550,9 +544,7 @@ def _convert_special_params(self, atoms=None):
                 converted_sparc_params["KPOINT_GRID"] = kpts
             else:
                 # TODO: customize error
-                raise ValueError(
-                    f"Input parameter kpts has invalid value {kpts}"
-                )
+                raise ValueError(f"Input parameter kpts has invalid value {kpts}")
 
         # nbands
         if "nbands" in params:
@@ -563,9 +555,7 @@ def _convert_special_params(self, atoms=None):
                 converted_sparc_params["NSTATES"] = nbands
             else:
                 # TODO: customize error
-                raise ValueError(
-                    f"Input parameter nbands has invalid value {nbands}"
-                )
+                raise ValueError(f"Input parameter nbands has invalid value {nbands}")
 
         # convergence is a dict
         if "convergence" in params:
@@ -573,10 +563,10 @@ def _convert_special_params(self, atoms=None):
             tol_e = convergence.get("energy", None)
             if tol_e:
                 # TOL SCF: Ha / atom <--> energy tol: eV / atom
-                converted_sparc_params["SCF_ENERGY_ACC"] = tol_e / Hartree
+                converted_sparc_params["TOL_SCF"] = tol_e / Hartree
 
             # TODO: per AJ's suggestion, better change forces to relaxation
-            tol_f = convergence.get("forces", None)
+            tol_f = convergence.get("relax", None)
             if tol_f:
                 # TOL SCF: Ha / Bohr <--> energy tol: Ha / Bohr
                 converted_sparc_params["TOL_RELAX"] = tol_f / Hartree * Bohr
@@ -620,9 +610,7 @@ def interpret_grid_input(self, atoms, **kwargs):
     def interpret_kpoint_input(self, atoms, **kwargs):
         return None
 
-    @deprecated(
-        "Please use SPARC.set instead for setting downsampling parameter"
-    )
+    @deprecated("Please use SPARC.set instead for setting downsampling parameter")
     def interpret_downsampling_input(self, atoms, **kwargs):
         return None
 
diff --git a/sparc/docparser.py b/sparc/docparser.py
index 868de7e1..b52dac0a 100644
--- a/sparc/docparser.py
+++ b/sparc/docparser.py
@@ -22,7 +22,7 @@ class SPARCDocParser(object):
     def __init__(
         self,
         directory=".",
-        main_file="Manual.tex",
+        main_file="Manual*.tex",
         intro_file="Introduction.tex",
         params_from_intro=True,
         parse_version=True,
@@ -34,7 +34,9 @@ def __init__(
             |---- Manual.tex
                   |---- Introduction.tex
                         |---- {Section}.tex
-        TODO: include the parameters for SQ / HT calculations
+
+        For parameters additional to the standard SPARC options, such as the SQ / cyclix
+        options, we merge the dict from the sub-dirs
 
         Arguments:
         `doc_root`: root directory to the LaTeX files, may look like `SPARC/doc/.LaTeX`
@@ -44,19 +46,24 @@ def __init__(
         `parse_date`: get the SPARC version by date
         """
         self.root = Path(directory)
-        self.main_file = self.root / main_file
-        if not self.main_file.is_file():
-            raise FileNotFoundError(f"Main file {main_file} is missing!")
+        self.main_file = self.find_main_file(main_file)
         self.intro_file = self.root / intro_file
         if not self.intro_file.is_file():
-            raise FileNotFoundError(
-                f"Introduction file {intro_file} is missing!"
-            )
+            raise FileNotFoundError(f"Introduction file {intro_file} is missing!")
         self.include_files = self.get_include_files()
         self.params_from_intro = params_from_intro
         self.parse_version(parse_version)
         self.parse_parameters()
 
+    def find_main_file(self, main_file_pattern):
+        """Find the matching name for the main-file, e.g. Manual.tex or Manual_cyclix.tex"""
+        candidates = list(self.root.glob(main_file_pattern))
+        if len(candidates) != 1:
+            raise FileNotFoundError(
+                f"Main file {main_file_pattern} is missing or more than 1 exists!"
+            )
+        return candidates[0]
+
     def get_include_files(self):
         """Get a list of included LaTeX files from Manual.tex"""
         pattern = r"\\begin\{document\}(.*?)\\end\{document\}"
@@ -101,9 +108,7 @@ def parse_version(self, parse=True):
             self.version = None
             return
         date_str = match[0].strip().replace(",", " ")
-        date_version = datetime.strptime(date_str, "%b %d %Y").strftime(
-            "%Y.%m.%d"
-        )
+        date_version = datetime.strptime(date_str, "%b %d %Y").strftime("%Y.%m.%d")
         self.version = date_version
         return
 
@@ -176,9 +181,9 @@ def __parse_intro_file(self):
         )
         pattern_block = r"\\begin\{block\}\{(.*?)\}([\s\S]*?)\\end\{block\}"
         pattern_line = r"\\hyperlink\{(.*?)\}{\\texttt\{(.*?)\}\}"
-        text_params = re.findall(
-            pattern_params, text_intro, re.DOTALL | re.MULTILINE
-        )[0]
+        text_params = re.findall(pattern_params, text_intro, re.DOTALL | re.MULTILINE)[
+            0
+        ]
         parameter_categories = []
         parameter_dict = {}
         for match in re.findall(pattern_block, text_params):
@@ -202,9 +207,7 @@ def __parse_intro_file(self):
                     label, symbol = match[0].strip(), convert_tex_parameter(
                         match[1].strip()
                     )
-                    parameter_dict[cat].append(
-                        {"label": label, "symbol": symbol}
-                    )
+                    parameter_dict[cat].append({"label": label, "symbol": symbol})
         return parameter_categories, parameter_dict
 
     def __parse_all_included_files(self):
@@ -247,7 +250,7 @@ def parse_parameters(self):
 
         return
 
-    def to_json(self, indent=False):
+    def to_dict(self):
         """Output a json string from current document parser
 
         Arguments:
@@ -260,15 +263,29 @@ def to_json(self, indent=False):
         doc["other_parameters"] = {
             k: v for k, v in sorted(self.other_parameters.items())
         }
-        doc["data_types"] = list(
-            set([p["type"] for p in self.parameters.values()])
-        )
-        json_string = json.dumps(doc, indent=indent)
-        return json_string
+        doc["data_types"] = list(set([p["type"] for p in self.parameters.values()]))
+        # json_string = json.dumps(doc, indent=indent)
+        return doc
 
     @classmethod
-    def from_directory(cls, directory=".", **kwargs):
-        return cls(directory=directory, **kwargs)
+    def json_from_directory(cls, directory=".", include_subdirs=True, **kwargs):
+        """Recursively add parameters from all Manual files"""
+        root_dict = cls(directory=directory, **kwargs).to_dict()
+        if include_subdirs:
+            for sub_manual_tex in directory.glob("*/Manual*.tex"):
+                subdir = sub_manual_tex.parent
+                try:
+                    sub_dict = cls(directory=subdir, parse_version=False).to_dict()
+                except FileNotFoundError:
+                    print(subdir)
+                    continue
+                # We only merge the parameters that have not appeared in the main
+                # manual. TODO: maybe assign repeating parameters to a different section?
+                for param, param_desc in sub_dict["parameters"].items():
+                    if param not in root_dict["parameters"]:
+                        root_dict["parameters"][param] = param_desc
+        json_string = json.dumps(root_dict, indent=True)
+        return json_string
 
 
 def convert_tex_parameter(text):
@@ -323,9 +340,7 @@ def convert_tex_default(text, desired_type=None):
         text = text.replace(m, r)
     text = re.sub(r"\n+", "\n", text)
     # Remove all comment lines
-    text = "\n".join(
-        [l for l in text.splitlines() if not l.lstrip().startswith("%")]
-    )
+    text = "\n".join([l for l in text.splitlines() if not l.lstrip().startswith("%")])
 
     # print(text)
     converted = None
@@ -367,9 +382,7 @@ def convert_comment(text):
         text = text.replace(m, r)
     text = re.sub(r"\n+", "\n", text)
     # Remove all comment lines
-    text = "\n".join(
-        [l for l in text.splitlines() if not l.lstrip().startswith("%")]
-    )
+    text = "\n".join([l for l in text.splitlines() if not l.lstrip().startswith("%")])
     return text
 
 
@@ -468,9 +481,7 @@ def sanitize_default(param_dict):
     sanitized_dict = param_dict.copy()
     original_default = sanitized_dict["default"]
     sanitized_dict["default_remark"] = original_default
-    converted_default = convert_tex_default(
-        original_default, param_dict["type"]
-    )
+    converted_default = convert_tex_default(original_default, param_dict["type"])
     sanitized_dict["default"] = converted_default
     return sanitized_dict
 
@@ -503,9 +514,7 @@ def sanitize_type(param_dict):
         sanitized_type = origin_type
 
     # Pass 2, test if int values are arrays
-    if (origin_type in ["int", "integer", "double"]) and (
-        sanitized_type is None
-    ):
+    if (origin_type in ["int", "integer", "double"]) and (sanitized_type is None):
         if "int" in origin_type:
             origin_type = "integer"
         # Test if the value from example is a single value or array
@@ -513,9 +522,7 @@ def sanitize_type(param_dict):
             example_value = param_dict["example"].split(":")[1]
             default = param_dict["default"]
             _array_test = is_array(example_value)
-            _bool_test = contain_only_bool(example_value) and contain_only_bool(
-                default
-            )
+            _bool_test = contain_only_bool(example_value) and contain_only_bool(default)
         except Exception as e:
             warn(
                 f"Array conversion failed for {example_value}, ignore."
@@ -535,9 +542,7 @@ def sanitize_type(param_dict):
     if sanitized_type is None:
         # Currently there is only one NPT_NH_QMASS has this type
         # TODO: think of a way to format a mixed array?
-        warn(
-            f"Type of {symbol} if not standard digit or array, mark as others."
-        )
+        warn(f"Type of {symbol} if not standard digit or array, mark as others.")
         sanitized_type = "other"
         # TODO: how about provide a true / false type?
     sanitized_dict["type"] = sanitized_type
@@ -555,13 +560,19 @@ def sanitize_type(param_dict):
         default="parameters.json",
         help="Output file name (json-formatted)",
     )
+    argp.add_argument(
+        "--include-subdirs",
+        action="store_true",
+        help="Parse manual parameters from subdirs",
+    )
     argp.add_argument(
         "root", help="Root directory of the latex files"
     )  # root directory of the LaTeX files
     args = argp.parse_args()
     output = Path(args.output).with_suffix(".json")
-    parser = SPARCDocParser(directory=Path(args.root))
-    json_string = parser.to_json(indent=True)
+    json_string = SPARCDocParser.json_from_directory(
+        directory=Path(args.root), include_subdirs=args.include_subdirs
+    )
     with open(output, "w", encoding="utf8") as fd:
         fd.write(json_string)
     print(f"SPARC parameter specifications written to {output}!")
diff --git a/sparc/download_data.py b/sparc/download_data.py
index 363fe0a2..7ccd4025 100644
--- a/sparc/download_data.py
+++ b/sparc/download_data.py
@@ -39,18 +39,14 @@ def download_psp(sparc_tag=sparc_tag, psp_dir=psp_dir):
                 source_dir = next(tmpdir.glob("SPARC-*/psps"))
                 print(f"Found source_dir at {source_dir}")
                 if not source_dir.is_dir():
-                    raise FileNotFoundError(
-                        "Error downloading or extracting zip"
-                    )
+                    raise FileNotFoundError("Error downloading or extracting zip")
                 print(f"Moving psp files to {psp_dir}")
                 for ext in ("*.psp8", "*.psp", "*.pot"):
                     for pspf in source_dir.glob(ext):
                         print(f"Found {pspf} --> {psp_dir}")
                         shutil.copy(pspf, psp_dir)
     if not is_psp_download_complete(psp_dir):
-        raise RuntimeError(
-            f"Files downloaded to {psp_dir} have different checksums!"
-        )
+        raise RuntimeError(f"Files downloaded to {psp_dir} have different checksums!")
     return True
 
 
diff --git a/sparc/io.py b/sparc/io.py
index 693512a4..d630812d 100644
--- a/sparc/io.py
+++ b/sparc/io.py
@@ -27,7 +27,7 @@
 from .sparc_parsers.pseudopotential import copy_psp_file, parse_psp8_header
 from .common import psp_dir as default_psp_dir
 from .download_data import is_psp_download_complete
-from .utils import string2index
+from .utils import string2index, deprecated
 
 from ase.calculators.singlepoint import SinglePointDFTCalculator
 from ase.units import Hartree
@@ -270,9 +270,7 @@ def read_raw_results(self, include_all_files=False):
         """
         # Find the max output index
         # TODO: move this into another function
-        last_out = sorted(
-            self.directory.glob(f"{self.label}.out*"), reverse=True
-        )
+        last_out = sorted(self.directory.glob(f"{self.label}.out*"), reverse=True)
         # No output file, only ion / inpt
         if len(last_out) == 0:
             self.last_image = -1
@@ -343,9 +341,7 @@ def _read_results_from_index(self, index, d_format="{:02d}"):
                 self.sorting = sorting
             else:
                 # Compare stored sorting
-                assert (
-                    tuple(self.sorting["sort"]) == tuple(sorting["sort"])
-                ) and (
+                assert (tuple(self.sorting["sort"]) == tuple(sorting["sort"])) and (
                     tuple(self.sorting["resort"]) == tuple(sorting["resort"])
                 ), "Sorting information changed!"
         return results_dict
@@ -367,17 +363,11 @@ def convert_to_ase(self, index=-1, include_all_files=False, **kwargs):
         # print("RAW RES: ", raw_results)
         for entry in raw_results:
             if "static" in entry:
-                calc_results, images = self._extract_static_results(
-                    entry, index=":"
-                )
+                calc_results, images = self._extract_static_results(entry, index=":")
             elif "geopt" in entry:
-                calc_results, images = self._extract_geopt_results(
-                    entry, index=":"
-                )
+                calc_results, images = self._extract_geopt_results(entry, index=":")
             elif "aimd" in entry:
-                calc_results, images = self._extract_aimd_results(
-                    entry, index=":"
-                )
+                calc_results, images = self._extract_aimd_results(entry, index=":")
             else:
                 calc_results, images = None, [self.init_atoms.copy()]
 
@@ -406,9 +396,7 @@ def _make_singlepoint(self, calc_results, images, raw_results):
             # Simply copy the results?
             sp.results.update(res)
             sp.name = "sparc"
-            sp.kpts = (
-                raw_results["inpt"].get("params", {}).get("KPOINT_GRID", None)
-            )
+            sp.kpts = raw_results["inpt"].get("params", {}).get("KPOINT_GRID", None)
             # There may be a better way handling the parameters...
             sp.parameters = raw_results["inpt"].get("params", {})
             sp.raw_parameters = {
@@ -442,6 +430,9 @@ def _extract_static_results(self, raw_results, index=":"):
         if "stress" in static_results:
             calc_results["stress"] = static_results["stress"]
 
+        if "stress_equiv" in static_results:
+            calc_results["stress_equiv"] = static_results["stress_equiv"]
+
         atoms = self.init_atoms.copy()
         if "atoms" in static_results:
             atoms_dict = static_results["atoms"]
@@ -449,9 +440,7 @@ def _extract_static_results(self, raw_results, index=":"):
             # TODO: Check naming, is it coord_frac or scaled_positions?
             if "coord_frac" in atoms_dict:
                 # TODO: check if set_scaled_positions requires constraint?
-                atoms.set_scaled_positions(
-                    atoms_dict["coord_frac"][self.resort]
-                )
+                atoms.set_scaled_positions(atoms_dict["coord_frac"][self.resort])
             elif "coord" in atoms_dict:
                 atoms.set_positions(
                     atoms_dict["coord"][self.resort], apply_constraint=False
@@ -552,13 +541,11 @@ def _extract_aimd_results(self, raw_results, index=":"):
             partial_result = {}
             atoms = self.init_atoms.copy()
             if "total energy per atom" in result:
-                partial_result["energy"] = result[
-                    "total energy per atom"
-                ] * len(atoms)
+                partial_result["energy"] = result["total energy per atom"] * len(atoms)
             if "free energy per atom" in result:
-                partial_result["free energy"] = result[
-                    "free energy per atom"
-                ] * len(atoms)
+                partial_result["free energy"] = result["free energy per atom"] * len(
+                    atoms
+                )
 
             if "forces" in result:
                 # The forces are already re-sorted!
@@ -566,9 +553,7 @@ def _extract_aimd_results(self, raw_results, index=":"):
 
             # Modify the atoms in-place
             if "positions" not in result:
-                raise ValueError(
-                    "Cannot have aimd without positions information!"
-                )
+                raise ValueError("Cannot have aimd without positions information!")
 
             atoms.set_positions(
                 result["positions"][self.resort], apply_constraint=False
@@ -679,15 +664,73 @@ def write_sparc(filename, images, **kwargs):
         atoms = images
     elif isinstance(images, list):
         if len(images) > 1:
-            raise ValueError(
-                "SPARC format only supports writing one atoms object!"
-            )
+            raise ValueError("SPARC format only supports writing one atoms object!")
         atoms = images[0]
     sb = SparcBundle(directory=filename, mode="w")
     sb._write_ion_and_inpt(atoms, **kwargs)
     return
 
 
+@deprecated(
+    "Reading individual .ion is not recommended. Please use read_sparc instead."
+)
+def read_ion(filename, **kwargs):
+    """Parse an .ion file inside the SPARC bundle using a wrapper around SparcBundle
+    The reader works only when other files (.inpt) exist.
+
+    The returned Atoms object of read_ion method only contains the initial positions
+    """
+    parent_dir = Path(filename).parent
+    sb = SparcBundle(directory=parent_dir)
+    atoms = sb._read_ion_and_inpt()
+    return atoms
+
+
+@deprecated(
+    "Writing individual .ion file is not recommended. Please use write_sparc instead."
+)
+def write_ion(filename, atoms, **kwargs):
+    """Write .ion and .inpt files using the SparcBundle wrapper.
+
+    This is only for backward compatibility
+    """
+    label = Path(filename).with_suffix("").name
+    parent_dir = Path(filename).parent
+    sb = SparcBundle(directory=parent_dir, label=label, mode="w")
+    sb._write_ion_and_inpt(atoms, **kwargs)
+    return atoms
+
+
+def read_static(filename, index=-1, **kwargs):
+    """Parse a .static file bundle using a wrapper around SparcBundle
+    The reader works only when other files (.ion, .inpt) exist.
+    """
+    parent_dir = Path(filename).parent
+    sb = SparcBundle(directory=parent_dir)
+    atoms_or_images = sb.convert_to_ase(index=index, **kwargs)
+    return atoms_or_images
+
+
+def read_geopt(filename, index=-1, **kwargs):
+    """Parse a .geopt file bundle using a wrapper around SparcBundle
+    The reader works only when other files (.ion, .inpt) exist.
+    """
+    parent_dir = Path(filename).parent
+    sb = SparcBundle(directory=parent_dir)
+    atoms_or_images = sb.convert_to_ase(index=index, **kwargs)
+    return atoms_or_images
+
+
+def read_aimd(filename, index=-1, **kwargs):
+    """Parse a .static file bundle using a wrapper around SparcBundle
+    The reader works only when other files (.ion, .inpt) exist.
+    """
+    parent_dir = Path(filename).parent
+    sb = SparcBundle(directory=parent_dir)
+    atoms_or_images = sb.convert_to_ase(index=index, **kwargs)
+    return atoms_or_images
+
+
 def register_ase_io_sparc(name="sparc"):
     """
     Monkey patching the ase.io and ase.io.formats
@@ -735,9 +778,7 @@ def _new_filetype(filename, read=True, guess=True):
     # Step 1: patch the ase.io.sparc module
     try:
         entry_points = next(
-            ep
-            for ep in pkg_resources.iter_entry_points("ase.io")
-            if ep.name == "sparc"
+            ep for ep in pkg_resources.iter_entry_points("ase.io") if ep.name == "sparc"
         )
         _monkey_mod = entry_points.load()
     except Exception as e:
@@ -791,6 +832,29 @@ def _new_filetype(filename, read=True, guess=True):
                     "Atomatic format inference for sparc is not correctly registered. "
                     "You may need to use format=sparc in ase.io.read and ase.io.write. "
                 )
+    # Add additional formats including .ion (r/w), .static, .geopt, .aimd
+    F(
+        "ion",
+        desc="SPARC .ion file",
+        module="sparc",
+        code="1S",
+        ext="ion",
+    )
+    F(
+        "static",
+        desc="SPARC single point results",
+        module="sparc",
+        code="+S",
+        ext="static",
+    )
+    F(
+        "geopt",
+        desc="SPARC geometric optimization results",
+        module="sparc",
+        code="+S",
+        ext="geopt",
+    )
+    F("aimd", desc="SPARC AIMD results", module="sparc", code="+S", ext="aimd")
 
     # TODO: remove print options as it may be redundant
-    print("Successfully registered sparc format with ase.io!")
+    print("Successfully registered sparc formats with ase.io!")
diff --git a/sparc/sparc_json_api/parameters b/sparc/sparc_json_api/parameters
new file mode 100644
index 00000000..7b12ffd8
--- /dev/null
+++ b/sparc/sparc_json_api/parameters
@@ -0,0 +1,2035 @@
+{
+ "sparc_version": "2023.09.05",
+ "categories": [
+  "system",
+  "scf",
+  "electrostatics",
+  "stress calculation",
+  "qmd",
+  "structural relaxation",
+  "print options",
+  "parallelization options"
+ ],
+ "parameters": {
+  "ACE_FLAG": {
+   "symbol": "ACE_FLAG",
+   "label": "ACE_FLAG",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "ACE_FLAG: 0",
+   "description": "Use ACE operator to accelarte the hybrid calculation.",
+   "remark": "Without ACE operator, the hybrid calculation will be way slower than with it on depending on the system size.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Use ACE operator to accelarte the hybrid calculation.",
+   "remark_raw": "Without ACE operator, the hybrid calculation will be way slower than with it on depending on the system size.",
+   "category": "scf"
+  },
+  "ATOM_TYPE": {
+   "symbol": "ATOM_TYPE",
+   "label": "ATOM_TYPE",
+   "type": "string",
+   "default": null,
+   "unit": "No unit",
+   "example": "ATOM_TYPE: Fe",
+   "description": "Atomic type symbol.",
+   "remark": "The atomic type symbol can be attached with a number, e.g., Fe1 and Fe2. This feature is useful if one needs to provide two different pseudopotential files (PSEUDO_POT) for the same element.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Atomic type symbol.",
+   "remark_raw": "The atomic type symbol can be attached with a number, e.g., Fe1 and Fe2. This feature is useful if one needs to provide two different pseudopotential files (\\hyperlink{PSEUDO_POT}{\\texttt{PSEUDO\\_POT}}) for the same element.",
+   "category": "system"
+  },
+  "BC": {
+   "symbol": "BC",
+   "label": "BC",
+   "type": "string",
+   "default": null,
+   "unit": "No unit",
+   "example": "BC: P D D",
+   "description": "A set of three whitespace delimited characters specifying the boundary conditions in the lattice vector directions, respectively. P represents periodic boundary conditions and D represents Dirichlet boundary conditions.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "A set of three whitespace delimited characters specifying the boundary conditions in the lattice vector directions, respectively. \\texttt{P} represents periodic boundary conditions and \\texttt{D} represents Dirichlet boundary conditions.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "system"
+  },
+  "CALC_PRES": {
+   "symbol": "CALC_PRES",
+   "label": "CALC_PRES",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "CALC_PRES: 1",
+   "description": "Flag for calculation of the pressure.",
+   "remark": "Pressure is directly calculated, without calculation of the stress tensor.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for calculation of the pressure.",
+   "remark_raw": "Pressure is directly calculated, without calculation of the stress tensor.",
+   "category": "stress calculation"
+  },
+  "CALC_STRESS": {
+   "symbol": "CALC_STRESS",
+   "label": "CALC_STRESS",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "CALC_STRESS: 1",
+   "description": "Flag for calculation of the Hellmann-Feynman stress tensor (in cartesian coordinates).",
+   "remark": "",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for calculation of the Hellmann-Feynman stress tensor (in cartesian coordinates).",
+   "remark_raw": "%",
+   "category": "stress calculation"
+  },
+  "CELL": {
+   "symbol": "CELL",
+   "label": "CELL",
+   "type": "double array",
+   "default": null,
+   "unit": "Bohr",
+   "example": "CELL: 10.20 11.21 7.58",
+   "description": "A set of three whitespace delimited values specifying the cell lengths in the lattice vector (LATVEC) directions, respectively.",
+   "remark": "Note that CELL ignores the lengths of the lattice vectors specified in the .inpt file and only treats them as unit vectors. LATVEC_SCALE and CELL cannot be specified simultaneously.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "A set of three whitespace delimited values specifying the cell lengths in the lattice vector (\\hyperlink{LATVEC}{\\texttt{LATVEC}}) directions, respectively.",
+   "remark_raw": "Note that \\hyperlink{CELL}{\\texttt{CELL}} ignores the lengths of the lattice vectors specified in the \\texttt{.inpt} file and only treats them as unit vectors. \\hyperlink{LATVEC_SCALE}{\\texttt{LATVEC\\_SCALE}} and \\hyperlink{CELL}{\\texttt{CELL}} cannot be specified simultaneously.",
+   "category": "system"
+  },
+  "CHEB_DEGREE": {
+   "symbol": "CHEB_DEGREE",
+   "label": "CHEB_DEGREE",
+   "type": "integer",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "CHEB_DEGEE: 25",
+   "description": "Degree of polynomial used for Chebyshev filtering.",
+   "remark": "For larger mesh-sizes, smaller values of CHEB_DEGREE are generally more efficient, and vice-versa.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically set.",
+   "description_raw": "Degree of polynomial used for Chebyshev filtering.",
+   "remark_raw": "For larger mesh-sizes, smaller values of \\texttt{CHEB\\_DEGREE} are generally more efficient, and vice-versa.",
+   "category": "scf"
+  },
+  "CHEFSI_BOUND_FLAG": {
+   "symbol": "CHEFSI_BOUND_FLAG",
+   "label": "CHEFSI_BOUND_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "CHEFSI_BOUND_FLAG: 1",
+   "description": "Flag to recalculate the bounds for Chebyshev filtering. If set to 0, then only for the very first SCF will the upper bound be evaluated based on the maximum eigenvalue using Lanczos algorithm, and the upper bound will be the same for the rest steps. If set to 1, the upper bound will be reevaluated for every SCF.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag to recalculate the bounds for Chebyshev filtering. If set to $0$, then only for the very first SCF will the upper bound be evaluated based on the maximum eigenvalue using Lanczos algorithm, and the upper bound will be the same for the rest steps. If set to $1$, the upper bound will be reevaluated for every SCF.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "COORD": {
+   "symbol": "COORD",
+   "label": "COORD",
+   "type": "double array",
+   "default": null,
+   "unit": "Bohr",
+   "example": "COORD: 0.0 0.0 0.0 \n2.5 2.5 2.5",
+   "description": "The Cartesian coordinates of atoms of a ATOM_TYPE specified before this variable. If the coordinates are outside the fundamental domain (see CELL and LATVEC) in the periodic directions (see BC), it will be automatically mapped back to the domain.",
+   "remark": "For a system with different types of atoms, one has to specify the coordinates for every ATOM_TYPE. One can also specify the coordinates of the atoms using COORD_FRAC.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "The Cartesian coordinates of atoms of a \\hyperlink{ATOM_TYPE}{\\texttt{ATOM\\_TYPE}} specified before this variable. If the coordinates are outside the fundamental domain (see \\hyperlink{CELL}{\\texttt{CELL}} and \\hyperlink{LATVEC}{\\texttt{LATVEC}}) in the periodic directions (see \\hyperlink{BC}{\\texttt{BC}}), it will be automatically mapped back to the domain.",
+   "remark_raw": "For a system with different types of atoms, one has to specify the coordinates for every \\hyperlink{ATOM_TYPE}{\\texttt{ATOM\\_TYPE}}. One can also specify the coordinates of the atoms using \\hyperlink{COORD_FRAC}{\\texttt{COORD\\_FRAC}}.",
+   "category": "system"
+  },
+  "COORD_FRAC": {
+   "symbol": "COORD_FRAC",
+   "label": "COORD_FRAC",
+   "type": "double array",
+   "default": null,
+   "unit": "None",
+   "example": "COORD_FRAC: 0.5 0.5 0.0 \n0.0 0.5 0.5",
+   "description": "The fractional coordinates of atoms of a ATOM_TYPE specified before this variable. COORD_FRAC(i,j) x CELL(j), (j = 1,2,3) gives the coordinate of the i^th atom along the j^th LATVEC direction. If the coordinates are outside the fundamental domain (see CELL and LATVEC) in the periodic directions (see BC), it will be automatically mapped back to the domain.",
+   "remark": "For a system with different types of atoms, one has to specify the coordinates for every ATOM_TYPE. One can also specify the coordinates of the atoms using COORD.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "The fractional coordinates of atoms of a \\hyperlink{ATOM_TYPE}{\\texttt{ATOM\\_TYPE}} specified before this variable. \\texttt{COORD\\_FRAC}$(i,j)$ $\\times$ \\hyperlink{CELL}{\\texttt{CELL}}$(j)$, $(j = 1,2,3)$ gives the coordinate of the $i^{th}$ atom along the $j^{th}$ \\hyperlink{LATVEC}{\\texttt{LATVEC}} direction. If the coordinates are outside the fundamental domain (see \\hyperlink{CELL}{\\texttt{CELL}} and \\hyperlink{LATVEC}{\\texttt{LATVEC}}) in the periodic directions (see \\hyperlink{BC}{\\texttt{BC}}), it will be automatically mapped back to the domain.",
+   "remark_raw": "For a system with different types of atoms, one has to specify the coordinates for every \\hyperlink{ATOM_TYPE}{\\texttt{ATOM\\_TYPE}}. One can also specify the coordinates of the atoms using \\hyperlink{COORD}{\\texttt{COORD}}.",
+   "category": "system"
+  },
+  "D3_CN_THR": {
+   "symbol": "D3_CN_THR",
+   "label": "D3_CN_THR",
+   "type": "double",
+   "default": 625.0,
+   "unit": "Bohr$^2$",
+   "example": "D3_CN_THR: 1600",
+   "description": "Square of cut-off radius for calculating CN value of every atom and DFT-D3 correction between three atoms",
+   "remark": "Only applicable when DFT-D3 correction D3_FLAG is used.\nD3_CN_THR should be smaller or equal to cutoff radius of DFT-D3 correction cutoff radius, D3_RTHR.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{625}",
+   "description_raw": "Square of cut-off radius for calculating CN value of every atom and DFT-D3 correction between three atoms",
+   "remark_raw": "Only applicable when DFT-D3 correction \\hyperlink{D3_FLAG}{\\texttt{D3\\_FLAG}} is used.\n\nD3\\_CN\\_THR should be smaller or equal to cutoff radius of DFT-D3 correction cutoff radius, \\hyperlink{D3_RTHR}{\\texttt{D3\\_RTHR}}.",
+   "category": "system"
+  },
+  "D3_FLAG": {
+   "symbol": "D3_FLAG",
+   "label": "D3_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "D3_FLAG: 1",
+   "description": "Flag for adding Grimme\u2019s DFT-D3 correction on the result",
+   "remark": "Only active when using GGA-PBE, GGA-RPBE and GGA-PBEsol.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{0}",
+   "description_raw": "Flag for adding Grimme\u2019s DFT-D3 correction on the result",
+   "remark_raw": "Only active when using GGA-PBE, GGA-RPBE and GGA-PBEsol.",
+   "category": "system"
+  },
+  "D3_RTHR": {
+   "symbol": "D3_RTHR",
+   "label": "D3_RTHR",
+   "type": "double",
+   "default": 1600.0,
+   "unit": "Bohr$^2$",
+   "example": "D3_RTHR: 9000",
+   "description": "Square of cut-off radius for calculating DFT-D3 correction between two atoms",
+   "remark": "Only applicable when DFT-D3 correction D3_FLAG is used.\nD3_RTHR should be larger or equal to cutoff radius of CN coefficient, D3_CN_THR.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{1600}",
+   "description_raw": "Square of cut-off radius for calculating DFT-D3 correction between two atoms",
+   "remark_raw": "Only applicable when DFT-D3 correction \\hyperlink{D3_FLAG}{\\texttt{D3\\_FLAG}} is used.\n\nD3\\_RTHR should be larger or equal to cutoff radius of CN coefficient, \\hyperlink{D3_CN_THR}{\\texttt{D3\\_CN\\_THR}}.",
+   "category": "system"
+  },
+  "ECUT": {
+   "symbol": "ECUT",
+   "label": "ECUT",
+   "type": "double",
+   "default": null,
+   "unit": "Ha",
+   "example": "ECUT: 30",
+   "description": "Equivalent plane-wave energy cutoff, based on which MESH_SPACING will be automatically calculated.",
+   "remark": "This is not exact, but rather an estimate. ECUT, MESH_SPACING, FD_GRID cannot be specified simultaneously.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Equivalent plane-wave energy cutoff, based on which \\hyperlink{MESH_SPACING}{\\texttt{MESH\\_SPACING}} will be automatically calculated.",
+   "remark_raw": "This is not exact, but rather an estimate. \\hyperlink{ECUT}{\\texttt{ECUT}}, \\hyperlink{MESH_SPACING}{\\texttt{MESH\\_SPACING}}, \\hyperlink{FD_GRID}{\\texttt{FD\\_GRID}} cannot be specified simultaneously.",
+   "category": "system"
+  },
+  "EIG_PARAL_BLKSZ": {
+   "symbol": "EIG_PARAL_BLKSZ",
+   "label": "EIG_PARAL_BLKSZ",
+   "type": "integer",
+   "default": 128,
+   "unit": "No unit",
+   "example": "EIG_PARAL_BLKSZ: 64",
+   "description": "Block size for the distribution of matrix in block-cyclic format in a parallel algorithm for solving the subspace eigenproblem.",
+   "allow_bool_input": false,
+   "default_remark": "128",
+   "description_raw": "Block size for the distribution of matrix in block-cyclic format in a parallel algorithm for solving the subspace eigenproblem.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "parallelization options"
+  },
+  "EIG_PARAL_MAXNP": {
+   "symbol": "EIG_PARAL_MAXNP",
+   "label": "EIG_PARAL_MAXNP",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "EIG_PARAL_MAXNP: 36",
+   "description": "Maximum number of processors used in parallel eigensolver. The number is machine dependent. Users could provide their own value for best performance.",
+   "allow_bool_input": false,
+   "default_remark": "Generated by a linear model.",
+   "description_raw": "Maximum number of processors used in parallel eigensolver. The number is machine dependent. Users could provide their own value for best performance.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "parallelization options"
+  },
+  "EIG_PARAL_ORFAC": {
+   "symbol": "EIG_PARAL_ORFAC",
+   "label": "EIG_PARAL_ORFAC",
+   "type": "double",
+   "default": 0.0,
+   "unit": "No unit",
+   "example": "EIG_PARAL_ORFAC: 0.001",
+   "description": "Specifies which eigenvectors should be reorthogonalized when using the parallel eigensolver \np?syevx \nor \np?sygvx \nfor solving the subspace eigenproblem. \nThe parallel eigensolvers can be turned on using the EIG_SERIAL_MAXNS flag.\nNo reorthogonalization will be done if EIG_PARAL_ORFAC equals zero. A default value of 0.001 is used if EIG_PARAL_ORFAC is negative. \nNote that reorthogonalization of eigenvectors is extremely time-consuming.",
+   "allow_bool_input": false,
+   "default_remark": "0.0",
+   "description_raw": "Specifies which eigenvectors should be reorthogonalized when using the parallel eigensolver \n\\href{https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/scalapack-routines/scalapack-driver-routines/p-syevx.html}{\\texttt{p?syevx}} \nor \n\\href{https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/scalapack-routines/scalapack-driver-routines/p-sygvx.html}{\\texttt{p?sygvx}} \nfor solving the subspace eigenproblem. \nThe parallel eigensolvers can be turned on using the \\hyperlink{EIG_SERIAL_MAXNS}{\\texttt{EIG\\_SERIAL\\_MAXNS}} flag.\nNo reorthogonalization will be done if \\texttt{EIG\\_PARAL\\_ORFAC} equals zero. A default value of $0.001$ is used if \\texttt{EIG\\_PARAL\\_ORFAC} is negative. \nNote that reorthogonalization of eigenvectors is extremely time-consuming.\n% See the corresponding documentation for more details.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "parallelization options"
+  },
+  "EIG_SERIAL_MAXNS": {
+   "symbol": "EIG_SERIAL_MAXNS",
+   "label": "EIG_SERIAL_MAXNS",
+   "type": "integer",
+   "default": 2000,
+   "unit": "No unit",
+   "example": "EIG_SERIAL_MAXNS: 1000",
+   "description": "Maximum NSTATES value up to which a serial algorithm will be used to solve the subspace eigenproblem.",
+   "remark": "If one wants to use a parallel algorithm to solve the subspace eigenproblem for all cases, simply set EIG_SERIAL_MAXNS to 0. Alternatively, set EIG_SERIAL_MAXNS to a very large value to always use serial algorithm.",
+   "allow_bool_input": false,
+   "default_remark": "2000",
+   "description_raw": "Maximum \\hyperlink{NSTATES}{\\texttt{NSTATES}} value up to which a serial algorithm will be used to solve the subspace eigenproblem.",
+   "remark_raw": "If one wants to use a parallel algorithm to solve the subspace eigenproblem for all cases, simply set \\texttt{EIG\\_SERIAL\\_MAXNS} to $0$. Alternatively, set \\texttt{EIG\\_SERIAL\\_MAXNS} to a very large value to always use serial algorithm.",
+   "category": "parallelization options"
+  },
+  "ELEC_TEMP": {
+   "symbol": "ELEC_TEMP",
+   "label": "ELEC_TEMP",
+   "type": "double",
+   "default": null,
+   "unit": "Kelvin",
+   "example": "ELEC_TEMP: 315.775",
+   "description": "Electronic temperature.",
+   "remark": "This is equivalent to setting SMEARING (0.001 Ha = 315.775 Kelvin).",
+   "allow_bool_input": false,
+   "default_remark": "2320.904 for \\texttt{gaussian} \\\\\n1160.452 for \\texttt{fermi-dirac}",
+   "description_raw": "Electronic temperature.",
+   "remark_raw": "This is equivalent to setting \\hyperlink{SMEARING}{\\texttt{SMEARING}} (0.001 Ha = 315.775 Kelvin).",
+   "category": "system"
+  },
+  "ELEC_TEMP_TYPE": {
+   "symbol": "ELEC_TEMP_TYPE",
+   "label": "ELEC_TEMP_TYPE",
+   "type": "string",
+   "default": "gaussian",
+   "unit": "No unit",
+   "example": "ELEC_TEMP_TYPE: fd",
+   "description": "Function used for the smearing (electronic temperature). Options are: fermi-dirac (or fd), gaussian.",
+   "remark": "Use ELEC_TEMP or SMEARING to set smearing value.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{gaussian}",
+   "description_raw": "Function used for the smearing (electronic temperature). Options are: \\texttt{fermi-dirac} (or \\texttt{fd}), \\texttt{gaussian}.",
+   "remark_raw": "Use \\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}} or \\hyperlink{SMEARING}{\\texttt{SMEARING}} to set smearing value.",
+   "category": "system"
+  },
+  "EXCHANGE_CORRELATION": {
+   "symbol": "EXCHANGE_CORRELATION",
+   "label": "EXCHANGE_CORRELATION",
+   "type": "string",
+   "default": null,
+   "unit": "No unit",
+   "example": "EXCHANGE_CORRELATION: LDA_PW",
+   "description": "Choice of exchange-correlation functional. Options are LDA_PW (Perdew-Wang LDA), LDA_PZ (Purdew-Zunger LDA), \nGGA_PBE (PBE GGA), GGA_RPBE (revised PBE GGA), and GGA_PBEsol (PBE GGA revised for solids), \nPBE0, HF (Hartree-Fock), HSE,\nvdWDF1 (van der Waals Density Functional developed by Dion et al.), vdWDF2 (vdW Density Functional modified by Lee et al), SCAN (SCAN metaGGA), RSCAN (rSCAN metaGGA), and R2SCAN (r2SCAN metaGGA).",
+   "remark": "For spin-polarized calculation (SPIN_TYP = 1), LDA_PZ is not available.\nCurrently SCAN, RSCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
+   "allow_bool_input": false,
+   "default_remark": "No Default",
+   "description_raw": "Choice of exchange-correlation functional. Options are \\texttt{LDA\\_PW} (Perdew-Wang LDA), \\texttt{LDA\\_PZ} (Purdew-Zunger LDA), \n\\texttt{GGA\\_PBE} (PBE GGA), \\texttt{GGA\\_RPBE} (revised PBE GGA), and \\texttt{GGA\\_PBEsol} (PBE GGA revised for solids), \n\\texttt{PBE0}, \\texttt{HF} (Hartree-Fock), \\texttt{HSE},\n\\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \\texttt{vdWDF2} (vdW Density Functional modified by Lee et al), \\texttt{SCAN} (SCAN metaGGA), \\texttt{RSCAN} (rSCAN metaGGA), and \\texttt{R2SCAN} (r2SCAN metaGGA).",
+   "remark_raw": "For spin-polarized calculation (\\hyperlink{SPIN_TYP}{\\texttt{SPIN\\_TYP}} = 1), \\texttt{LDA\\_PZ} is not available.\n\nCurrently SCAN, RSCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
+   "category": "system"
+  },
+  "EXX_ACE_VALENCE_STATES": {
+   "symbol": "EXX_ACE_VALENCE_STATES",
+   "label": "EXX_ACE_VALENCE_STATES",
+   "type": "integer",
+   "default": 3,
+   "unit": "No unit",
+   "example": "EXX_ACE_VALENCE_STATES: 1",
+   "description": "Control of number of unoccupied states used to construct ACE operator.",
+   "remark": "Only active when using hybrid functionals with ACE operator.",
+   "allow_bool_input": false,
+   "default_remark": "3",
+   "description_raw": "Control of number of unoccupied states used to construct ACE operator.",
+   "remark_raw": "Only active when using hybrid functionals with ACE operator.",
+   "category": "scf"
+  },
+  "EXX_DIVERGENCE": {
+   "symbol": "EXX_DIVERGENCE",
+   "label": "EXX_DIVERGENCE",
+   "type": "string",
+   "default": "SPHERICAL",
+   "unit": "No unit",
+   "example": "EXX_DIVERGENCE: AUXILIARY",
+   "description": "Treatment of divergence in exact exchange. Options are SPHERICAL (spherical truncation), \nAUXILIARY (auxiliary function method) and ERFC (erfc screening).",
+   "remark": "For systems with cube-like geometry, both methods converge fast. For slab and wire, auxiliary function method is a better option. \nERFC screening is the default option for HSE in bulk and molecule simulation.",
+   "allow_bool_input": false,
+   "default_remark": "SPHERICAL",
+   "description_raw": "Treatment of divergence in exact exchange. Options are \\texttt{SPHERICAL} (spherical truncation), \n\\texttt{AUXILIARY} (auxiliary function method) and \\texttt{ERFC} (erfc screening).",
+   "remark_raw": "For systems with cube-like geometry, both methods converge fast. For slab and wire, auxiliary function method is a better option. \nERFC screening is the default option for HSE in bulk and molecule simulation.",
+   "category": "scf"
+  },
+  "EXX_DOWNSAMPLING": {
+   "symbol": "EXX_DOWNSAMPLING",
+   "label": "EXX_DOWNSAMPLING",
+   "type": "integer array",
+   "default": [
+    1,
+    1,
+    1
+   ],
+   "unit": "No unit",
+   "example": "EXX_DOWNSAMPLING: 1 2 3",
+   "description": "Down-sampling of k-points grids. There should be 3 nonnegative integers. 0 means using 0 k-point in that direction, \nrequiring 0 is one of the k-point after time-reversal symmetry in that direction. \nPositive value should be a factor of the number of grid points in that direction.",
+   "allow_bool_input": false,
+   "default_remark": "1 1 1",
+   "description_raw": "Down-sampling of k-points grids. There should be 3 nonnegative integers. 0 means using 0 k-point in that direction, \nrequiring 0 is one of the k-point after time-reversal symmetry in that direction. \nPositive value should be a factor of the number of grid points in that direction.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "EXX_FRAC": {
+   "symbol": "EXX_FRAC",
+   "label": "EXX_FRAC",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "EXX_FRAC: 0.3",
+   "description": "Fraction of exact exchange in hybrid functional, e.g. PBE0 and HSE, while the fraction of PBE is 1-EXX_FRAC",
+   "allow_bool_input": false,
+   "default_remark": "0.25 for PBE0 and HSE",
+   "description_raw": "Fraction of exact exchange in hybrid functional, e.g. PBE0 and HSE, while the fraction of PBE is 1-\\texttt{EXX\\_FRAC}",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "EXX_MEM": {
+   "symbol": "EXX_MEM",
+   "label": "EXX_MEM",
+   "type": "integer",
+   "default": 20,
+   "unit": "No unit",
+   "example": "EXX_MEM: 0",
+   "description": "Number of Poisson's equations to be solved in each process at a time when creating exact exchange operator or ACE operator. Typically, when EXX_MEM is larger than 20, the speed of code is barely affected. When it is 0, all Poisson's equations are solved together and it hits the fastest speed but largest memory requirement.",
+   "allow_bool_input": false,
+   "default_remark": "20",
+   "description_raw": "Number of Poisson's equations to be solved in each process at a time when creating exact exchange operator or ACE operator. Typically, when \\texttt{EXX\\_MEM} is larger than 20, the speed of code is barely affected. When it is 0, all Poisson's equations are solved together and it hits the fastest speed but largest memory requirement.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "EXX_METHOD": {
+   "symbol": "EXX_METHOD",
+   "label": "EXX_METHOD",
+   "type": "string",
+   "default": "FOURIER_SPACE",
+   "unit": "No unit",
+   "example": "EXX_METHOD: REAL_SPACE",
+   "description": "Methods to solve Poisson's equation in Exact Exchange part. Options include using FFT to solve it in Fourier space and using linear solver, like CG, to solve in Real space.",
+   "remark": "Only active when using hybrid functionals for molecule simulation, like PBE0 and HSE. \nFOURIER_SPACE method is much faster than REAL_SPACE method.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{FOURIER\\_SPACE}",
+   "description_raw": "Methods to solve Poisson's equation in Exact Exchange part. Options include using FFT to solve it in Fourier space and using linear solver, like CG, to solve in Real space.",
+   "remark_raw": "Only active when using hybrid functionals for molecule simulation, like PBE0 and HSE. \n\\texttt{FOURIER\\_SPACE} method is much faster than \\texttt{REAL\\_SPACE} method.",
+   "category": "scf"
+  },
+  "EXX_RANGE_FOCK": {
+   "symbol": "EXX_RANGE_FOCK",
+   "label": "EXX_RANGE_FOCK",
+   "type": "double",
+   "default": 0.1587,
+   "unit": "1/Bohr",
+   "example": "EXX_RANGE_FOCK: 0.106",
+   "description": "Short range screen parameter of hartree-fock operator in HSE functional.",
+   "remark": "Default is using VASP's HSE03 value. Different code has different parameters. Be careful with the results.",
+   "allow_bool_input": false,
+   "default_remark": "0.1587",
+   "description_raw": "Short range screen parameter of hartree-fock operator in HSE functional.",
+   "remark_raw": "Default is using VASP's HSE03 value. Different code has different parameters. Be careful with the results.",
+   "category": "system"
+  },
+  "EXX_RANGE_PBE": {
+   "symbol": "EXX_RANGE_PBE",
+   "label": "EXX_RANGE_PBE",
+   "type": "double",
+   "default": 0.1587,
+   "unit": "1/Bohr",
+   "example": "EXX_RANGE_PBE: 0.106",
+   "description": "Short range screen parameter of PBE in HSE functional.",
+   "remark": "Default is using VASP's HSE03 value. Different code has different parameters. Be careful with the results.",
+   "allow_bool_input": false,
+   "default_remark": "0.1587",
+   "description_raw": "Short range screen parameter of PBE in HSE functional.",
+   "remark_raw": "Default is using VASP's HSE03 value. Different code has different parameters. Be careful with the results.",
+   "category": "system"
+  },
+  "FD_GRID": {
+   "symbol": "FD_GRID",
+   "label": "FD_GRID",
+   "type": "integer array",
+   "default": null,
+   "unit": "No unit",
+   "example": "FD_GRID: 26 26 30",
+   "description": "A set of three whitespace delimited values specifying the number of finite-difference intervals in the lattice vector (LATVEC) directions, respectively.",
+   "remark": "The convergence of results with respect to spatial discretization needs to be verified. ECUT, MESH_SPACING, FD_GRID cannot be specified simultaneously.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "A set of three whitespace delimited values specifying the number of finite-difference intervals in the lattice vector (\\hyperlink{LATVEC}{\\texttt{LATVEC}}) directions, respectively.",
+   "remark_raw": "The convergence of results with respect to spatial discretization needs to be verified. \\hyperlink{ECUT}{\\texttt{ECUT}}, \\hyperlink{MESH_SPACING}{\\texttt{MESH\\_SPACING}}, \\hyperlink{FD_GRID}{\\texttt{FD\\_GRID}} cannot be specified simultaneously.",
+   "category": "system"
+  },
+  "FD_ORDER": {
+   "symbol": "FD_ORDER",
+   "label": "FD_ORDER",
+   "type": "integer",
+   "default": 12,
+   "unit": "No unit",
+   "example": "FD_ORDER: 12",
+   "description": "Order of the finite-difference approximation.",
+   "remark": "Restricted to even integers since central finite-differences are employed. The default value of 12 has been found to be an efficient choice for most systems.",
+   "allow_bool_input": false,
+   "default_remark": "12",
+   "description_raw": "Order of the finite-difference approximation.",
+   "remark_raw": "Restricted to even integers since central finite-differences are employed. The default value of 12 has been found to be an efficient choice for most systems.",
+   "category": "system"
+  },
+  "FIRE_DT": {
+   "symbol": "FIRE_DT",
+   "label": "FIRE_DT",
+   "type": "double",
+   "default": 1.0,
+   "unit": "Femto second",
+   "example": "FIRE_DT: 0.1",
+   "description": "Time step used in FIRE (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "Time step used in FIRE (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "FIRE_MASS": {
+   "symbol": "FIRE_MASS",
+   "label": "FIRE_MASS",
+   "type": "double",
+   "default": 1.0,
+   "unit": "Atomic mass unit",
+   "example": "FIRE_MASS: 2.5",
+   "description": "Pseudomass used in FIRE (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "1.0",
+   "description_raw": "Pseudomass used in FIRE (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "FIRE_MAXMOV": {
+   "symbol": "FIRE_MAXMOV",
+   "label": "FIRE_MAXMOV",
+   "type": "double",
+   "default": 0.2,
+   "unit": "Bohr",
+   "example": "FIRE_MAXMOV: 1.0",
+   "description": "Maximum movement for any atom in FIRE (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "0.2",
+   "description_raw": "Maximum movement for any atom in FIRE (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "FIX_RAND": {
+   "symbol": "FIX_RAND",
+   "label": "FIX_RAND",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "FIX_RAND: 1",
+   "description": "Flag to fix the random seeds for setting initial guesses. Once set to 1, the random seeds will be fixed for different runs and for different numbers of processors. This option will make sure the answers will be exactly the same (up to machine precision) when SPARC is executed with different numbers of processors.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag to fix the random seeds for setting initial guesses. Once set to $1$, the random seeds will be fixed for different runs and for different numbers of processors. This option will make sure the answers will be exactly the same (up to machine precision) when SPARC is executed with different numbers of processors.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "ION_TEMP": {
+   "symbol": "ION_TEMP",
+   "label": "ION_TEMP",
+   "type": "double",
+   "default": null,
+   "unit": "Kelvin",
+   "example": "ION_TEMP: 315",
+   "description": "Starting ionic temperature in QMD, used to generate initial velocity distribution.",
+   "remark": "Must be specified if MD_FLAG is set to 1. It is also the target temperature in MD_METHOD NPT_NH and NPT_NP.",
+   "allow_bool_input": false,
+   "default_remark": "No Default",
+   "description_raw": "Starting ionic temperature in QMD, used to generate initial velocity distribution.",
+   "remark_raw": "Must be specified if \\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} is set to $1$. It is also the target temperature in \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} NPT\\_NH and NPT\\_NP.",
+   "category": "qmd"
+  },
+  "ION_TEMP_END": {
+   "symbol": "ION_TEMP_END",
+   "label": "ION_TEMP_END",
+   "type": "double",
+   "default": null,
+   "unit": "Kelvin",
+   "example": "ION_TEMP_END: 100",
+   "description": "Specifies the final temperature of the thermostat. Thermostat temperature is varied linearly from ION_TEMP to ION_TEMP_END with respect to time.",
+   "remark": "Available for NVT_NH quantum molecular dynamics only. Not supported in NPT_NH and NPT_NP.",
+   "allow_bool_input": false,
+   "default_remark": "ION\\_TEMP",
+   "description_raw": "Specifies the final temperature of the thermostat. Thermostat temperature is varied linearly from \\hyperlink{ION_TEMP}{\\texttt{ION\\_TEMP}} to \\hyperlink{ION_TEMP_END}{\\texttt{ION\\_TEMP\\_END}} with respect to time.",
+   "remark_raw": "Available for NVT\\_NH quantum molecular dynamics only. Not supported in NPT\\_NH and NPT\\_NP.",
+   "category": "qmd"
+  },
+  "ION_VEL_DSTR": {
+   "symbol": "ION_VEL_DSTR",
+   "label": "ION_VEL_DSTR",
+   "type": "integer",
+   "default": 2,
+   "unit": "No unit",
+   "example": "ION_VEL_DSTR: 1",
+   "description": "Specifies the type of distribution for the initial velocity of atoms based on their initial temperature. Choose 1 for uniform velocity distribution and 2 for Maxwell-Boltzmann distribution.",
+   "remark": "Currently, the code supports only two options for the variable.",
+   "allow_bool_input": false,
+   "default_remark": "2",
+   "description_raw": "Specifies the type of distribution for the initial velocity of atoms based on their initial temperature. Choose $1$ for uniform velocity distribution and $2$ for Maxwell-Boltzmann distribution.",
+   "remark_raw": "Currently, the code supports only two options for the variable.",
+   "category": "qmd"
+  },
+  "ION_VEL_DSTR_RAND": {
+   "symbol": "ION_VEL_DSTR_RAND",
+   "label": "ION_VEL_DSTR_RAND",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "ION_VEL_DSTR_RAND: 1",
+   "description": "Flag to reseed the initial velocities of atoms in a QMD simulation. Set this flag to 1 to shuffle (change the random seed for) the initial velocities for different runs. Set this flag to 0 to maintain the same initial velocities.",
+   "remark": "This option is convenient for parallel statistics calculations.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag to reseed the initial velocities of atoms in a QMD simulation. Set this flag to $1$ to shuffle (change the random seed for) the initial velocities for different runs. Set this flag to $0$ to maintain the same initial velocities.",
+   "remark_raw": "This option is convenient for parallel statistics calculations.",
+   "category": "qmd"
+  },
+  "KPOINT_GRID": {
+   "symbol": "KPOINT_GRID",
+   "label": "KPOINT_GRID",
+   "type": "integer array",
+   "default": [
+    1,
+    1,
+    1
+   ],
+   "unit": "No unit",
+   "example": "KPOINT_GRID: 2 3 4",
+   "description": "Number of k-points in each direction of the Monkhorst-Pack grid for Brillouin zone integration.",
+   "remark": "Time-reversal symmetry is assumed to hold.",
+   "allow_bool_input": false,
+   "default_remark": "1 1 1",
+   "description_raw": "Number of k-points in each direction of the Monkhorst-Pack grid for Brillouin zone integration.",
+   "remark_raw": "Time-reversal symmetry is assumed to hold.",
+   "category": "system"
+  },
+  "KPOINT_SHIFT": {
+   "symbol": "KPOINT_SHIFT",
+   "label": "KPOINT_SHIFT",
+   "type": "double array",
+   "default": null,
+   "unit": "No unit",
+   "example": "KPOINT_SHIFT: 0.5 0.5 0.5",
+   "description": "Shift of k-points with respect to k-point grid containing \\Gamma-point.",
+   "remark": "The shift is in reduced coordinates.",
+   "allow_bool_input": false,
+   "default_remark": "0.0 for odd k-point mesh\\\\\n 0.5 for even k-point mesh",
+   "description_raw": "Shift of k-points with respect to k-point grid containing $\\Gamma$-point.",
+   "remark_raw": "The shift is in reduced coordinates.",
+   "category": "system"
+  },
+  "LATVEC": {
+   "symbol": "LATVEC",
+   "label": "LATVEC",
+   "type": "double array",
+   "default": [
+    [
+     1.0,
+     0.0,
+     0.0
+    ],
+    [
+     0.0,
+     1.0,
+     0.0
+    ],
+    [
+     0.0,
+     0.0,
+     1.0
+    ]
+   ],
+   "unit": "No unit",
+   "example": "LATVEC: 0.5 0.5 0.0\n0.0 0.5 0.5\n0.5 0.0 0.5",
+   "description": "A set of three vectors in row major order specifying the lattice vectors of the simulation domain (CELL).",
+   "remark": "",
+   "allow_bool_input": false,
+   "default_remark": "1.0 0.0 0.0\\\\\n0.0 1.0 0.0\\\\\n0.0 0.0 1.0\\\\",
+   "description_raw": "A set of three vectors in row major order specifying the lattice vectors of the simulation domain (\\hyperlink{CELL}{\\texttt{CELL}}).",
+   "remark_raw": "%Lattice vectors need not be unit vectors.\n%",
+   "category": "system"
+  },
+  "LATVEC_SCALE": {
+   "symbol": "LATVEC_SCALE",
+   "label": "LATVEC_SCALE",
+   "type": "double array",
+   "default": null,
+   "unit": "Bohr",
+   "example": "LATVEC_SCALE: 10.20 11.21 7.58",
+   "description": "A set of three whitespace delimited values specifying the scaling factors in the lattice vectors (LATVEC), respectively.",
+   "remark": "The difference between LATVEC_SCALE and CELL is that CELL treats the lattice vectors as unit vectors, whereas LATVEC_SCALE scales the lattice vectors directly as specified by the user. LATVEC_SCALE and CELL cannot be specified simultaneously.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "A set of three whitespace delimited values specifying the scaling factors in the lattice vectors (\\hyperlink{LATVEC}{\\texttt{LATVEC}}), respectively.",
+   "remark_raw": "The difference between \\hyperlink{LATVEC_SCALE}{\\texttt{LATVEC\\_SCALE}} and \\hyperlink{CELL}{\\texttt{CELL}} is that \\hyperlink{CELL}{\\texttt{CELL}} treats the lattice vectors as unit vectors, whereas \\hyperlink{LATVEC_SCALE}{\\texttt{LATVEC\\_SCALE}} scales the lattice vectors directly as specified by the user. \\hyperlink{LATVEC_SCALE}{\\texttt{LATVEC\\_SCALE}} and \\hyperlink{CELL}{\\texttt{CELL}} cannot be specified simultaneously.",
+   "category": "system"
+  },
+  "L_AUTOSCALE": {
+   "symbol": "L_AUTOSCALE",
+   "label": "L_AUTOSCALE",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "L_AUTOSCALE: 0",
+   "description": "Flag for automatically determining the inverse curvature that is used to determine the direction for next iteration in LBFGS (RELAX_METHOD).",
+   "remark": "Default works well in most cases.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for automatically determining the inverse curvature that is used to determine the direction for next iteration in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "L_FINIT_STP": {
+   "symbol": "L_FINIT_STP",
+   "label": "L_FINIT_STP",
+   "type": "double",
+   "default": 0.005,
+   "unit": "Bohr",
+   "example": "L_FINIT_STP: 0.01",
+   "description": "Step length for line optimizer in LBFGS (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "5e-3",
+   "description_raw": "Step length for line optimizer in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "L_HISTORY": {
+   "symbol": "L_HISTORY",
+   "label": "L_HISTORY",
+   "type": "integer",
+   "default": 20,
+   "unit": "No unit",
+   "example": "L_HISTORY: 15",
+   "description": "Size of history in LBFGS (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "20",
+   "description_raw": "Size of history in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "L_ICURV": {
+   "symbol": "L_ICURV",
+   "label": "L_ICURV",
+   "type": "double",
+   "default": 1.0,
+   "unit": "No unit",
+   "example": "L_ICURV: 0.1",
+   "description": "Initial inverse curvature, used to construct the inverse Hessian matrix in LBFGS (RELAX_METHOD).",
+   "remark": "Needed only if  L_AUTOSCALE is 0. Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "1.0",
+   "description_raw": "Initial inverse curvature, used to construct the inverse Hessian matrix in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Needed only if  \\hyperlink{L_AUTOSCALE}{\\texttt{L\\_AUTOSCALE}} is $0$. Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "L_LINEOPT": {
+   "symbol": "L_LINEOPT",
+   "label": "L_LINEOPT",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "L_LINEOPT: 0",
+   "description": "Flag for atomic force based line minimization in LBFGS (RELAX_METHOD).",
+   "remark": "Required only if L_AUTOSCALE is 0.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for atomic force based line minimization in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Required only if \\hyperlink{L_AUTOSCALE}{\\texttt{L\\_AUTOSCALE}} is $0$.",
+   "category": "structural relaxation"
+  },
+  "L_MAXMOV": {
+   "symbol": "L_MAXMOV",
+   "label": "L_MAXMOV",
+   "type": "double",
+   "default": 0.2,
+   "unit": "Bohr",
+   "example": "L_MAXMOV: 1.0",
+   "description": "The maximum allowed step size in LBFGS (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "0.2",
+   "description_raw": "The maximum allowed step size in LBFGS (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "MAXIT_FOCK": {
+   "symbol": "MAXIT_FOCK",
+   "label": "MAXIT_FOCK",
+   "type": "integer",
+   "default": 20,
+   "unit": "No unit",
+   "example": "MAXIT_FOCK: 50",
+   "description": "The maximum number of iterations for Hartree-Fock outer loop.",
+   "remark": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "allow_bool_input": false,
+   "default_remark": "20",
+   "description_raw": "The maximum number of iterations for Hartree-Fock outer loop.",
+   "remark_raw": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "category": "scf"
+  },
+  "MAXIT_POISSON": {
+   "symbol": "MAXIT_POISSON",
+   "label": "MAXIT_POISSON",
+   "type": "integer",
+   "default": 3000,
+   "unit": "No unit",
+   "example": "MAXIT_POISSON: 1000",
+   "description": "The maximum number of iterations for solving the Poisson equation using an iterative linear solver.",
+   "allow_bool_input": false,
+   "default_remark": "$3000$",
+   "description_raw": "The maximum number of iterations for solving the Poisson equation using an iterative linear solver.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "electrostatics"
+  },
+  "MAXIT_SCF": {
+   "symbol": "{MAXIT_SCF",
+   "label": "MAXIT_SCF",
+   "type": "integer",
+   "default": 100,
+   "unit": "No unit",
+   "example": "MAXIT_SCF: 50",
+   "description": "Maximum number of SCF iterations.",
+   "remark": "Larger values than the default of 100 may be required for highly inhomogeneous systems, particularly when small values of SMEARING/ELEC_TEMP are chosen.",
+   "allow_bool_input": false,
+   "default_remark": "100",
+   "description_raw": "Maximum number of SCF iterations.",
+   "remark_raw": "Larger values than the default of 100 may be required for highly inhomogeneous systems, particularly when small values of \\hyperlink{SMEARING}{\\texttt{SMEARING}}/\\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}} are chosen.",
+   "category": "scf"
+  },
+  "MD_FLAG": {
+   "symbol": "MD_FLAG",
+   "label": "MD_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "MD_FLAG: 1",
+   "description": "QMD simulations are performed if the flag is set to 1.",
+   "remark": "MD_FLAG and RELAX_FLAG both cannot be set to a value greater than 0.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "QMD simulations are performed if the flag is set to 1.",
+   "remark_raw": "\\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} and \\hyperlink{RELAX_FLAG}{\\texttt{RELAX\\_FLAG}} both cannot be set to a value greater than 0.",
+   "category": "qmd"
+  },
+  "MD_METHOD": {
+   "symbol": "MD_METHOD",
+   "label": "MD_METHOD",
+   "type": "string",
+   "default": "NVT_NH",
+   "unit": "No unit",
+   "example": "MD_METHOD: NVE",
+   "description": "Type of QMD to be performed. Currently, NVE (microcanonical ensemble), NVT_NH (canonical ensemble with Nose-Hoover thermostat), NVK_G (isokinetic ensemble with Gaussian thermostat), NPT_NH (isothermal-isobaric ensemble with Nose-Hoover thermostat) and NPT_NP (isothermal-isobaric ensemble with Nose-Poincare thermostat) are supported",
+   "allow_bool_input": false,
+   "default_remark": "NVT\\_NH",
+   "description_raw": "Type of QMD to be performed. Currently, NVE (microcanonical ensemble), NVT\\_NH (canonical ensemble with Nose-Hoover thermostat), NVK\\_G (isokinetic ensemble with Gaussian thermostat), NPT\\_NH (isothermal-isobaric ensemble with Nose-Hoover thermostat) and NPT\\_NP (isothermal-isobaric ensemble with Nose-Poincare thermostat) are supported",
+   "remark_raw": "",
+   "remark": "",
+   "category": "qmd"
+  },
+  "MD_NSTEP": {
+   "symbol": "MD_NSTEP",
+   "label": "MD_NSTEP",
+   "type": "integer",
+   "default": 10000000,
+   "unit": "No unit",
+   "example": "MD_NSTEP: 100",
+   "description": "Specifies the number of QMD steps.",
+   "remark": "If MD_NSTEP = N, the QMD runs from 0 to (N-1) x MD_TIMESTEP fs.",
+   "allow_bool_input": false,
+   "default_remark": "1e7",
+   "description_raw": "Specifies the number of QMD steps.",
+   "remark_raw": "If MD\\_NSTEP $= N$, the QMD runs from $0$ to $(N-1) \\times$ \\hyperlink{MD_TIMESTEP}{\\texttt{MD\\_TIMESTEP}} fs.",
+   "category": "qmd"
+  },
+  "MD_TIMESTEP": {
+   "symbol": "MD_TIMESTEP",
+   "label": "MD_TIMESTEP",
+   "type": "double",
+   "default": 1.0,
+   "unit": "Femtosecond",
+   "example": "MD_TIMESTEP: 0.1",
+   "description": "QMD time step.",
+   "remark": "Total QMD time is given by:  MD_TIMESTEP x (MD_NSTEP - 1).",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "QMD time step.",
+   "remark_raw": "Total QMD time is given by:  \\hyperlink{MD_TIMESTEP}{\\texttt{MD\\_TIMESTEP}} $\\times$ (\\hyperlink{MD_NSTEP}{\\texttt{MD\\_NSTEP}} - 1).",
+   "category": "qmd"
+  },
+  "MESH_SPACING": {
+   "symbol": "MESH_SPACING",
+   "label": "MESH_SPACING",
+   "type": "double",
+   "default": null,
+   "unit": "Bohr",
+   "example": "MESH_SPACING: 0.4",
+   "description": "Mesh spacing of the finite-difference grid.",
+   "remark": "The exact mesh-size will be determined by the size of the domain.  MESH_SPACING, FD_GRID, ECUT cannot be specified simultaneously.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Mesh spacing of the finite-difference grid.",
+   "remark_raw": "The exact mesh-size will be determined by the size of the domain.  \\hyperlink{MESH_SPACING}{\\texttt{MESH\\_SPACING}}, \\hyperlink{FD_GRID}{\\texttt{FD\\_GRID}}, \\hyperlink{ECUT}{\\texttt{ECUT}} cannot be specified simultaneously.",
+   "category": "system"
+  },
+  "MINIT_FOCK": {
+   "symbol": "MINIT_FOCK",
+   "label": "MINIT_FOCK",
+   "type": "integer",
+   "default": 2,
+   "unit": "No unit",
+   "example": "MINIT_FOCK: 3",
+   "description": "The minimum number of iterations for Hartree-Fock outer loop.",
+   "remark": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "allow_bool_input": false,
+   "default_remark": "2",
+   "description_raw": "The minimum number of iterations for Hartree-Fock outer loop.",
+   "remark_raw": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "category": "scf"
+  },
+  "MINIT_SCF": {
+   "symbol": "{MINIT_SCF",
+   "label": "MINIT_SCF",
+   "type": "integer",
+   "default": 3,
+   "unit": "No unit",
+   "example": "MAXIT_SCF: 5",
+   "description": "Minimum number of SCF iterations.",
+   "allow_bool_input": false,
+   "default_remark": "3",
+   "description_raw": "Minimum number of SCF iterations.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "MIXING_HISTORY": {
+   "symbol": "MIXING_HISTORY",
+   "label": "MIXING_HISTORY",
+   "type": "integer",
+   "default": 7,
+   "unit": "No unit",
+   "example": "MIXING_HISTORY: 40",
+   "description": "The mixing history used in Pulay mixing.",
+   "remark": "Too small values of MIXING_HISTORY can result in poor SCF convergence.",
+   "allow_bool_input": false,
+   "default_remark": "7",
+   "description_raw": "The mixing history used in Pulay mixing.",
+   "remark_raw": "Too small values of \\hyperlink{MIXING_HISTORY}{\\texttt{MIXING\\_HISTORY}} can result in poor SCF convergence.",
+   "category": "scf"
+  },
+  "MIXING_PARAMETER": {
+   "symbol": "MIXING_PARAMETER",
+   "label": "MIXING_PARAMETER",
+   "type": "double",
+   "default": 0.3,
+   "unit": "No unit",
+   "example": "MIXING_PARAMETER: 0.1",
+   "description": "The value of the relaxation parameter used in Pulay/simple mixing.",
+   "remark": "Values larger than the default value of 0.3 can be used for insulating systems, whereas smaller values are generally required for metallic  systems, particularly at small values of SMEARING or ELEC_TEMP.",
+   "allow_bool_input": false,
+   "default_remark": "0.3",
+   "description_raw": "The value of the relaxation parameter used in Pulay/simple mixing.",
+   "remark_raw": "Values larger than the default value of 0.3 can be used for insulating systems, whereas smaller values are generally required for metallic  systems, particularly at small values of \\hyperlink{SMEARING}{\\texttt{SMEARING}} or \\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}}.",
+   "category": "scf"
+  },
+  "MIXING_PARAMETER_MAG": {
+   "symbol": "MIXING_PARAMETER_MAG",
+   "label": "MIXING_PARAMETER_MAG",
+   "type": "double",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "MIXING_PARAMETER_MAG: 4.0",
+   "description": "The mixing parameter for the magnetization density in Pulay mixing for spin-polarized calculations.",
+   "remark": "For spin-polarized calculations, when SCF has difficulty to converge, increasing the mixing parameter to magnetization density might help. For example, setting it to 4.0, while turning off the preconditioner applied to the magnetization density (by setting MIXING_PRECOND_MAG to `none') is a good choice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically set to the same as \\hyperlink{MIXING_PARAMETER}{\\texttt{MIXING\\_PARAMETER}}.",
+   "description_raw": "The mixing parameter for the magnetization density in Pulay mixing for spin-polarized calculations.",
+   "remark_raw": "For spin-polarized calculations, when SCF has difficulty to converge, increasing the mixing parameter to magnetization density might help. For example, setting it to 4.0, while turning off the preconditioner applied to the magnetization density (by setting \\hyperlink{MIXING_PRECOND_MAG}{\\texttt{MIXING\\_PRECOND\\_MAG}} to `\\texttt{none}') is a good choice.",
+   "category": "scf"
+  },
+  "MIXING_PARAMETER_SIMPLE": {
+   "symbol": "MIXING_PARAMETER_SIMPLE",
+   "label": "MIXING_PARAMETER_SIMPLE",
+   "type": "double",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "MIXING_PARAMETER_SIMPLE: 0.1",
+   "description": "The value of the relaxation parameter used in the simple mixing step in the periodic Pulay scheme.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically set to the same as \\hyperlink{MIXING_PARAMETER}{\\texttt{MIXING\\_PARAMETER}}",
+   "description_raw": "The value of the relaxation parameter used in the simple mixing step in the periodic Pulay scheme.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "MIXING_PARAMETER_SIMPLE_MAG": {
+   "symbol": "MIXING_PARAMETER_SIMPLE_MAG",
+   "label": "MIXING_PARAMETER_SIMPLE_MAG",
+   "type": "double",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "MIXING_PARAMETER_SIMPLE_MAG: 4.0",
+   "description": "The value of the relaxation parameter for the magnetization density used in the simple mixing step in the periodic Pulay scheme for spin-polarized calculations.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically set to the same as \\hyperlink{MIXING_PARAMETER_MAG}{\\texttt{MIXING\\_PARAMETER\\_MAG}}",
+   "description_raw": "The value of the relaxation parameter for the magnetization density used in the simple mixing step in the periodic Pulay scheme for spin-polarized calculations.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "MIXING_PRECOND": {
+   "symbol": "MIXING_PRECOND",
+   "label": "MIXING_PRECOND",
+   "type": "string",
+   "default": "kerker",
+   "unit": "No unit",
+   "example": "MIXING_PRECOND: none",
+   "description": "This specifies the preconditioner used in the SCF iteration. Available options are: none, kerker.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{kerker}",
+   "description_raw": "This specifies the preconditioner used in the SCF iteration. Available options are: \\texttt{none}, \\texttt{kerker}.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "MIXING_PRECOND_MAG": {
+   "symbol": "MIXING_PRECOND_MAG",
+   "label": "MIXING_PRECOND_MAG",
+   "type": "string",
+   "default": null,
+   "unit": "No unit",
+   "example": "MIXING_PRECOND_MAG: kerker",
+   "description": "This specifies the preconditioner used for the magnetization density in the SCF iteration for spin-polarized calculations. Available options are: none, kerker.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{none}",
+   "description_raw": "This specifies the preconditioner used for the magnetization density in the SCF iteration for spin-polarized calculations. Available options are: \\texttt{none}, \\texttt{kerker}.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "MIXING_VARIABLE": {
+   "symbol": "MIXING_VARIABLE",
+   "label": "MIXING_VARIABLE",
+   "type": "string",
+   "default": "density",
+   "unit": "No unit",
+   "example": "MIXING_VARIABLE: potential",
+   "description": "This specifies whether potential or density mixing is performed in the SCF iteration. Available options are: potential and density.",
+   "allow_bool_input": false,
+   "default_remark": "\\texttt{density}",
+   "description_raw": "This specifies whether potential or density mixing is performed in the SCF iteration. Available options are: \\texttt{potential} and \\texttt{density}.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "NLCG_SIGMA": {
+   "symbol": "NLCG_SIGMA",
+   "label": "NLCG_SIGMA",
+   "type": "double",
+   "default": 0.5,
+   "unit": "No unit",
+   "example": "NLCG_SIGMA: 1",
+   "description": "Parameter in the secant method used to control the step length in NLCG (RELAX_METHOD).",
+   "remark": "Default value works well in most cases.",
+   "allow_bool_input": false,
+   "default_remark": "0.5",
+   "description_raw": "Parameter in the secant method used to control the step length in NLCG (\\hyperlink{RELAX_METHOD}{\\texttt{RELAX\\_METHOD}}).",
+   "remark_raw": "Default value works well in most cases.",
+   "category": "structural relaxation"
+  },
+  "NPT_NH_BMASS": {
+   "symbol": "NPT_NH_BMASS",
+   "label": "NPT_NH_BMASS",
+   "type": "double",
+   "default": null,
+   "unit": "atomic unit",
+   "example": "NPT_NH_BMASS: 5000",
+   "description": "Gives the inertia mass for the barostat variable in NPT_NH.",
+   "remark": "Applicable to NPT_NH MD_METHOD only.\nProgram will exit if NPT_NH is selected but NPT_NH_BMASS is not input",
+   "allow_bool_input": false,
+   "default_remark": "No default value",
+   "description_raw": "Gives the inertia mass for the barostat variable in NPT\\_NH.",
+   "remark_raw": "Applicable to NPT\\_NH \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.\nProgram will exit if NPT\\_NH is selected but NPT\\_NH\\_BMASS is not input",
+   "category": "qmd"
+  },
+  "NPT_NH_QMASS": {
+   "symbol": "NPT_NH_QMASS",
+   "label": "NPT_NH_QMASS",
+   "type": "other",
+   "default": null,
+   "unit": "atomic unit",
+   "example": "NPT_NH_QMASS: 2\n;700.0\n;700.0",
+   "description": "Gives the amount (first number) and inertia masses (others) of thermostats in NPT_NH. The maximum amount of thermostat variables of the Nose-Hoover chain is 60",
+   "remark": "Applicable to NPT_NH MD_METHOD only.\nProgram will exit if NPT_NH is selected but NPT_NH_QMASS is not input",
+   "allow_bool_input": false,
+   "default_remark": "No default value",
+   "description_raw": "Gives the amount (first number) and inertia masses (others) of thermostats in NPT\\_NH. The maximum amount of thermostat variables of the Nose-Hoover chain is 60",
+   "remark_raw": "Applicable to NPT\\_NH \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.\nProgram will exit if NPT\\_NH is selected but NPT\\_NH\\_QMASS is not input",
+   "category": "qmd"
+  },
+  "NPT_NP_BMASS": {
+   "symbol": "NPT_NP_BMASS",
+   "label": "NPT_NP_BMASS",
+   "type": "double",
+   "default": null,
+   "unit": "atomic unit",
+   "example": "NPT_NP_BMASS: 20",
+   "description": "Gives the inertia mass for the barostat variable in NPT_NP.",
+   "remark": "Applicable to NPT_NP MD_METHOD only.\nProgram will exit if NPT_NP is selected but NPT_NP_BMASS is not input",
+   "allow_bool_input": false,
+   "default_remark": "No default value",
+   "description_raw": "Gives the inertia mass for the barostat variable in NPT\\_NP.",
+   "remark_raw": "Applicable to NPT\\_NP \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.\nProgram will exit if NPT\\_NP is selected but NPT\\_NP\\_BMASS is not input",
+   "category": "qmd"
+  },
+  "NPT_NP_QMASS": {
+   "symbol": "NPT_NP_QMASS",
+   "label": "NPT_NP_QMASS",
+   "type": "double",
+   "default": null,
+   "unit": "atomic unit",
+   "example": "NPT_NP_QMASS: 100",
+   "description": "Gives the inertia mass for the thermostat variable in NPT_NP.",
+   "remark": "Applicable to NPT_NP MD_METHOD only.\nProgram will exit if NPT_NP is selected but NPT_NP_BMASS is not input",
+   "allow_bool_input": false,
+   "default_remark": "No default value",
+   "description_raw": "Gives the inertia mass for the thermostat variable in NPT\\_NP.",
+   "remark_raw": "Applicable to NPT\\_NP \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.\nProgram will exit if NPT\\_NP is selected but NPT\\_NP\\_BMASS is not input",
+   "category": "qmd"
+  },
+  "NPT_SCALE_VECS": {
+   "symbol": "NPT_SCALE_VECS",
+   "label": "NPT_SCALE_VECS",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "NPT_SCALE_VECS: 2",
+   "description": "Specify which lattice vectors can be rescaled in NPT-NH",
+   "remark": "Only numbers 1, 2 and 3 can be accepted. For example, if `` 2 3\" is the input, the cell will only expand or shrink in the directions of lattice vector 2 and lattice vector 3. Only NPT-NH can specify the rescaled vector",
+   "allow_bool_input": false,
+   "default_remark": "1 2 3",
+   "description_raw": "Specify which lattice vectors can be rescaled in NPT-NH",
+   "remark_raw": "Only numbers 1, 2 and 3 can be accepted. For example, if `` 2 3\" is the input, the cell will only expand or shrink in the directions of lattice vector 2 and lattice vector 3. Only NPT-NH can specify the rescaled vector",
+   "category": "qmd"
+  },
+  "NP_BAND_PARAL": {
+   "symbol": "NP_BAND_PARAL",
+   "label": "NP_BAND_PARAL",
+   "type": "integer",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_BAND_PARAL: 5",
+   "description": "Number of band groups.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Number of band groups.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "parallelization options"
+  },
+  "NP_DOMAIN_PARAL": {
+   "symbol": "NP_DOMAIN_PARAL",
+   "label": "NP_DOMAIN_PARAL",
+   "type": "integer array",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_DOMAIN_PARAL: 3 3 2",
+   "description": "Dimensions of the 3D Cartesian topology embedded in each band group.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Dimensions of the 3D Cartesian topology embedded in each band group.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "parallelization options"
+  },
+  "NP_DOMAIN_PHI_PARAL": {
+   "symbol": "NP_DOMAIN_PHI_PARAL",
+   "label": "NP_DOMAIN_PHI_PARAL",
+   "type": "integer array",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_DOMAIN_PHI_PARAL: 1 1 2",
+   "description": "Dimensions of the 3D Cartesian topology embedded in the global communicator.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Dimensions of the 3D Cartesian topology embedded in the global communicator.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "parallelization options"
+  },
+  "NP_KPOINT_PARAL": {
+   "symbol": "NP_KPOINT_PARAL",
+   "label": "NP_KPOINT_PARAL",
+   "type": "integer",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_KPOINT_PARAL: 5",
+   "description": "Number of k-point groups.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Number of k-point groups.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "parallelization options"
+  },
+  "NP_SPIN_PARAL": {
+   "symbol": "NP_SPIN_PARAL",
+   "label": "NP_SPIN_PARAL",
+   "type": "integer",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_SPIN_PARAL: 2",
+   "description": "Number of spin groups.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Number of spin groups.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "parallelization options"
+  },
+  "NSTATES": {
+   "symbol": "NSTATES",
+   "label": "NSTATES",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "NSTATES: 24",
+   "description": "The number of Kohn-Sham states/orbitals.",
+   "remark": "This number should not be smaller than half of\nthe total number of valence electrons (N_e) in the system. Note that the number of additional states required increases with increasing values of ELEC_TEMP/SMEARING.",
+   "allow_bool_input": false,
+   "default_remark": "$N_e/2 \\times 1.2 + 5$",
+   "description_raw": "The number of Kohn-Sham states/orbitals.",
+   "remark_raw": "This number should not be smaller than half of\nthe total number of valence electrons ($N_e$) in the system. Note that the number of additional states required increases with increasing values of \\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}}/\\hyperlink{SMEARING}{\\texttt{SMEARING}}.",
+   "category": "system"
+  },
+  "NUM_CHEFSI": {
+   "symbol": "NUM_CHEFSI",
+   "label": "NUM_CHEFSI",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "NUM_CHEFSI: 2",
+   "description": "The number of times ChefSI algorithm is repeated in SCF iteration except the first one, which is controlled by RHO_TRIGGER.",
+   "remark": "For non-collinear spin calculation, it might helped SCF convergence in some cases.",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "The number of times ChefSI algorithm is repeated in SCF iteration except the first one, which is controlled by \\texttt{RHO\\_TRIGGER}.",
+   "remark_raw": "For non-collinear spin calculation, it might helped SCF convergence in some cases.",
+   "category": "scf"
+  },
+  "N_TYPE_ATOM": {
+   "symbol": "N_TYPE_ATOM",
+   "label": "N_TYPE_ATOM",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "N_TYPE_ATOM: 2",
+   "description": "The number of atoms of a ATOM_TYPE specified immediately before this variable.",
+   "remark": "For a system with different types of atoms, one has to specify the number of atoms for every type.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "The number of atoms of a \\hyperlink{ATOM_TYPE}{\\texttt{ATOM\\_TYPE}} specified immediately before this variable.",
+   "remark_raw": "For a system with different types of atoms, one has to specify the number of atoms for every type.",
+   "category": "system"
+  },
+  "OUTPUT_FILE": {
+   "symbol": "OUTPUT_FILE",
+   "label": "OUTPUT_FILE",
+   "type": "string",
+   "default": "Same as the input file name",
+   "unit": "No unit",
+   "example": "OUTPUT_FILE: myfname",
+   "description": "The name of the output files. The output files are attached with a suffix (`.out',`.static',`.geopt' and `.aimd').",
+   "remark": "If an output file with the same name already exist, the results will be written to a file with a number attached, e.g., `myfname.out_1'. The maximum number of output files with the same name allowed is 100. After that the output files will be overwritten in succession.",
+   "allow_bool_input": false,
+   "default_remark": "Same as the input file name",
+   "description_raw": "The name of the output files. The output files are attached with a suffix (`.out',`.static',`.geopt' and `.aimd').",
+   "remark_raw": "If an output file with the same name already exist, the results will be written to a file with a number attached, e.g., `myfname.out\\_1'. The maximum number of output files with the same name allowed is 100. After that the output files will be overwritten in succession.",
+   "category": "print options"
+  },
+  "PRECOND_KERKER_KTF": {
+   "symbol": "PRECOND_KERKER_KTF",
+   "label": "PRECOND_KERKER_KTF",
+   "type": "double",
+   "default": 1.0,
+   "unit": "$\\textrm{Bohr}^{-1}$",
+   "example": "PRECOND_KERKER_KTF: 0.8",
+   "description": "The Thomas-Fermi screening length appearing in the kerker preconditioner (MIXING_PRECOND).",
+   "allow_bool_input": false,
+   "default_remark": "1.0",
+   "description_raw": "The Thomas-Fermi screening length appearing in the \\texttt{kerker} preconditioner (\\hyperlink{MIXING_PRECOND}{\\texttt{MIXING\\_PRECOND}}).",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "PRECOND_KERKER_KTF_MAG": {
+   "symbol": "PRECOND_KERKER_KTF_MAG",
+   "label": "PRECOND_KERKER_KTF_MAG",
+   "type": "double",
+   "default": 1.0,
+   "unit": "$\\textrm{Bohr}^{-1}$",
+   "example": "PRECOND_KERKER_KTF_MAG: 0.8",
+   "description": "The Thomas-Fermi screening length appearing in the kerker preconditioner for the magnetization density (MIXING_PRECOND_MAG).",
+   "allow_bool_input": false,
+   "default_remark": "1.0",
+   "description_raw": "The Thomas-Fermi screening length appearing in the \\texttt{kerker} preconditioner for the magnetization density (\\hyperlink{MIXING_PRECOND_MAG}{\\texttt{MIXING\\_PRECOND\\_MAG}}).",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "PRECOND_KERKER_THRESH": {
+   "symbol": "PRECOND_KERKER_THRESH",
+   "label": "PRECOND_KERKER_THRESH",
+   "type": "double",
+   "default": 0.1,
+   "unit": "No unit",
+   "example": "PRECOND_KERKER_THRESH: 0.0",
+   "description": "The threshold for the kerker preconditioner (MIXING_PRECOND).",
+   "remark": "This threshold will be scaled by the MIXING_PARAMETER. If the threshold is set to 0, the original kerker preconditioner is recovered.",
+   "allow_bool_input": false,
+   "default_remark": "0.1",
+   "description_raw": "The threshold for the \\texttt{kerker} preconditioner (\\hyperlink{MIXING_PRECOND}{\\texttt{MIXING\\_PRECOND}}).",
+   "remark_raw": "This threshold will be scaled by the \\hyperlink{MIXING_PARAMETER}{\\texttt{MIXING\\_PARAMETER}}. If the threshold is set to 0, the original \\texttt{kerker} preconditioner is recovered.",
+   "category": "scf"
+  },
+  "PRECOND_KERKER_THRESH_MAG": {
+   "symbol": "PRECOND_KERKER_THRESH_MAG",
+   "label": "PRECOND_KERKER_THRESH_MAG",
+   "type": "double",
+   "default": 0.1,
+   "unit": "No unit",
+   "example": "PRECOND_KERKER_THRESH_MAG: 0.0",
+   "description": "The threshold for the kerker preconditioner the magnetization density (MIXING_PRECOND_MAG).",
+   "remark": "This threshold will be scaled by the MIXING_PARAMETER_MAG. If the threshold is set to 0, the original kerker preconditioner is recovered.",
+   "allow_bool_input": false,
+   "default_remark": "0.1",
+   "description_raw": "The threshold for the \\texttt{kerker} preconditioner the magnetization density (\\hyperlink{MIXING_PRECOND_MAG}{\\texttt{MIXING\\_PRECOND\\_MAG}}).",
+   "remark_raw": "This threshold will be scaled by the \\hyperlink{MIXING_PARAMETER_MAG}{\\texttt{MIXING\\_PARAMETER\\_MAG}}. If the threshold is set to 0, the original \\texttt{kerker} preconditioner is recovered.",
+   "category": "scf"
+  },
+  "PRINT_ATOMS": {
+   "symbol": "PRINT_ATOMS",
+   "label": "PRINT_ATOMS",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_ATOMS: 1",
+   "description": "Flag for writing the atomic positions. For ground-state calculations, atom positions are printed to a `.static' output file. For structural relaxation calculations, atom positions are printed to a `.geopt' file. For QMD calculations, atom positions are printed to a `.aimd' file.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing the atomic positions. For ground-state calculations, atom positions are printed to a `.static' output file. For structural relaxation calculations, atom positions are printed to a `.geopt' file. For QMD calculations, atom positions are printed to a `.aimd' file.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "print options"
+  },
+  "PRINT_DENSITY": {
+   "symbol": "PRINT_DENSITY",
+   "label": "PRINT_DENSITY",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_DENSITY: 1",
+   "description": "Flag for writing electron density into cube format. For spin-unpolarized calculation, electron density is printed into .dens file. For collinear spin calculation, total, spin-up and spin-down electron density are printed into .dens, .densUp and .densDwn file, respectively.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing electron density into cube format. For spin-unpolarized calculation, electron density is printed into .dens file. For collinear spin calculation, total, spin-up and spin-down electron density are printed into .dens, .densUp and .densDwn file, respectively.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "print options"
+  },
+  "PRINT_EIGEN": {
+   "symbol": "PRINT_EIGEN",
+   "label": "PRINT_EIGEN",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_EIGEN: 1",
+   "description": "Flag for writing eigenvalues and occupations into .eigen file.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing eigenvalues and occupations into .eigen file.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "print options"
+  },
+  "PRINT_ENERGY_DENSITY": {
+   "symbol": "PRINT_ENERGY_DENSITY",
+   "label": "PRINT_ENERGY_DENSITY",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_ENERGY_DENSITY: 1",
+   "description": "Flag for writing a few energy densities into cube format. Currently, only kinetic energy density, exchange correlation energy density (without exact exchange contribution) and exact exchange energy density (if any) are implemented.",
+   "remark": "For spin-unpolarized calculation, kinetic energy density is written into .kedens, exchange correlation energy density is written into .xcedens, and exact exchange energy density is written into .exxedens. For collinear spin calculation, total, spin-up, spin-down kinetic energy density are written into .kedens, kedensUp, kedensDwn files, total, spin-up, spin-down exact exchange energy density are writted into .exxedens, .exxedensUp, .exxedensDwn files, respectively.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing a few energy densities into cube format. Currently, only kinetic energy density, exchange correlation energy density (without exact exchange contribution) and exact exchange energy density (if any) are implemented.",
+   "remark_raw": "For spin-unpolarized calculation, kinetic energy density is written into .kedens, exchange correlation energy density is written into .xcedens, and exact exchange energy density is written into .exxedens. For collinear spin calculation, total, spin-up, spin-down kinetic energy density are written into .kedens, kedensUp, kedensDwn files, total, spin-up, spin-down exact exchange energy density are writted into .exxedens, .exxedensUp, .exxedensDwn files, respectively.",
+   "category": "print options"
+  },
+  "PRINT_FORCES": {
+   "symbol": "PRINT_FORCES",
+   "label": "PRINT_FORCES",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_FORCES: 1",
+   "description": "Flag for writing the atomic forces. For ground-state calculations, forces are printed to a `.static' output file. For structural relaxation calculations, forces are printed to a `.geopt' file. For QMD calculations, forces are printed to a `.aimd' file.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing the atomic forces. For ground-state calculations, forces are printed to a `.static' output file. For structural relaxation calculations, forces are printed to a `.geopt' file. For QMD calculations, forces are printed to a `.aimd' file.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "print options"
+  },
+  "PRINT_MDOUT": {
+   "symbol": "PRINT_MDOUT",
+   "label": "PRINT_MDOUT",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PRINT_MDOUT: 0",
+   "description": "Flag for printing the the QMD output into the .aimd file.",
+   "remark": "",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for printing the the QMD output into the .aimd file.",
+   "remark_raw": "%",
+   "category": "print options"
+  },
+  "PRINT_ORBITAL": {
+   "symbol": "PRINT_ORBITAL",
+   "label": "PRINT_ORBITAL",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PRINT_ORBITAL: 1",
+   "description": "Flag for writing Kohn-Sham orbitals into a binary file.",
+   "remark": "It consists of headers with system information and the orbitals. \nFirst define a few variables and their types. \n\\begintable[]\n\\begintabular|m1.8cm|m2cm|m6.3cm|\n\\hline\nname         & type, length       & description                                                                                                   \n \\hline\nNx Ny Nz     & int, 1             & Number of FD nodes in x,y,z directions                                                                        \n \\hline\nNd           & int, 1             & Total number of FD nodes                                                                                      \n \\hline\ndx, dy,dz    & double, 1         & mesh size in x,y,z directions                                                                                 \n \\hline\ndV           & double, 1          & unit Volume                                                                                                   \n \\hline\nNspinor      & int, 2             & Number of spinor in orbitals                                                                                  \n \\hline\nisReal       & int, 1             & Flag for orbitals being real or complex                                                                       \n \\hline\nnspin        & int, 1             & Number of spin channel printed                                                                                \n \\hline\nnkpt         & int, 1             & Number of k-point printed                                                                                     \n \\hline\n\\endtabular\n\\endtable",
+   "remark -- cont.": "\\begin{algorithm}[H]\n\\begin{algorithmic}\n\\For{ispin = 1:nspin}\n\\For{ikpt = 1:nkpt}\n\\For{iband = 1:nband}\n    \\State spin\\_index, kpt\\_index, kpt\\_vec, band\\_index\n    \\If{isReal == 1}\n        \\State psi\\_real\n    \\Else\n        \\State psi\\_complex\n    \\EndIf\n\\EndFor\n\\EndFor\n\\EndFor\n\\end{algorithmic}\n\\end{algorithm}",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for writing Kohn-Sham orbitals into a binary file.",
+   "remark_raw": "It consists of headers with system information and the orbitals. \nFirst define a few variables and their types. \n\\begin{table}[]\n\\begin{tabular}{|m{1.8cm}|m{2cm}|m{6.3cm}|}\n\\hline\nname         & type, length       & description                                                                                                   \\\\ \\hline\nNx Ny Nz     & int, 1             & Number of FD nodes in x,y,z directions                                                                        \\\\ \\hline\nNd           & int, 1             & Total number of FD nodes                                                                                      \\\\ \\hline\ndx, dy,dz    & double, 1         & mesh size in x,y,z directions                                                                                 \\\\ \\hline\ndV           & double, 1          & unit Volume                                                                                                   \\\\ \\hline\nNspinor      & int, 2             & Number of spinor in orbitals                                                                                  \\\\ \\hline\nisReal       & int, 1             & Flag for orbitals being real or complex                                                                       \\\\ \\hline\nnspin        & int, 1             & Number of spin channel printed                                                                                \\\\ \\hline\nnkpt         & int, 1             & Number of k-point printed                                                                                     \\\\ \\hline\n\\end{tabular}\n\\end{table}",
+   "category": "print options"
+  },
+  "PRINT_RELAXOUT": {
+   "symbol": "PRINT_RELAXOUT",
+   "label": "PRINT_RELAXOUT",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PRINT_RELAXOUT: 0",
+   "description": "Flag for printing the structural relaxation data in a .geopt file.",
+   "remark": "Required only if RELAX_FLAG is greater than 0.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for printing the structural relaxation data in a .geopt file.",
+   "remark_raw": "Required only if \\hyperlink{RELAX_FLAG}{\\texttt{RELAX\\_FLAG}} is greater than $0$.",
+   "category": "print options"
+  },
+  "PRINT_RESTART": {
+   "symbol": "PRINT_RESTART",
+   "label": "PRINT_RESTART",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PRINT_RESTART: 0",
+   "description": "Flag for writing the .restart file, used to restart QMD and structural relaxation simulations.",
+   "remark": "Relevant only if either MD_FLAG is 1 or RELAX_FLAG is 1.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for writing the .restart file, used to restart QMD and structural relaxation simulations.",
+   "remark_raw": "Relevant only if either \\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} is $1$ or \\hyperlink{RELAX_FLAG}{\\texttt{RELAX\\_FLAG}} is $1$.",
+   "category": "print options"
+  },
+  "PRINT_RESTART_FQ": {
+   "symbol": "PRINT_RESTART_FQ",
+   "label": "PRINT_RESTART_FQ",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PRINT_RESTART_FQ: 10",
+   "description": "Frequency at which .restart file is written in QMD and structural optimization simulations.",
+   "remark": "Relevant only if either MD_FLAG is 1 or RELAX_FLAG is 1.",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "Frequency at which .restart file is written in QMD and structural optimization simulations.",
+   "remark_raw": "Relevant only if either \\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} is $1$ or \\hyperlink{RELAX_FLAG}{\\texttt{RELAX\\_FLAG}} is $1$.",
+   "category": "print options"
+  },
+  "PRINT_VELS": {
+   "symbol": "PRINT_VELS",
+   "label": "PRINT_VELS",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PRINT_VELS: 0",
+   "description": "Flag for printing the ion velocities in an QMD simulation into the .aimd file.",
+   "remark": "Relevant only if MD_FLAG is set to 1.",
+   "allow_bool_input": true,
+   "default_remark": "1",
+   "description_raw": "Flag for printing the ion velocities in an QMD simulation into the .aimd file.",
+   "remark_raw": "Relevant only if \\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} is set to $1$.",
+   "category": "print options"
+  },
+  "PSEUDO_POT": {
+   "symbol": "PSEUDO_POT",
+   "label": "PSEUDO_POT",
+   "type": "string",
+   "default": null,
+   "unit": "No unit",
+   "example": "PSEUDO_POT: ../psp/Fe.psp8",
+   "description": "Path to the pseudopotential file.",
+   "remark": "The default directory for the pseudopotential files is the same as the input files. For example, if a pseudopotential Fe.psp8 is put in the same directory as the input files, one can simply specify PSEUDO_POT: Fe.psp8.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Path to the pseudopotential file.",
+   "remark_raw": "The default directory for the pseudopotential files is the same as the input files. For example, if a pseudopotential Fe.psp8 is put in the same directory as the input files, one can simply specify \\texttt{PSEUDO\\_POT: Fe.psp8}.",
+   "category": "system"
+  },
+  "PULAY_FREQUENCY": {
+   "symbol": "PULAY_FREQUENCY",
+   "label": "PULAY_FREQUENCY",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "PULAY_FREQUENCY: 4",
+   "description": "The frequency of Pulay mixing in Periodic Pulay.",
+   "remark": "The default value of 1 corresponds to Pulay mixing.",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "The frequency of Pulay mixing in Periodic Pulay.",
+   "remark_raw": "The default value of 1 corresponds to Pulay mixing.",
+   "category": "scf"
+  },
+  "PULAY_RESTART": {
+   "symbol": "PULAY_RESTART",
+   "label": "PULAY_RESTART",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "PULAY_RESTART: 1",
+   "description": "The flag for restarting the `Periodic Pulay' mixing. If set to 0, the restarted Pulay method is turned off.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "The flag for restarting the `Periodic Pulay' mixing. If set to 0, the restarted Pulay method is turned off.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "scf"
+  },
+  "QMASS": {
+   "symbol": "QMASS",
+   "label": "QMASS",
+   "type": "double",
+   "default": 1653.65493345972,
+   "unit": "atomic unit",
+   "example": "QMASS: 100000",
+   "description": "Gives the inertia factor for Nose Hoover thermostat.",
+   "remark": "Applicable to NVT_NH MD_METHOD only.",
+   "allow_bool_input": false,
+   "default_remark": "1653.654933459720",
+   "description_raw": "Gives the inertia factor for Nose Hoover thermostat.",
+   "remark_raw": "Applicable to NVT\\_NH \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.",
+   "category": "qmd"
+  },
+  "REFERENCE_CUTOFF": {
+   "symbol": "REFERENCE_CUTOFF",
+   "label": "REFERENCE_CUTOFF",
+   "type": "double",
+   "default": 0.5,
+   "unit": "Bohr",
+   "example": "REFERENCE_CUTOFF: 1.0",
+   "description": "The cutoff radius of the reference potential used for evaluating the electrostatic correction arising from overlapping pseudocharge densities.",
+   "remark": "This number should be smaller than half the smallest interatomic distance.",
+   "allow_bool_input": false,
+   "default_remark": "0.5",
+   "description_raw": "The cutoff radius of the reference potential used for evaluating the electrostatic correction arising from overlapping pseudocharge densities.",
+   "remark_raw": "This number should be smaller than half the smallest interatomic distance.",
+   "category": "electrostatics"
+  },
+  "RELAX": {
+   "symbol": "RELAX",
+   "label": "RELAX",
+   "type": "integer array",
+   "default": [
+    1,
+    1,
+    1
+   ],
+   "unit": "No unit",
+   "example": "RELAX: 1 0 1 \n0 1 0",
+   "description": "Atomic coordinate with the corresponding RELAX value 0 is held fixed during relaxation/QMD.",
+   "allow_bool_input": true,
+   "default_remark": "1 1 1",
+   "description_raw": "Atomic coordinate with the corresponding \\texttt{RELAX} value 0 is held fixed during relaxation/QMD.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "system"
+  },
+  "RELAX_FLAG": {
+   "symbol": "RELAX_FLAG",
+   "label": "RELAX_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "RELAX_FLAG: 1",
+   "description": "Flag for performing structural relaxation. 0 means no structural relaxation. 1 represents relaxation of atom positions. 2 represents optimization of volume with the fractional coordinates of the atoms fixed. 3 represents full optimization of the cell i.e., both atoms and cell volume are relaxed",
+   "remark": "This flag should not be specified if MD_FLAG is set to 1.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for performing structural relaxation. $0$ means no structural relaxation. $1$ represents relaxation of atom positions. $2$ represents optimization of volume with the fractional coordinates of the atoms fixed. $3$ represents full optimization of the cell i.e., both atoms and cell volume are relaxed",
+   "remark_raw": "This flag should not be specified if \\hyperlink{MD_FLAG}{\\texttt{MD\\_FLAG}} is set to $1$.",
+   "category": "structural relaxation"
+  },
+  "RELAX_MAXDILAT": {
+   "symbol": "RELAX_MAXDILAT",
+   "label": "RELAX_MAXDILAT",
+   "type": "double",
+   "default": 1.06,
+   "unit": "No unit",
+   "example": "RELAX_MAXDILAT: 1.4",
+   "description": "The maximum scaling of the volume allowed with respect to the initial volume defined by CELL and LATVEC. This will determine the upper-bound and lower-bound in the bisection method (Brent's method) for the volume optimization.",
+   "remark": "",
+   "allow_bool_input": false,
+   "default_remark": "1.06",
+   "description_raw": "The maximum scaling of the volume allowed with respect to the initial volume defined by \\hyperlink{CELL}{\\texttt{CELL}} and \\hyperlink{LATVEC}{\\texttt{LATVEC}}. This will determine the upper-bound and lower-bound in the bisection method (Brent's method) for the volume optimization.",
+   "remark_raw": "%",
+   "category": "structural relaxation"
+  },
+  "RELAX_METHOD": {
+   "symbol": "RELAX_METHOD",
+   "label": "RELAX_METHOD",
+   "type": "string",
+   "default": "LBFGS",
+   "unit": "No unit",
+   "example": "RELAX_METHOD: NLCG",
+   "description": "Specifies the algorithm for structural relaxation. The choices are `LBFGS' (limited-memory BFGS), `NLCG' (Non-linear conjugate gradient), and `FIRE' (Fast inertial relaxation engine).",
+   "remark": "LBFGS is typically the best choice.",
+   "allow_bool_input": false,
+   "default_remark": "LBFGS",
+   "description_raw": "Specifies the algorithm for structural relaxation. The choices are `LBFGS' (limited-memory BFGS), `NLCG' (Non-linear conjugate gradient), and `FIRE' (Fast inertial relaxation engine).",
+   "remark_raw": "LBFGS is typically the best choice.",
+   "category": "structural relaxation"
+  },
+  "RELAX_NITER": {
+   "symbol": "RELAX_NITER",
+   "label": "RELAX_NITER",
+   "type": "integer",
+   "default": 300,
+   "unit": "No unit",
+   "example": "RELAX_NITER: 25",
+   "description": "Specifies the maximum number of iterations for the structural relaxation (RELAX_FLAG).",
+   "remark": "If RESTART_FLAG is set to 1, then relaxation will restart from the last atomic configuration and run for maximum of RELAX_NITER iterations.",
+   "allow_bool_input": false,
+   "default_remark": "300",
+   "description_raw": "Specifies the maximum number of iterations for the structural relaxation (\\hyperlink{RELAX_FLAG}{\\texttt{RELAX\\_FLAG}}).",
+   "remark_raw": "If \\hyperlink{RESTART_FLAG}{\\texttt{RESTART\\_FLAG}} is set to $1$, then relaxation will restart from the last atomic configuration and run for maximum of \\hyperlink{RELAX_NITER}{\\texttt{RELAX\\_NITER}} iterations.",
+   "category": "structural relaxation"
+  },
+  "RESTART_FLAG": {
+   "symbol": "RESTART_FLAG",
+   "label": "RESTART_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "RESTART_FLAG: 0",
+   "description": "Flag for restarting quantum molecular dynamics and structural relaxation. Stores last three histories for quantum molecular dynamics simulations in .restart, .restart-0 and .restart-1 files, respectively.",
+   "remark": "Restarts from the previous configuration which is stored in a .restart file. Currently, code provides restart feature for atomic relaxation and QMD only.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag for restarting quantum molecular dynamics and structural relaxation. Stores last three histories for quantum molecular dynamics simulations in .restart, .restart-$0$ and .restart-$1$ files, respectively.",
+   "remark_raw": "Restarts from the previous configuration which is stored in a .restart file. Currently, code provides restart feature for atomic relaxation and QMD only.",
+   "category": "qmd"
+  },
+  "RHO_TRIGGER": {
+   "symbol": "RHO_TRIGGER",
+   "label": "RHO_TRIGGER",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "RHO_TRIGGER: 5",
+   "description": "The number of times Chebyshev filtering is repeated before updating the electron density in the very first SCF iteration.",
+   "remark": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations. For non-collinear spin calculation, default is 6 otherwise 4.",
+   "allow_bool_input": false,
+   "default_remark": "4 or 6",
+   "description_raw": "The number of times Chebyshev filtering is repeated before updating the electron density in the very first SCF iteration.",
+   "remark_raw": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations. For non-collinear spin calculation, default is 6 otherwise 4.",
+   "category": "scf"
+  },
+  "SCF_ENERGY_ACC": {
+   "symbol": "SCF_ENERGY_ACC",
+   "label": "SCF_ENERGY_ACC",
+   "type": "double",
+   "default": null,
+   "unit": "Ha/atom",
+   "example": "SCF_ENERGY_ACC: 1e-5",
+   "description": "The tolerance on the free energy for the convergence of the SCF iteration.",
+   "remark": "Only one of SCF_ENERGY_ACC, SCF_FORCE_ACC, or TOL_SCF can be specified.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "The tolerance on the free energy for the convergence of the SCF iteration.",
+   "remark_raw": "Only one of \\hyperlink{SCF_ENERGY_ACC}{\\texttt{SCF\\_ENERGY\\_ACC}}, \\hyperlink{SCF_FORCE_ACC}{\\texttt{SCF\\_FORCE\\_ACC}}, or \\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}} can be specified.",
+   "category": "scf"
+  },
+  "SCF_FORCE_ACC": {
+   "symbol": "SCF_FORCE_ACC",
+   "label": "SCF_FORCE_ACC",
+   "type": "double",
+   "default": null,
+   "unit": "Ha/Bohr",
+   "example": "SCF_FORCE_ACC: 1e-4",
+   "description": "The tolerance on the atomic forces for convergence of the SCF iteration.",
+   "remark": "Only one of SCF_FORCE_ACC, TOL_SCF or SCF_ENERGY_ACC can be specified.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "The tolerance on the atomic forces for convergence of the SCF iteration.",
+   "remark_raw": "Only one of \\hyperlink{SCF_FORCE_ACC}{\\texttt{SCF\\_FORCE\\_ACC}}, \\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}} or \\hyperlink{SCF_ENERGY_ACC}{\\texttt{SCF\\_ENERGY\\_ACC}} can be specified.",
+   "category": "scf"
+  },
+  "SMEARING": {
+   "symbol": "SMEARING",
+   "label": "SMEARING",
+   "type": "double",
+   "default": null,
+   "unit": "Ha",
+   "example": "SMEARING: 0.001",
+   "description": "Value of smearing.",
+   "remark": "Equivalent to setting ELEC_TEMP (0.001 Ha = 315.775 Kelvin).",
+   "allow_bool_input": false,
+   "default_remark": "0.007350 for \\texttt{gaussian} \\\\\n0.003675 for \\texttt{fermi-dirac}",
+   "description_raw": "Value of smearing.",
+   "remark_raw": "Equivalent to setting \\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}} (0.001 Ha = 315.775 Kelvin).",
+   "category": "system"
+  },
+  "SPIN": {
+   "symbol": "SPIN",
+   "label": "SPIN",
+   "type": "double array",
+   "default": [
+    0.0
+   ],
+   "unit": "No unit",
+   "example": "SPIN: 0 0 1.0 \n0 0 -1.0",
+   "description": "Specifies the net initial spin on each atom for a spin-polarized calculation. If collinear spin used, user could use either 1 column of data for z-direction of each atom, or 3 columns of data with 0 on the first 2 columns (x,y-directions). For noncollinear spin, use need to use 3 columns of data for all directions.",
+   "allow_bool_input": false,
+   "default_remark": "0.0",
+   "description_raw": "Specifies the net initial spin on each atom for a spin-polarized calculation. If collinear spin used, user could use either 1 column of data for z-direction of each atom, or 3 columns of data with 0 on the first 2 columns (x,y-directions). For noncollinear spin, use need to use 3 columns of data for all directions.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "system"
+  },
+  "SPIN_TYP": {
+   "symbol": "SPIN_TYP",
+   "label": "SPIN_TYP",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "SPIN_TYP: 1",
+   "description": "SPIN_TYP: 0 performs spin unpolarized calculation. \nSPIN_TYP: 1 performs unconstrained collinear spin-polarized calculation.  \nSPIN_TYP: 1 performs unconstrained noncollinear spin-polarized calculation.",
+   "remark": "SPIN_TYP can only take values 0, 1, 2. For collinear calculation, non-relativistic pseudopotential need to be used. For noncollinear calculation, fully relativistic pseudopotentiail need to be used.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "\\texttt{SPIN\\_TYP}: 0 performs spin unpolarized calculation. \\\\\n\\texttt{SPIN\\_TYP}: 1 performs unconstrained collinear spin-polarized calculation.  \\\\\n\\texttt{SPIN\\_TYP}: 1 performs unconstrained noncollinear spin-polarized calculation.",
+   "remark_raw": "\\texttt{SPIN\\_TYP} can only take values 0, 1, 2. For collinear calculation, non-relativistic pseudopotential need to be used. For noncollinear calculation, fully relativistic pseudopotentiail need to be used.",
+   "category": "system"
+  },
+  "TARGET_PRESSURE": {
+   "symbol": "TARGET_PRESSURE",
+   "label": "TARGET_PRESSURE",
+   "type": "double",
+   "default": 0.0,
+   "unit": "GPa",
+   "example": "TARGET_PRESSURE: 40.9611",
+   "description": "Gives the outer pressure in NPT_NH and NPT_NP.",
+   "remark": "Applicable to NPT_NH and NPT_NP MD_METHOD only.",
+   "allow_bool_input": false,
+   "default_remark": "0.0",
+   "description_raw": "Gives the outer pressure in NPT\\_NH and NPT\\_NP.",
+   "remark_raw": "Applicable to NPT\\_NH and NPT\\_NP \\hyperlink{MD_METHOD}{\\texttt{MD\\_METHOD}} only.",
+   "category": "qmd"
+  },
+  "TOL_FOCK": {
+   "symbol": "TOL_FOCK",
+   "label": "TOL_FOCK",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_FOCK: 1e-6",
+   "description": "The tolerance on the Hartree-Fock outer loop, measured by the exact exchange energy difference per atom in 2 consecutive outer loops.",
+   "remark": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "allow_bool_input": false,
+   "default_remark": "$0.2*$\\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}}",
+   "description_raw": "The tolerance on the Hartree-Fock outer loop, measured by the exact exchange energy difference per atom in 2 consecutive outer loops.",
+   "remark_raw": "Only active when using hybrid functionals, like PBE0 and HSE.",
+   "category": "scf"
+  },
+  "TOL_LANCZOS": {
+   "symbol": "TOL_LANCZOS",
+   "label": "TOL_LANCZOS",
+   "type": "double",
+   "default": 0.01,
+   "unit": "No unit",
+   "example": "TOL_LANCZOS: 1e-3",
+   "description": "The tolerance within the Lanczos algorithm for calculating the extremal eigenvalues of the Hamiltonian, required as part of the CheFSI method.",
+   "remark": "Typically, the Lanczos tolerance does not need to be very strict.",
+   "allow_bool_input": false,
+   "default_remark": "1e-2",
+   "description_raw": "The tolerance within the Lanczos algorithm for calculating the extremal eigenvalues of the Hamiltonian, required as part of the CheFSI method.",
+   "remark_raw": "Typically, the Lanczos tolerance does not need to be very strict.",
+   "category": "scf"
+  },
+  "TOL_POISSON": {
+   "symbol": "TOL_POISSON",
+   "label": "TOL_POISSON",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_POISSON: 1e-6",
+   "description": "The tolerance on the norm of the relative residual for the Poisson equation.",
+   "remark": "The tolerance for poisson solver should not be worse than TOL_SCF, otherwise it might seriously affect  the convergence of the SCF iteration.",
+   "allow_bool_input": false,
+   "default_remark": "\\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}}$\\times 0.01$",
+   "description_raw": "The tolerance on the norm of the relative residual for the Poisson equation.",
+   "remark_raw": "The tolerance for poisson solver should not be worse than \\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}}, otherwise it might seriously affect  the convergence of the SCF iteration.",
+   "category": "electrostatics"
+  },
+  "TOL_PRECOND": {
+   "symbol": "TOL_PRECOND",
+   "label": "TOL_PRECOND",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_PRECOND: 1e-4",
+   "description": "The tolerance on the relative residual for the linear systems arising during the real-space preconditioning of the SCF.",
+   "remark": "The linear systems do not need to be solved very accurately. h is the mesh spacing.",
+   "allow_bool_input": false,
+   "default_remark": "$h^2\\times0.001$",
+   "description_raw": "The tolerance on the relative residual for the linear systems arising during the real-space preconditioning of the SCF.",
+   "remark_raw": "The linear systems do not need to be solved very accurately. $h$ is the mesh spacing.",
+   "category": "scf"
+  },
+  "TOL_PSEUDOCHARGE": {
+   "symbol": "TOL_PSEUDOCHARGE",
+   "label": "TOL_PSEUDOCHARGE",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_PSEUDOCHARGE: 1e-6",
+   "description": "The error in the net enclosed charge for the pseudocharge density of each atom.",
+   "allow_bool_input": false,
+   "default_remark": "\\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}}$\\times 0.001$",
+   "description_raw": "The error in the net enclosed charge for the pseudocharge density of each atom.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "electrostatics"
+  },
+  "TOL_RELAX": {
+   "symbol": "TOL_RELAX",
+   "label": "TOL_RELAX",
+   "type": "double",
+   "default": 0.0005,
+   "unit": "Ha/Bohr",
+   "example": "TOL_RELAX: 1e-3",
+   "description": "Specifies the tolerance for termination of the structural relaxation. The tolerance is defined on the maximum force component (in absolute sense) over all atoms.",
+   "remark": "",
+   "allow_bool_input": false,
+   "default_remark": "5e-4",
+   "description_raw": "Specifies the tolerance for termination of the structural relaxation. The tolerance is defined on the maximum force component (in absolute sense) over all atoms.",
+   "remark_raw": "%",
+   "category": "structural relaxation"
+  },
+  "TOL_RELAX_CELL": {
+   "symbol": "TOL_RELAX_CELL",
+   "label": "TOL_RELAX_CELL",
+   "type": "double",
+   "default": 0.01,
+   "unit": "GPa",
+   "example": "TOL_RELAX: 1e-3",
+   "description": "Specifies the tolerance for termination of the cell relaxation. The tolerance is defined on the maximum principle stress component.",
+   "remark": "",
+   "allow_bool_input": false,
+   "default_remark": "1e-2",
+   "description_raw": "Specifies the tolerance for termination of the cell relaxation. The tolerance is defined on the maximum principle stress component.",
+   "remark_raw": "%",
+   "category": "structural relaxation"
+  },
+  "TOL_SCF": {
+   "symbol": "TOL_SCF",
+   "label": "TOL_SCF",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_SCF: 1e-5",
+   "description": "In case of single point calculation, TOL_SCF is set for 10^-5 Ha/atom energy accuracy.\nIn case of MD, TOL_SCF is set for 10^-3 Ha/Bohr force accuracy.\nIn case of relaxation, TOL_SCF is set for TOL_RELAX/5 Ha/Bohr force accuracy. \nThe tolerance on the normalized residual of the effective potential or the electron density for convergence of the SCF iteration.",
+   "remark": "Only one of TOL_SCF, SCF_ENERGY_ACC, or SCF_FORCE_ACC can be specified.",
+   "allow_bool_input": false,
+   "default_remark": "see description",
+   "description_raw": "In case of single point calculation, \\texttt{TOL\\_SCF} is set for $10^{-5}$ Ha/atom energy accuracy.\nIn case of MD, \\texttt{TOL\\_SCF} is set for $10^{-3}$ Ha/Bohr force accuracy.\nIn case of relaxation, \\texttt{TOL\\_SCF} is set for \\hyperlink{TOL_RELAX}{\\texttt{TOL\\_RELAX}}/5 Ha/Bohr force accuracy. \\\\\nThe tolerance on the normalized residual of the effective potential or the electron density for convergence of the SCF iteration.",
+   "remark_raw": "Only one of \\hyperlink{TOL_SCF}{\\texttt{TOL\\_SCF}}, \\hyperlink{SCF_ENERGY_ACC}{\\texttt{SCF\\_ENERGY\\_ACC}}, or \\hyperlink{SCF_FORCE_ACC}{\\texttt{SCF\\_FORCE\\_ACC}} can be specified.",
+   "category": "scf"
+  },
+  "TOL_SCF_INIT": {
+   "symbol": "TOL_SCF_INIT",
+   "label": "TOL_SCF_INIT",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "TOL_SCF_INIT: 1e-6",
+   "description": "The initial SCF tolerance for PBE iteration when using hybrid functionals.",
+   "remark": "Only active when using hybrid functionals, like PBE0 and HSE. Change the TOL_SCF_INIT to change the initial guess for Hartree Fock outer loop.",
+   "allow_bool_input": false,
+   "default_remark": "$\\max($\\hyperlink{TOL_FOCK}{\\texttt{TOL\\_FOCK}}$\\times 10,0.001)$",
+   "description_raw": "The initial SCF tolerance for PBE iteration when using hybrid functionals.",
+   "remark_raw": "Only active when using hybrid functionals, like PBE0 and HSE. Change the \\texttt{TOL\\_SCF\\_INIT} to change the initial guess for Hartree Fock outer loop.",
+   "category": "scf"
+  },
+  "TWTIME": {
+   "symbol": "TWTIME",
+   "label": "TWTIME",
+   "type": "double",
+   "default": 1000000000.0,
+   "unit": "min",
+   "example": "TWTIME: 1000",
+   "description": "Gives the upper bound on the wall time for quantum molecular dynamics.",
+   "allow_bool_input": false,
+   "default_remark": "1e9",
+   "description_raw": "Gives the upper bound on the wall time for quantum molecular dynamics.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "qmd"
+  },
+  "TWIST_ANGLE": {
+   "symbol": "TWIST_ANGLE",
+   "label": "TWIST_ANGLE",
+   "type": "double",
+   "default": 0.0,
+   "unit": "rad/Bohr",
+   "example": "TWIST_ANGLE: 0.0045",
+   "description": "External twist per unit length applied on the nanotube.",
+   "remark": "If using helical symmetry (D C H), we also have to add the intrinsic twist.",
+   "allow_bool_input": false,
+   "default_remark": "0",
+   "description_raw": "External twist per unit length applied on the nanotube.",
+   "remark_raw": "If using helical symmetry (D C H), we also have to add the intrinsic twist.",
+   "category": "cyclix"
+  },
+  "NP_DOMAIN_SQ_PARAL": {
+   "symbol": "NP_DOMAIN_SQ_PARAL",
+   "label": "NP_DOMAIN_SQ_PARAL",
+   "type": "integer array",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_DOMAIN_SQ_PARAL: 3 3 2",
+   "description": "Dimensions of the 3D Cartesian topology for SQ method.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Dimensions of the 3D Cartesian topology for SQ method.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "spectral quadrature"
+  },
+  "SQ_FLAG": {
+   "symbol": "SQ_FLAG",
+   "label": "SQ_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "SQ_FLAG: 1",
+   "description": "Flag to turn on SQ method",
+   "remark": "SQ method can not be turned on simultaneously with CS, SQ3, hybrid functionals.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag to turn on SQ method",
+   "remark_raw": "SQ method can not be turned on simultaneously with CS, SQ3, hybrid functionals.",
+   "category": "spectral quadrature"
+  },
+  "SQ_GAUSS_MEM": {
+   "symbol": "SQ_GAUSS_MEM",
+   "label": "SQ_GAUSS_MEM",
+   "type": "string",
+   "default": "LOW",
+   "unit": "No unit",
+   "example": "SQ_GAUSS_MEM: HIGH",
+   "description": "Flag for memory option when using Gauss quadrature for density matrix.",
+   "allow_bool_input": false,
+   "default_remark": "LOW",
+   "description_raw": "Flag for memory option when using Gauss quadrature for density matrix.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "spectral quadrature"
+  },
+  "SQ_NPL_G": {
+   "symbol": "SQ_NPL_G",
+   "label": "SQ_NPL_G",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "SQ_NPL_G: 24",
+   "description": "Degree of polynomial for Gauss Quadrature.",
+   "remark": "SQ_NPL_G must be specified if SQ is turned on.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Degree of polynomial for Gauss Quadrature.",
+   "remark_raw": "\\texttt{SQ\\_NPL\\_G} must be specified if SQ is turned on.",
+   "category": "spectral quadrature"
+  },
+  "SQ_RCUT": {
+   "symbol": "{SQ_RCUT",
+   "label": "SQ_RCUT",
+   "type": "double",
+   "default": null,
+   "unit": "Bohr",
+   "example": "SQ_RCUT: 2.0",
+   "description": "Truncation or localization radius",
+   "remark": "SQ_RCUT must be specified if SQ is turned on.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Truncation or localization radius",
+   "remark_raw": "\\texttt{{SQ\\_RCUT}} must be specified if SQ is turned on.",
+   "category": "spectral quadrature"
+  },
+  "SQ_TOL_OCC": {
+   "symbol": "SQ_TOL_OCC",
+   "label": "SQ_TOL_OCC",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "SQ_TOL_OCC: 1E-5",
+   "description": "Tolerance for occupation corresponding to maximum eigenvalue.",
+   "allow_bool_input": false,
+   "default_remark": "$10^{-6}$",
+   "description_raw": "Tolerance for occupation corresponding to maximum eigenvalue.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "spectral quadrature"
+  }
+ },
+ "other_parameters": {},
+ "data_types": [
+  "integer array",
+  "string",
+  "double",
+  "other",
+  "integer",
+  "double array"
+ ]
+}
diff --git a/sparc/sparc_json_api/parameters.json b/sparc/sparc_json_api/parameters.json
index 42ff9158..1f44fb0d 100644
--- a/sparc/sparc_json_api/parameters.json
+++ b/sparc/sparc_json_api/parameters.json
@@ -1,5 +1,5 @@
 {
- "sparc_version": "2023.07.10",
+ "sparc_version": "2023.09.05",
  "categories": [
   "system",
   "scf",
@@ -318,12 +318,12 @@
    "default": null,
    "unit": "No unit",
    "example": "EXCHANGE_CORRELATION: LDA_PW",
-   "description": "Choice of exchange-correlation functional. Options are LDA_PW (Perdew-Wang LDA), LDA_PZ (Purdew-Zunger LDA), \nGGA_PBE (PBE GGA), GGA_RPBE (revised PBE GGA), and GGA_PBEsol (PBE GGA revised for solids), \nPBE0, HF (Hartree-Fock), HSE,\nvdWDF1 (van der Waals Density Functional developed by Dion et al.), vdWDF2 (vdW Density Functional modified by Lee et al), and SCAN (SCAN metaGGA), and R2SCAN (r2SCAN metaGGA).",
-   "remark": "For spin-polarized calculation (SPIN_TYP = 1), LDA_PZ is not available.\nCurrently SCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
+   "description": "Choice of exchange-correlation functional. Options are LDA_PW (Perdew-Wang LDA), LDA_PZ (Purdew-Zunger LDA), \nGGA_PBE (PBE GGA), GGA_RPBE (revised PBE GGA), and GGA_PBEsol (PBE GGA revised for solids), \nPBE0, HF (Hartree-Fock), HSE,\nvdWDF1 (van der Waals Density Functional developed by Dion et al.), vdWDF2 (vdW Density Functional modified by Lee et al), SCAN (SCAN metaGGA), RSCAN (rSCAN metaGGA), and R2SCAN (r2SCAN metaGGA).",
+   "remark": "For spin-polarized calculation (SPIN_TYP = 1), LDA_PZ is not available.\nCurrently SCAN, RSCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
    "allow_bool_input": false,
    "default_remark": "No Default",
-   "description_raw": "Choice of exchange-correlation functional. Options are \\texttt{LDA\\_PW} (Perdew-Wang LDA), \\texttt{LDA\\_PZ} (Purdew-Zunger LDA), \n\\texttt{GGA\\_PBE} (PBE GGA), \\texttt{GGA\\_RPBE} (revised PBE GGA), and \\texttt{GGA\\_PBEsol} (PBE GGA revised for solids), \n\\texttt{PBE0}, \\texttt{HF} (Hartree-Fock), \\texttt{HSE},\n\\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \\texttt{vdWDF2} (vdW Density Functional modified by Lee et al), and \\texttt{SCAN} (SCAN metaGGA), and \\texttt{R2SCAN} (r2SCAN metaGGA).",
-   "remark_raw": "For spin-polarized calculation (\\hyperlink{SPIN_TYP}{\\texttt{SPIN\\_TYP}} = 1), \\texttt{LDA\\_PZ} is not available.\n\nCurrently SCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
+   "description_raw": "Choice of exchange-correlation functional. Options are \\texttt{LDA\\_PW} (Perdew-Wang LDA), \\texttt{LDA\\_PZ} (Purdew-Zunger LDA), \n\\texttt{GGA\\_PBE} (PBE GGA), \\texttt{GGA\\_RPBE} (revised PBE GGA), and \\texttt{GGA\\_PBEsol} (PBE GGA revised for solids), \n\\texttt{PBE0}, \\texttt{HF} (Hartree-Fock), \\texttt{HSE},\n\\texttt{vdWDF1} (van der Waals Density Functional developed by Dion et al.), \\texttt{vdWDF2} (vdW Density Functional modified by Lee et al), \\texttt{SCAN} (SCAN metaGGA), \\texttt{RSCAN} (rSCAN metaGGA), and \\texttt{R2SCAN} (r2SCAN metaGGA).",
+   "remark_raw": "For spin-polarized calculation (\\hyperlink{SPIN_TYP}{\\texttt{SPIN\\_TYP}} = 1), \\texttt{LDA\\_PZ} is not available.\n\nCurrently SCAN, RSCAN and R2SCAN does not support nonlinear core correction pseudopotential.",
    "category": "system"
   },
   "EXX_ACE_VALENCE_STATES": {
@@ -1118,7 +1118,7 @@
   "NPT_SCALE_VECS": {
    "symbol": "NPT_SCALE_VECS",
    "label": "NPT_SCALE_VECS",
-   "type": "integer",
+   "type": "integer array",
    "default": 1,
    "unit": "No unit",
    "example": "NPT_SCALE_VECS: 2",
@@ -1220,6 +1220,21 @@
    "remark_raw": "This number should not be smaller than half of\nthe total number of valence electrons ($N_e$) in the system. Note that the number of additional states required increases with increasing values of \\hyperlink{ELEC_TEMP}{\\texttt{ELEC\\_TEMP}}/\\hyperlink{SMEARING}{\\texttt{SMEARING}}.",
    "category": "system"
   },
+  "NUM_CHEFSI": {
+   "symbol": "NUM_CHEFSI",
+   "label": "NUM_CHEFSI",
+   "type": "integer",
+   "default": 1,
+   "unit": "No unit",
+   "example": "NUM_CHEFSI: 2",
+   "description": "The number of times ChefSI algorithm is repeated in SCF iteration except the first one, which is controlled by RHO_TRIGGER.",
+   "remark": "For non-collinear spin calculation, it might helped SCF convergence in some cases.",
+   "allow_bool_input": false,
+   "default_remark": "1",
+   "description_raw": "The number of times ChefSI algorithm is repeated in SCF iteration except the first one, which is controlled by \\texttt{RHO\\_TRIGGER}.",
+   "remark_raw": "For non-collinear spin calculation, it might helped SCF convergence in some cases.",
+   "category": "scf"
+  },
   "N_TYPE_ATOM": {
    "symbol": "N_TYPE_ATOM",
    "label": "N_TYPE_ATOM",
@@ -1649,15 +1664,15 @@
    "symbol": "RHO_TRIGGER",
    "label": "RHO_TRIGGER",
    "type": "integer",
-   "default": 4,
+   "default": null,
    "unit": "No unit",
    "example": "RHO_TRIGGER: 5",
    "description": "The number of times Chebyshev filtering is repeated before updating the electron density in the very first SCF iteration.",
-   "remark": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations.",
+   "remark": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations. For non-collinear spin calculation, default is 6 otherwise 4.",
    "allow_bool_input": false,
-   "default_remark": "4",
+   "default_remark": "4 or 6",
    "description_raw": "The number of times Chebyshev filtering is repeated before updating the electron density in the very first SCF iteration.",
-   "remark_raw": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations.",
+   "remark_raw": "Values smaller than the default value of 4 can result in a significant increase in the number of SCF\niterations. Larger values can sometimes reduce the number of SCF iterations. For non-collinear spin calculation, default is 6 otherwise 4.",
    "category": "scf"
   },
   "SCF_ENERGY_ACC": {
@@ -1713,11 +1728,11 @@
     0.0
    ],
    "unit": "No unit",
-   "example": "SPIN: 1.0 \n-1.0",
-   "description": "Specifies the net initial spin on each atom for a spin-polarized calculation.",
+   "example": "SPIN: 0 0 1.0 \n0 0 -1.0",
+   "description": "Specifies the net initial spin on each atom for a spin-polarized calculation. If collinear spin used, user could use either 1 column of data for z-direction of each atom, or 3 columns of data with 0 on the first 2 columns (x,y-directions). For noncollinear spin, use need to use 3 columns of data for all directions.",
    "allow_bool_input": false,
    "default_remark": "0.0",
-   "description_raw": "Specifies the net initial spin on each atom for a spin-polarized calculation.",
+   "description_raw": "Specifies the net initial spin on each atom for a spin-polarized calculation. If collinear spin used, user could use either 1 column of data for z-direction of each atom, or 3 columns of data with 0 on the first 2 columns (x,y-directions). For noncollinear spin, use need to use 3 columns of data for all directions.",
    "remark_raw": "",
    "remark": "",
    "category": "system"
@@ -1729,12 +1744,12 @@
    "default": 0,
    "unit": "No unit",
    "example": "SPIN_TYP: 1",
-   "description": "SPIN_TYP: 0 performs spin unpolarized calculation. \nSPIN_TYP: 1 performs unconstrained collinear spin-polarized calculation.",
-   "remark": "SPIN_TYP can only take values 0 and 1.",
+   "description": "SPIN_TYP: 0 performs spin unpolarized calculation. \nSPIN_TYP: 1 performs unconstrained collinear spin-polarized calculation.  \nSPIN_TYP: 1 performs unconstrained noncollinear spin-polarized calculation.",
+   "remark": "SPIN_TYP can only take values 0, 1, 2. For collinear calculation, non-relativistic pseudopotential need to be used. For noncollinear calculation, fully relativistic pseudopotentiail need to be used.",
    "allow_bool_input": true,
    "default_remark": "0",
-   "description_raw": "\\texttt{SPIN\\_TYP}: 0 performs spin unpolarized calculation. \\\\\n\\texttt{SPIN\\_TYP}: 1 performs unconstrained collinear spin-polarized calculation.",
-   "remark_raw": "\\texttt{SPIN\\_TYP} can only take values 0 and 1.",
+   "description_raw": "\\texttt{SPIN\\_TYP}: 0 performs spin unpolarized calculation. \\\\\n\\texttt{SPIN\\_TYP}: 1 performs unconstrained collinear spin-polarized calculation.  \\\\\n\\texttt{SPIN\\_TYP}: 1 performs unconstrained noncollinear spin-polarized calculation.",
+   "remark_raw": "\\texttt{SPIN\\_TYP} can only take values 0, 1, 2. For collinear calculation, non-relativistic pseudopotential need to be used. For noncollinear calculation, fully relativistic pseudopotentiail need to be used.",
    "category": "system"
   },
   "TARGET_PRESSURE": {
@@ -1901,15 +1916,120 @@
    "remark_raw": "",
    "remark": "",
    "category": "qmd"
+  },
+  "TWIST_ANGLE": {
+   "symbol": "TWIST_ANGLE",
+   "label": "TWIST_ANGLE",
+   "type": "double",
+   "default": 0.0,
+   "unit": "rad/Bohr",
+   "example": "TWIST_ANGLE: 0.0045",
+   "description": "External twist per unit length applied on the nanotube.",
+   "remark": "If using helical symmetry (D C H), we also have to add the intrinsic twist.",
+   "allow_bool_input": false,
+   "default_remark": "0",
+   "description_raw": "External twist per unit length applied on the nanotube.",
+   "remark_raw": "If using helical symmetry (D C H), we also have to add the intrinsic twist.",
+   "category": "cyclix"
+  },
+  "NP_DOMAIN_SQ_PARAL": {
+   "symbol": "NP_DOMAIN_SQ_PARAL",
+   "label": "NP_DOMAIN_SQ_PARAL",
+   "type": "integer array",
+   "default": "auto",
+   "unit": "No unit",
+   "example": "NP_DOMAIN_SQ_PARAL: 3 3 2",
+   "description": "Dimensions of the 3D Cartesian topology for SQ method.",
+   "remark": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "allow_bool_input": false,
+   "default_remark": "Automatically optimized",
+   "description_raw": "Dimensions of the 3D Cartesian topology for SQ method.",
+   "remark_raw": "This option is for development purpose. It's better to let SPARC choose the parallization parameters in practice.",
+   "category": "spectral quadrature"
+  },
+  "SQ_FLAG": {
+   "symbol": "SQ_FLAG",
+   "label": "SQ_FLAG",
+   "type": "integer",
+   "default": 0,
+   "unit": "No unit",
+   "example": "SQ_FLAG: 1",
+   "description": "Flag to turn on SQ method",
+   "remark": "SQ method can not be turned on simultaneously with CS, SQ3, hybrid functionals.",
+   "allow_bool_input": true,
+   "default_remark": "0",
+   "description_raw": "Flag to turn on SQ method",
+   "remark_raw": "SQ method can not be turned on simultaneously with CS, SQ3, hybrid functionals.",
+   "category": "spectral quadrature"
+  },
+  "SQ_GAUSS_MEM": {
+   "symbol": "SQ_GAUSS_MEM",
+   "label": "SQ_GAUSS_MEM",
+   "type": "string",
+   "default": "LOW",
+   "unit": "No unit",
+   "example": "SQ_GAUSS_MEM: HIGH",
+   "description": "Flag for memory option when using Gauss quadrature for density matrix.",
+   "allow_bool_input": false,
+   "default_remark": "LOW",
+   "description_raw": "Flag for memory option when using Gauss quadrature for density matrix.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "spectral quadrature"
+  },
+  "SQ_NPL_G": {
+   "symbol": "SQ_NPL_G",
+   "label": "SQ_NPL_G",
+   "type": "integer",
+   "default": null,
+   "unit": "No unit",
+   "example": "SQ_NPL_G: 24",
+   "description": "Degree of polynomial for Gauss Quadrature.",
+   "remark": "SQ_NPL_G must be specified if SQ is turned on.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Degree of polynomial for Gauss Quadrature.",
+   "remark_raw": "\\texttt{SQ\\_NPL\\_G} must be specified if SQ is turned on.",
+   "category": "spectral quadrature"
+  },
+  "SQ_RCUT": {
+   "symbol": "{SQ_RCUT",
+   "label": "SQ_RCUT",
+   "type": "double",
+   "default": null,
+   "unit": "Bohr",
+   "example": "SQ_RCUT: 2.0",
+   "description": "Truncation or localization radius",
+   "remark": "SQ_RCUT must be specified if SQ is turned on.",
+   "allow_bool_input": false,
+   "default_remark": "None",
+   "description_raw": "Truncation or localization radius",
+   "remark_raw": "\\texttt{{SQ\\_RCUT}} must be specified if SQ is turned on.",
+   "category": "spectral quadrature"
+  },
+  "SQ_TOL_OCC": {
+   "symbol": "SQ_TOL_OCC",
+   "label": "SQ_TOL_OCC",
+   "type": "double",
+   "default": null,
+   "unit": "No unit",
+   "example": "SQ_TOL_OCC: 1E-5",
+   "description": "Tolerance for occupation corresponding to maximum eigenvalue.",
+   "allow_bool_input": false,
+   "default_remark": "$10^{-6}$",
+   "description_raw": "Tolerance for occupation corresponding to maximum eigenvalue.",
+   "remark_raw": "",
+   "remark": "",
+   "category": "spectral quadrature"
   }
  },
  "other_parameters": {},
  "data_types": [
+  "integer array",
   "string",
-  "other",
   "double",
+  "other",
   "integer",
-  "integer array",
   "double array"
  ]
-}
\ No newline at end of file
+}
diff --git a/sparc/sparc_parsers/aimd.py b/sparc/sparc_parsers/aimd.py
index 697daef9..5dd98a7a 100644
--- a/sparc/sparc_parsers/aimd.py
+++ b/sparc/sparc_parsers/aimd.py
@@ -41,9 +41,7 @@ def _read_aimd(fileobj):
 
     # find the index for all atom type lines. They should be at the
     # top of their block
-    step_bounds = [i for i, x in enumerate(data) if ":MDSTEP:" in x] + [
-        len(data)
-    ]
+    step_bounds = [i for i, x in enumerate(data) if ":MDSTEP:" in x] + [len(data)]
     raw_aimd_blocks = [
         data[start:end] for start, end in zip(step_bounds[:-1], step_bounds[1:])
     ]
@@ -157,10 +155,7 @@ def _read_aimd_step(raw_aimd_text):
             warn(f"MD output keyword {header_name} will not be parsed.")
             value = None
         else:
-            warn(
-                f"MD output keyword {header_name} not known to SPARC. "
-                "Ignore."
-            )
+            warn(f"MD output keyword {header_name} not known to SPARC. " "Ignore.")
             value = None
         if value is not None:
             data[name] = value
@@ -173,6 +168,4 @@ def _write_aimd(
     fileobj,
     data_dict,
 ):
-    raise NotImplementedError(
-        "Writing aimd file from python-api " "not supported!"
-    )
+    raise NotImplementedError("Writing aimd file from python-api " "not supported!")
diff --git a/sparc/sparc_parsers/atoms.py b/sparc/sparc_parsers/atoms.py
index 6266cef5..bda6cc57 100644
--- a/sparc/sparc_parsers/atoms.py
+++ b/sparc/sparc_parsers/atoms.py
@@ -2,21 +2,20 @@
 and vice versa
 """
 
-import numpy as np
+from copy import deepcopy
+from warnings import warn
 
-from ase import Atoms, Atom
+import numpy as np
+from ase import Atom, Atoms
+from ase.constraints import FixAtoms, FixedLine, FixedPlane
 from ase.units import Bohr
 
-# from .sparc_parsers.ion import read_ion, write_ion
-
-from .ion import _ion_coord_to_ase_pos
 from .inpt import _inpt_cell_to_ase_cell
+from .ion import _ion_coord_to_ase_pos
 from .pseudopotential import find_pseudo_path
 from .utils import make_reverse_mapping
-from ase.constraints import FixAtoms, FixedLine, FixedPlane
 
-from warnings import warn
-from copy import deepcopy
+# from .sparc_parsers.ion import read_ion, write_ion
 
 
 def atoms_to_dict(
@@ -117,6 +116,11 @@ def atoms_to_dict(
     cell_au = atoms.cell / Bohr
     inpt_blocks = {"LATVEC": cell_au, "LATVEC_SCALE": [1.0, 1.0, 1.0]}
 
+    # Step 5: retrieve boundary condition information
+    # TODO: have to use space to join the single keywords
+    # breakpoint()
+    inpt_blocks.update(atoms_bc_to_sparc(atoms))
+
     if not isinstance(comments, list):
         comments = comments.split("\n")
     ion_data = {
@@ -163,9 +167,7 @@ def dict_to_atoms(data_dict):
         atoms_count += len(positions)
 
     if "sorting" in data_dict["ion"]:
-        resort = data_dict["ion"]["sorting"].get(
-            "resort", np.arange(len(atoms))
-        )
+        resort = data_dict["ion"]["sorting"].get("resort", np.arange(len(atoms)))
         # Resort may be None
         if len(resort) == 0:
             resort = np.arange(len(atoms))
@@ -178,9 +180,9 @@ def dict_to_atoms(data_dict):
             "Length of resort mapping is different from the number of atoms!"
         )
     # TODO: check if this mapping is correct
-    print(relax_dict)
+    # print(relax_dict)
     sort = make_reverse_mapping(resort)
-    print(resort, sort)
+    # print(resort, sort)
     sorted_relax_dict = {sort[i]: r for i, r in relax_dict.items()}
     # Now we do a sort on the atom indices. The atom positions read from
     # .ion correspond to the `sort` and we use `resort` to transform
@@ -191,8 +193,13 @@ def dict_to_atoms(data_dict):
     constraints = constraints_from_relax(sorted_relax_dict)
     atoms.constraints = constraints
 
-    # TODO: set pbc and relax
-    atoms.pbc = True
+    # @2023.08.31 add support for PBC
+    # TODO: move to a more modular function
+    # TODO: Datatype for BC in the API, should it be string, or string array?
+    sparc_bc = new_data_dict["inpt"]["params"].get("BC", "P P P").split()
+    twist_angle = float(new_data_dict["inpt"]["params"].get("TWIST_ANGLE", 0))
+    modify_atoms_bc(atoms, sparc_bc, twist_angle)
+
     return atoms
 
 
@@ -261,15 +268,11 @@ def constraints_from_relax(relax_dict):
         # DegreeF == 1 --> move along line, fix line
         elif degree_freedom == 1:
             for ind in indices:
-                cons_list.append(
-                    FixedLine(ind, np.array(relax_type).astype(int))
-                )
+                cons_list.append(FixedLine(ind, np.array(relax_type).astype(int)))
         # DegreeF == 1 --> move along line, fix plane
         elif degree_freedom == 2:
             for ind in indices:
-                cons_list.append(
-                    FixedPlane(ind, (~np.array(relax_type)).astype(int))
-                )
+                cons_list.append(FixedPlane(ind, (~np.array(relax_type)).astype(int)))
     return cons_list
 
 
@@ -323,3 +326,71 @@ def relax_from_all_constraints(constraints, natoms):
             # always make it more constrained
             relax[atom_index] = list(np.bitwise_and(relax[atom_index], rdims))
     return relax
+
+
+def modify_atoms_bc(atoms, sparc_bc, twist_angle=0):
+    """Modify the atoms boundary conditions in-place from the bc information
+    sparc_bc is a keyword from inpt
+    twist_angle is the helix twist angle in inpt
+
+    conversion rules:
+    BC: P --> pbc=True
+    BC: D, H, C --> pbc=False
+    """
+    ase_bc = []
+    # print(sparc_bc, type(sparc_bc))
+    for bc_ in sparc_bc:
+        if bc_.upper() in ["C", "H"]:
+            warn(
+                (
+                    "Parsing SPARC's helix or cyclic boundary conditions"
+                    " into ASE atoms is only partially supported. "
+                    "Saving the atoms object into other format may cause "
+                    "data-loss of the SPARC-specific BC information."
+                )
+            )
+            pbc = (
+                False  # Do not confuse ase-gui, we'll manually handle the visualization
+            )
+        elif bc_.upper() == "D":
+            pbc = False
+        elif bc_.upper() == "P":
+            pbc = True
+        else:
+            raise ValueError("Unknown BC keyword values!")
+        ase_bc.append(pbc)
+    atoms.info["sparc_bc"] = [bc_.upper() for bc_ in sparc_bc]
+    if twist_angle != 0:
+        atoms.info["twist_angle (rad/Bohr)"] = twist_angle
+    atoms.pbc = ase_bc
+    return
+
+
+def atoms_bc_to_sparc(atoms):
+    """Use atoms' internal pbc and info to construct inpt blocks
+
+    Returns:
+    a dict containing 'BC' or 'TWIST_ANGLE'
+    """
+    sparc_bc = ["P" if bc_ else "D" for bc_ in atoms.pbc]
+
+    # If "sparc_bc" info is stored in the atoms object, convert again
+    if "sparc_bc" in atoms.info.keys():
+        converted_bc = []
+        stored_sparc_bc = atoms.info["sparc_bc"]
+        for bc1, bc2 in zip(sparc_bc, stored_sparc_bc):
+            # We store helix and cyclic BCs as non-periodic in ase-atoms
+            print(bc1, bc2)
+            if ((bc1 == "D") and (bc2 != "P")) or ((bc1 == "P") and (bc2 == "P")):
+                converted_bc.append(bc2)
+            else:
+                raise ValueError(
+                    "Boundary conditions stored in ASE "
+                    "atoms.pbc and atoms.info['sparc_bc'] "
+                    "are different!"
+                )
+        sparc_bc = converted_bc
+    block = {"BC": " ".join(sparc_bc)}
+    if "twist_angle" in atoms.info.keys():
+        block["TWIST_ANGLE"] = atoms.info["twist_angle (rad/Bohr)"]
+    return block
diff --git a/sparc/sparc_parsers/geopt.py b/sparc/sparc_parsers/geopt.py
index 21182f49..10c9fb38 100644
--- a/sparc/sparc_parsers/geopt.py
+++ b/sparc/sparc_parsers/geopt.py
@@ -40,9 +40,7 @@ def _read_geopt(fileobj):
 
     # find the index for all atom type lines. They should be at the
     # top of their block
-    step_bounds = [i for i, x in enumerate(data) if ":RELAXSTEP:" in x] + [
-        len(data)
-    ]
+    step_bounds = [i for i, x in enumerate(data) if ":RELAXSTEP:" in x] + [len(data)]
     raw_geopt_blocks = [
         data[start:end] for start, end in zip(step_bounds[:-1], step_bounds[1:])
     ]
@@ -59,9 +57,7 @@ def _read_geopt_step(raw_step_text):
     """
     header, body = raw_step_text[0], raw_step_text[1:]
     if ":RELAXSTEP:" not in header:
-        raise ValueError(
-            "Wrong geopt format! The :RELAXSTEP: label is missing."
-        )
+        raise ValueError("Wrong geopt format! The :RELAXSTEP: label is missing.")
     # Geopt file uses 1-indexed step names, convert to 0-indexed
     step = int(header.split(":RELAXSTEP:")[-1]) - 1
     print("Step ", step)
@@ -138,6 +134,4 @@ def _write_geopt(
     fileobj,
     data_dict,
 ):
-    raise NotImplementedError(
-        "Writing geopt file from python-api not supported!"
-    )
+    raise NotImplementedError("Writing geopt file from python-api not supported!")
diff --git a/sparc/sparc_parsers/inpt.py b/sparc/sparc_parsers/inpt.py
index cf6b64ad..196fe448 100644
--- a/sparc/sparc_parsers/inpt.py
+++ b/sparc/sparc_parsers/inpt.py
@@ -32,9 +32,7 @@ def _write_inpt(fileobj, data_dict):
     inpt_dict = data_dict["inpt"]
 
     if "params" not in inpt_dict:
-        raise ValueError(
-            "Input dict for inpt file does not have `params` field!"
-        )
+        raise ValueError("Input dict for inpt file does not have `params` field!")
 
     comments = inpt_dict.get("comments", [])
     banner = "Input File Generated By SPARC ASE Calculator"
@@ -76,27 +74,32 @@ def _inpt_cell_to_ase_cell(data_dict):
     if ("CELL" in inpt_blocks) and ("LATVEC_SCALE" in inpt_blocks):
         # TODO: customize the exception class
         # TODO: how do we convert the rule from doc?
-        raise ValueError(
-            "LATVEC_SCALE and CELL cannot be specified simultaneously!"
-        )
+        raise ValueError("LATVEC_SCALE and CELL cannot be specified simultaneously!")
 
     # if "CELL" in inpt_blocks:
     #     cell = np.eye(inpt_blocks["CELL"]) * Bohr
     if "LATVEC" not in inpt_blocks:
-        raise KeyError("LATVEC is missing in inpt file!")
+        if ("CELL" in inpt_blocks) or ("LATVEC_SCALE" in inpt_blocks):
+            lat_array = np.eye(3) * Bohr
+        else:
+            raise KeyError(
+                "LATVEC is missing in inpt file and no CELL / LATVEC_SCALE provided!"
+            )
+    else:
+        lat_array = np.array(inpt_blocks["LATVEC"]) * Bohr
 
-    lat_array = np.array(inpt_blocks["LATVEC"]) * Bohr
     # LATVEC_SCALE: just multiplies
     if "LATVEC_SCALE" in inpt_blocks:
-        scale = inpt_blocks["LATVEC_SCALE"]
-        # TODO: Is this correct?
+        scale = np.array(inpt_blocks["LATVEC_SCALE"])
         cell = (lat_array.T * scale).T
 
     # CELL: the lengths are in the LATVEC directions
+    # TODO: the documentation about CELL is a bit messy. Is CELL always orthogonal?
+    # Anyway the lat_array when CELL is none should be ok
     elif "CELL" in inpt_blocks:
-        scale = np.array(inpt_blocks["CELL"]) * Bohr
-        unit_lat_array = lat_array / np.linalg.norm(
-            lat_array, axis=1, keepdims=True
+        scale = np.array(inpt_blocks["CELL"])
+        unit_lat_array = (
+            lat_array / np.linalg.norm(lat_array, axis=1, keepdims=True) * Bohr
         )
         cell = (unit_lat_array.T * scale).T
     else:
diff --git a/sparc/sparc_parsers/ion.py b/sparc/sparc_parsers/ion.py
index dbf83433..dd9531de 100644
--- a/sparc/sparc_parsers/ion.py
+++ b/sparc/sparc_parsers/ion.py
@@ -53,9 +53,7 @@ def _read_ion(fileobj):
     sort, resort, new_comments = _read_sort_comment(comments)
 
     # find the index for all atom type lines. They should be at the top of their block
-    atom_type_bounds = [i for i, x in enumerate(data) if "ATOM_TYPE" in x] + [
-        len(data)
-    ]
+    atom_type_bounds = [i for i, x in enumerate(data) if "ATOM_TYPE" in x] + [len(data)]
     atom_blocks = [
         read_block_input(data[start:end], validator=defaultAPI)
         for start, end in zip(atom_type_bounds[:-1], atom_type_bounds[1:])
@@ -131,9 +129,7 @@ def _write_ion(
             val = block.get(key, None)
             # print(key, val)
             if (key not in ["RELAX", "COORD", "COORD_FRAC"]) and (val is None):
-                raise ValueError(
-                    f"Key {key} is not provided! Abort writing ion file"
-                )
+                raise ValueError(f"Key {key} is not provided! Abort writing ion file")
             # TODO: change the API version
             if val is None:
                 continue
@@ -225,7 +221,5 @@ def _read_sort_comment(lines):
             "ASE atoms resort comment block is not properly formatted, this may cause data loss!"
         )
     sort = make_reverse_mapping(resort)
-    assert set(sort) == set(
-        resort
-    ), "Sort and resort info are of different length!"
+    assert set(sort) == set(resort), "Sort and resort info are of different length!"
     return sort, resort, new_lines
diff --git a/sparc/sparc_parsers/out.py b/sparc/sparc_parsers/out.py
index 60002b12..232ac6f1 100644
--- a/sparc/sparc_parsers/out.py
+++ b/sparc/sparc_parsers/out.py
@@ -66,14 +66,10 @@ def _read_sparc_version(header):
     date_str = match[0].strip().replace(",", " ")
     # Accept both abbreviate and full month name
     try:
-        date_version = datetime.strptime(date_str, "%B %d %Y").strftime(
-            "%Y.%m.%d"
-        )
+        date_version = datetime.strptime(date_str, "%B %d %Y").strftime("%Y.%m.%d")
     except ValueError:
         try:
-            date_version = datetime.strptime(date_str, "%b %d %Y").strftime(
-                "%Y.%m.%d"
-            )
+            date_version = datetime.strptime(date_str, "%b %d %Y").strftime("%Y.%m.%d")
         except ValueError:
             warn("Cannot fetch SPARC version information!")
             date_version = None
@@ -129,9 +125,8 @@ def _read_scfs(contents):
 
 
     """
-    convergence_info = _get_block_text(
-        contents, r"Self Consistent Field \(SCF.*?\)"
-    )
+    # import pdb; pdb.set_trace()
+    convergence_info = _get_block_text(contents, r"Self Consistent Field \(SCF.*?\)")
     results_info = _get_block_text(contents, "Energy and force calculation")
 
     if len(convergence_info) != len(results_info):
@@ -146,25 +141,49 @@ def _read_scfs(contents):
         conv, res = step
         # TODO: add support for convergence fields
         conv_lines = conv.splitlines()
+        # conv_header is normally 4-column table
         conv_header = re.split(r"\s{3,}", conv_lines[0])
-        # In some cases the ionic step ends with a warning message
-        # To be flexible, we only extract lines starting with a number
-        conv_array = np.genfromtxt(
-            [l for l in conv_lines if l.split()[0].isdigit()], dtype=float
-        )
-        # TODO: the meaning of the header should me split to the width
 
-        conv_dict = {}
-        for i, field in enumerate(conv_header):
-            field = field.split("(")[0].strip().lower()
-            value = conv_array[:, i]
-            if "free energy" in field:
-                value *= Hartree
-            conv_dict[field] = value
+        scf_sub_steps = []
+        # For ground-state calculations, the output will be only 1 block
+        # For hybrid (HSE/PBE0) calculations the EXX loops will also be included
+        # General rule: we search for the line "Total number of SCF: N", read back N(+1) lines
+        for lino, line in enumerate(conv_lines):
+            if "Total number of SCF:" not in line:
+                continue
+            # import pdb; pdb.set_trace()
+            scf_num = int(line.split(":")[-1])
+            conv_array = np.genfromtxt(
+                [
+                    l
+                    for l in conv_lines[lino - scf_num : lino]
+                    if l.split()[0].isdigit()
+                ],
+                dtype=float,
+                ndmin=2,
+            )
+            conv_dict = {}
+            for i, field in enumerate(conv_header):
+                field = field.strip()
+                value = conv_array[:, i]
+                # TODO: re-use the value conversion function in res part
+                if "Ha/atom" in field:
+                    value *= Hartree
+                    field.replace("Ha/atom", "eV/atom")
+                if "Iteration" in field:
+                    value = value.astype(int)
+                conv_dict[field] = value
+            # Determine if the current block is a ground-state or EXX
+            name_line = conv_lines[lino - scf_num - 1]
+            if "Iteration" in name_line:
+                name = "ground state"
+            else:
+                name = name_line
 
-        current_step["convergence"] = conv_dict
-        # TODO: what is this?
-        {"header": conv_header, "values": conv_array}
+            conv_dict["name"] = name
+            scf_sub_steps.append(conv_dict)
+
+        current_step["convergence"] = scf_sub_steps
 
         res = res.splitlines()
         for line in res:
@@ -227,6 +246,4 @@ def _write_out(
     fileobj,
     data_dict,
 ):
-    raise NotImplementedError(
-        "Writing output file from python-api not supported!"
-    )
+    raise NotImplementedError("Writing output file from python-api not supported!")
diff --git a/sparc/sparc_parsers/pseudopotential.py b/sparc/sparc_parsers/pseudopotential.py
index a8c00c68..095f07ea 100644
--- a/sparc/sparc_parsers/pseudopotential.py
+++ b/sparc/sparc_parsers/pseudopotential.py
@@ -46,9 +46,7 @@ def parse_psp8_header(text):
     mgroup = match.groupdict()
     psp8_data["symbol"] = mgroup["symbol"].strip()
     psp8_data["psp8ver"] = mgroup["psp8ver"].strip()
-    psp8_data["r_core"] = np.fromstring(
-        mgroup["r_core"].strip(), sep=" ", dtype=float
-    )
+    psp8_data["r_core"] = np.fromstring(mgroup["r_core"].strip(), sep=" ", dtype=float)
     # Line 2
     zatom, zion, pspd, *_ = header[1].split()
     psp8_data["zatom"] = float(zatom)
@@ -173,9 +171,7 @@ def find_pseudo_path(symbol, search_path=None, pseudopotential_mapping={}):
         str_psp = str(mapping_psp)
         mapping_psp = Path(mapping_psp)
         # if psp contains any path information (/, \\), treat is as a direct file
-        is_node_file_name = (mapping_psp.name == str_psp) and (
-            os.sep not in str_psp
-        )
+        is_node_file_name = (mapping_psp.name == str_psp) and (os.sep not in str_psp)
         if is_node_file_name:
             if search_path is None:
                 raise NoMatchingPseudopotential(
diff --git a/sparc/sparc_parsers/static.py b/sparc/sparc_parsers/static.py
index 6dcbbb96..1e2c227f 100644
--- a/sparc/sparc_parsers/static.py
+++ b/sparc/sparc_parsers/static.py
@@ -39,8 +39,7 @@ def _read_static(fileobj):
     block_bounds = [i for i, x in enumerate(data) if ":" in x] + [len(data)]
     # blocks = [read_static_block(data[start:end]) for start, end in zip(block_bounds[:-1], block_bounds[1:])]
     raw_blocks = [
-        data[start:end]
-        for start, end in zip(block_bounds[:-1], block_bounds[1:])
+        data[start:end] for start, end in zip(block_bounds[:-1], block_bounds[1:])
     ]
     static_dict = read_static_blocks(raw_blocks)
     return {"static": static_dict}
@@ -75,8 +74,14 @@ def _read_static_block(raw_block):
         name = "free energy"
     elif "Atomic forces" in header_name:
         name = "forces"
-    elif "Stress" in header_name:
+    elif "Stress (GPa)" in header_name:
         name = "stress"
+    elif "Stress equiv." in header_name:
+        name = "stress_equiv"
+    elif "Stress (Ha/Bohr)" in header_name:
+        name = "stress_1d"
+    elif "Stress (Ha/Bohr**2)" in header_name:
+        name = "stress_2d"
     elif "Fractional coordinates" in header_name:
         # Fractional coordinates of Si -- > name=coord_frac symbol="Si"
         name = "coord_frac"
@@ -126,6 +131,7 @@ def read_static_blocks(raw_blocks):
             value = raw_value * Hartree / Bohr
         elif name == "stress":
             # Stress is in eV/Ang^3, may need to convert to Virial later when cell is known
+            # For low-dimension stress info, use stress_equiv
             stress_ev_a3 = raw_value * GPa
             if stress_ev_a3.shape != (3, 3):
                 raise ValueError("Stress from static file is not a 3x3 matrix!")
@@ -141,6 +147,14 @@ def read_static_blocks(raw_blocks):
                     stress_ev_a3[0, 1],
                 ]
             )
+        elif name == "stress_equiv":
+            # Only store the size up to the max. periodic directions,
+            # let the atom parser decide how to transform the matrix
+            value = raw_value * GPa
+        elif name == "stress_1d":
+            value = raw_value * Hartree / Bohr
+        elif name == "stress_2d":
+            value = raw_value * Hartree / (Bohr**2)
 
         # Non-frac coord
         if value is not None:
@@ -176,6 +190,4 @@ def _write_static(
     fileobj,
     data_dict,
 ):
-    raise NotImplementedError(
-        "Writing static file from python-api not supported!"
-    )
+    raise NotImplementedError("Writing static file from python-api not supported!")
diff --git a/sparc/sparc_parsers/utils.py b/sparc/sparc_parsers/utils.py
index cd8b976a..6ee7d5ba 100644
--- a/sparc/sparc_parsers/utils.py
+++ b/sparc/sparc_parsers/utils.py
@@ -29,27 +29,44 @@ def bisect_and_strip(text, delimiter):
 
 
 def read_block_input(block, validator=None):
-    """Read blocks of inputs from ion or inpt file and convert with validator"""
+    """Read blocks of inputs from ion or inpt file and convert with validator
+
+    the following inputs are accepted:
+    1) single line input: KEY: VALUE
+    2) multiline input: KEY: VALUE1 \n VALUE2 --> (concanate the values)
+    3) multiline input w/ blank first line: KEY: \n VALUE1 \n VALUE2 --> (append the values)
+    """
     block_dict = {}
     multiline_key = ""
+    concat = False
     use_validator = True if validator else False
     for line in block:
         if ":" not in line:
-            # no key, assume multiline value
-            block_dict[multiline_key].append(line.strip())
+            # import pdb; pdb.set_trace()
+            # no key, assume multiline value.
+            # be careful not to add blank lines
+            if multiline_key:
+                if concat:
+                    block_dict[multiline_key] = (
+                        block_dict[multiline_key] + f" {line.strip()}"
+                    )
+                else:
+                    block_dict[multiline_key].append(line.strip())
             continue
         key, value = bisect_and_strip(line, ":")
         key = key.upper()
-        # print(key, value)
+
         if key and value:
             block_dict[key] = value
+            multiline_key = key
+            concat = True
         elif key:
             # no value, assume that this key has a list of values
             # in the following lines
             block_dict[key] = []
             multiline_key = key
+            concat = False
     for key, val in block_dict.items():
-        # print(key, val)
         _use_validator_this_key = use_validator
         if _use_validator_this_key:
             if key not in validator.parameters.keys():
diff --git a/sparc/utils.py b/sparc/utils.py
index a6e88571..71aa2905 100644
--- a/sparc/utils.py
+++ b/sparc/utils.py
@@ -3,9 +3,26 @@
 import os
 import shutil
 import numpy as np
+from warnings import warn
 from typing import Union, List, Optional
 
 
+def deprecated(message):
+    def decorator(func):
+        def new_func(*args, **kwargs):
+            warn(
+                "Function {} is deprecated sparc-dft-api >= v0.2! {}".format(
+                    func.__name__, message
+                ),
+                category=DeprecationWarning,
+            )
+            return func(*args, **kwargs)
+
+        return new_func
+
+    return decorator
+
+
 def string2index(string: str) -> Union[int, slice, str]:
     """Convert index string to either int or slice
     This method is a copy of ase.io.formats.string2index


Test patch:
diff --git a/.github/workflows/installation_test.yml b/.github/workflows/installation_test.yml
index ea7f6a85..0a65c243 100644
--- a/.github/workflows/installation_test.yml
+++ b/.github/workflows/installation_test.yml
@@ -30,22 +30,29 @@ jobs:
     - name: Install dependencies
       run: |
         # mamba install -c conda-forge ase>=3.22 pymatgen flake8 pytest
-        mamba install -c alchem0x2a sparc
+        mamba install -c conda-forge sparc-x
         # pip install pyfakefs
     - name: Install package
       run: |
         pip install -e ".[test]"
         # Download the external psp data
         python -m sparc.download_data
+    - name: Download SPARC output files to SPARC-master
+      run: |
+        wget https://github.com/SPARC-X/SPARC/archive/refs/heads/master.zip
+        unzip master.zip
     - name: Test with pytest
       run: |
         # python -m pytest -svv tests/ --cov=sparc --cov-report=json --cov-report=html
+        export SPARC_TESTS_DIR="./SPARC-master/tests"
+        export ASE_SPARC_COMMAND="mpirun -n 1 sparc"
         coverage run -a -m pytest -svv tests/
         coverage json --omit="tests/*.py"
         coverage html --omit="tests/*.py"
         COVERAGE=`cat coverage.json | jq .totals.percent_covered | xargs printf '%.*f' 0`
         echo "Current coverage is $COVERAGE"
         echo "COVPERCENT=$COVERAGE" >> $GITHUB_ENV
+        
     - name: Lint with flake8
       run: |
         echo $CONDA_PREFIX
@@ -58,8 +65,8 @@ jobs:
         # A coverage test badge
         anybadge --value=$COVPERCENT --file=badges/coverage.svg coverage
         # A version badge
-        PIPVER=`pip show sparc-dft-api | grep Version | cut -d ' ' -f2`
-        anybadge --value=$PIPVER --file=badges/package.svg -l sparc-dft-api
+        PIPVER=`pip show sparc-x-api | grep Version | cut -d ' ' -f2`
+        anybadge --value=$PIPVER --file=badges/package.svg -l sparc-x-api
     - name: Manually add git badges
       run: |
         # Assuming a badges branch already exists!
diff --git a/sparc/quicktest.py b/sparc/quicktest.py
index fe3423f6..e23fdc7d 100644
--- a/sparc/quicktest.py
+++ b/sparc/quicktest.py
@@ -120,9 +120,7 @@ def main():
     results["Pseudopotential"] = psp_test()
     results["JSON API"] = api_test()
     results["SPARC command"] = command_test()
-    results["Calculation"] = (
-        False if results["SPARC command"] is False else calc_test()
-    )
+    results["Calculation"] = False if results["SPARC command"] is False else calc_test()
 
     cprint(
         "\nSummary of test results",
diff --git a/tests/archive/test.py b/tests/archive/test.py
index 72518743..e7573101 100644
--- a/tests/archive/test.py
+++ b/tests/archive/test.py
@@ -16,9 +16,7 @@
 
 print("input writing functional")
 # check reading and writing .ion files
-atoms.set_constraint(
-    [FixAtoms([0]), FixedLine(1, [0, 1, 0]), FixedPlane(2, [1, 0, 0])]
-)
+atoms.set_constraint([FixAtoms([0]), FixedLine(1, [0, 1, 0]), FixedPlane(2, [1, 0, 0])])
 write_ion(open("in1.ion", "w"), atoms, pseudo_dir=".")
 recovered_atoms = read_ion(open("in1.ion", "r"))
 assert compare_atoms(atoms, recovered_atoms) == []
diff --git a/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.inpt b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.inpt
new file mode 100644
index 00000000..f0235cc5
--- /dev/null
+++ b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.inpt
@@ -0,0 +1,21 @@
+# nprocs: 1
+# Test: H2O sheet in x,y plane #
+LATVEC_SCALE: 24 5.67 5.67 
+LATVEC:
+1 0 0
+0 1 0
+0 0 1
+MESH_SPACING: 0.2
+FD_ORDER: 12
+BC: D P P
+EXCHANGE_CORRELATION: GGA_PBE
+ELEC_TEMP_TYPE: fd
+SMEARING: 0.003674932
+TOL_SCF: 1e-6
+MIXING_VARIABLE: potential
+MIXING_PRECOND: none
+NSTATES: 8
+PRINT_FORCES: 1
+PRINT_ATOMS: 1
+
+CALC_STRESS: 1
diff --git a/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.ion b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.ion
new file mode 100644
index 00000000..0eff5d88
--- /dev/null
+++ b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.ion
@@ -0,0 +1,29 @@
+#=========================
+# format of ion file
+#=========================
+# ATOM_TYPE: <atom type name> 
+# PSEUDO_POT: <path/to/pseudopotential/file>
+# N_TYPE_ATOM: <num of atoms of this type>
+# COORD:
+# <xcoord> <ycoord> <zcoord>
+# ...
+# RELAX:
+# <xrelax> <yrelax> <zrelax>
+# ...
+
+# Reminder: when changing number of atoms, change the RELAX flags accordingly
+#           as well.
+
+ATOM_TYPE: H                 # atom type
+PSEUDO_POT: ../../../psps/01_H_1_1.0_1.0_pbe_v1.0.psp8
+N_TYPE_ATOM: 2               # number of atoms of this type
+COORD:                       # coordinates follows
+   12.556010134214192  4.265111864110533   2.834589187490385  
+   12.556010134214192  1.404066510870237   2.834589187490385  
+
+ATOM_TYPE: O               # atom type 
+PSEUDO_POT: ../../../psps/08_O_6_1.2_1.4_pbe_n_v1.0.psp8
+N_TYPE_ATOM: 1               # number of atoms of this type
+COORD:                       # coordinates follows
+   11.448630624967949  2.834589187490385   2.834589187490385   
+
diff --git a/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.out b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.out
new file mode 100644
index 00000000..846fadb1
--- /dev/null
+++ b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.out
@@ -0,0 +1,138 @@
+***************************************************************************
+*                       SPARC (version Feb 03, 2023)                      *
+*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *
+*           Distributed under GNU General Public License 3 (GPL)          *
+*                   Start time: Sun Feb  5 13:17:48 2023                  *
+***************************************************************************
+                           Input parameters                                
+***************************************************************************
+LATVEC_SCALE: 24 5.67 5.67 
+LATVEC:
+1 0 0 
+0 1 0 
+0 0 1 
+FD_GRID: 120 29 29
+FD_ORDER: 12
+BC: D P P
+KPOINT_GRID: 1 1 1
+KPOINT_SHIFT: 0 0 0
+SPIN_TYP: 0
+ELEC_TEMP_TYPE: Fermi-Dirac
+SMEARING: 0.003674932
+EXCHANGE_CORRELATION: GGA_PBE
+NSTATES: 8
+CHEB_DEGREE: 35
+CHEFSI_BOUND_FLAG: 0
+CALC_STRESS: 1
+MAXIT_SCF: 100
+MINIT_SCF: 2
+MAXIT_POISSON: 3000
+TOL_SCF: 1.00E-06
+POISSON_SOLVER: AAR
+TOL_POISSON: 1.00E-08
+TOL_LANCZOS: 1.00E-02
+TOL_PSEUDOCHARGE: 1.00E-09
+MIXING_VARIABLE: potential
+MIXING_PRECOND: none
+MIXING_PARAMETER: 0.3
+MIXING_HISTORY: 7
+PULAY_FREQUENCY: 1
+PULAY_RESTART: 0
+REFERENCE_CUTOFF: 0.5
+RHO_TRIGGER: 4
+FIX_RAND: 0
+VERBOSITY: 1
+PRINT_FORCES: 1
+PRINT_ATOMS: 1
+PRINT_EIGEN: 0
+PRINT_DENSITY: 0
+PRINT_ENERGY_DENSITY: 0
+OUTPUT_FILE: H2O_sheet/temp_run3/H2O_sheet
+***************************************************************************
+                                Cell                                       
+***************************************************************************
+Lattice vectors (Bohr):
+24.000000000000000 0.000000000000000 0.000000000000000 
+0.000000000000000 5.670000000000000 0.000000000000000 
+0.000000000000000 0.000000000000000 5.670000000000000 
+Volume: 7.7157360000E+02 (Bohr^3)
+***************************************************************************
+                           Parallelization                                 
+***************************************************************************
+NP_SPIN_PARAL: 1
+NP_KPOINT_PARAL: 1
+NP_BAND_PARAL: 8
+NP_DOMAIN_PARAL: 6 1 1
+NP_DOMAIN_PHI_PARAL: 12 2 2
+EIG_SERIAL_MAXNS: 1500
+***************************************************************************
+                             Initialization                                
+***************************************************************************
+Number of processors               :  48
+Mesh spacing in x-direction        :  0.2 (Bohr)
+Mesh spacing in y-direction        :  0.195517 (Bohr)
+Mesh spacing in z-direction        :  0.195517 (Bohr)
+Number of symmetry adapted k-points:  1
+Output printed to                  :  H2O_sheet/temp_run3/H2O_sheet.out
+Total number of atom types         :  2
+Total number of atoms              :  3
+Total number of electrons          :  8
+Atom type 1  (valence electrons)   :  H 1
+Pseudopotential                    :  ../psps/01_H_1_1.0_1.0_pbe_v1.0.psp8
+Pseudocharge radii of atom type 1  :  4.20 4.30 4.30 (x, y, z dir)
+Number of atoms of type 1          :  2
+Atom type 2  (valence electrons)   :  O 6
+Pseudopotential                    :  ../psps/08_O_6_1.2_1.4_pbe_n_v1.0.psp8
+Pseudocharge radii of atom type 2  :  7.20 7.23 7.23 (x, y, z dir)
+Number of atoms of type 2          :  1
+Estimated total memory usage       :  59.00 MB
+Estimated memory per processor     :  1.23 MB
+===================================================================
+                    Self Consistent Field (SCF#1)                     
+===================================================================
+Iteration     Free Energy (Ha/atom)   SCF Error        Timing (sec)
+1            -5.8502567922E+00        4.363E-01        0.175
+2            -5.8736011131E+00        8.444E-02        0.081
+3            -5.8741485773E+00        4.037E-02        0.090
+4            -5.8742064274E+00        2.879E-02        0.076
+5            -5.8742349219E+00        9.541E-03        0.074
+6            -5.8742365448E+00        2.953E-03        0.064
+7            -5.8742369371E+00        1.289E-03        0.066
+8            -5.8742369712E+00        3.459E-04        0.058
+9            -5.8742369911E+00        1.962E-04        0.056
+10           -5.8742369936E+00        7.025E-05        0.039
+11           -5.8742369876E+00        9.685E-06        0.034
+12           -5.8742369950E+00        2.589E-06        0.031
+13           -5.8742369979E+00        7.358E-07        0.023
+Total number of SCF: 13    
+====================================================================
+                    Energy and force calculation                    
+====================================================================
+Free energy per atom               : -5.8742369979E+00 (Ha/atom)
+Total free energy                  : -1.7622710994E+01 (Ha)
+Band structure energy              : -4.0591778476E+00 (Ha)
+Exchange correlation energy        : -4.9328873198E+00 (Ha)
+Self and correction energy         : -2.6913043403E+01 (Ha)
+-Entropy*kb*T                      : -1.7023417939E-07 (Ha)
+Fermi level                        : -1.8886200437E-01 (Ha)
+RMS force                          :  2.4015614396E-02 (Ha/Bohr)
+Maximum force                      :  2.9014859269E-02 (Ha/Bohr)
+Time for force calculation         :  0.017 (sec)
+Maximum stress                     :  5.9325459799E-03 (Ha/Bohr**2)
+Maximum stress equiv. to periodic  :  7.2725636831E+00 (GPa)
+Time for stress calculation        :  0.031 (sec)
+***************************************************************************
+                               Timing info                                 
+***************************************************************************
+Total walltime                     :  1.390 sec
+___________________________________________________________________________
+
+***************************************************************************
+*             Material Physics & Mechanics Group, Georgia Tech            *
+*                       PI: Phanish Suryanarayana                         *
+*               List of contributors: See the documentation               *
+*         Citation: See README.md or the documentation for details        *
+*  Acknowledgements: U.S. DOE SC (DE-SC0019410), U.S. DOE NNSA (ASC)      *
+*      {Preliminary developments: U.S. NSF (1333500,1663244,1553212)}     *
+***************************************************************************
+                                                                           
diff --git a/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.static b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.static
new file mode 100644
index 00000000..665697cc
--- /dev/null
+++ b/tests/outputs/H2O_sheet_yz.sparc/H2O_sheet.static
@@ -0,0 +1,19 @@
+***************************************************************************
+                            Atom positions                                 
+***************************************************************************
+Fractional coordinates of H:
+      0.5231670889       0.7522243147       0.4999275463
+      0.5231670889       0.2476307779       0.4999275463
+Fractional coordinates of O:
+      0.4770262760       0.4999275463       0.4999275463
+Total free energy (Ha): -1.762271099378250E+01
+Atomic forces (Ha/Bohr):
+ -7.0091512819E-03  -2.8155137475E-02   4.7625587321E-08
+ -7.0083540340E-03   2.8155728229E-02   4.1091979535E-08
+  1.4017505316E-02  -5.9075452258E-07  -8.8717566856E-08
+Stress (Ha/Bohr**2): 
+ -5.9325459799E-03   3.7871815297E-11 
+  3.7871815297E-11  -5.4929931315E-03 
+Stress equiv. to all periodic (GPa): 
+ -7.2725636831E+00   4.6426136346E-08 
+  4.6426136346E-08  -6.7337265476E+00 
diff --git a/tests/outputs/H2O_sheet_yz.sparc/output.sparc b/tests/outputs/H2O_sheet_yz.sparc/output.sparc
new file mode 100644
index 00000000..e78d5d26
--- /dev/null
+++ b/tests/outputs/H2O_sheet_yz.sparc/output.sparc
@@ -0,0 +1,21 @@
+---------------------------------------
+Begin PBS Prologue Wed Oct 27 17:53:39 EDT 2021
+Job ID:     7040483.sched-hive.pace.gatech.edu
+User ID:    skumar416
+Job name:   testing_suite
+Queue:      hive-interact
+End PBS Prologue Wed Oct 27 17:53:39 EDT 2021
+---------------------------------------
+/storage/hive/scratch/5/skumar416/00merge_dev_sparc_10272021/sparc_public/tests_LHQ/H2O_sheet/temp_run3
+---------------------------------------
+Begin PBS Epilogue Wed Oct 27 17:54:05 EDT 2021
+Job ID:     7040483.sched-hive.pace.gatech.edu
+User ID:    skumar416
+Job name:   testing_suite
+Resources:  nodes=1:ppn=1,mem=10gb,walltime=01:00:00,neednodes=1:ppn=1
+Rsrc Used:  cput=00:00:26,vmem=0kb,walltime=00:00:26,mem=0kb,energy_used=0
+Queue:      hive-interact
+Nodes:     
+atl1-1-01-017-3-l.pace.gatech.edu
+End PBS Epilogue Wed Oct 27 17:54:05 EDT 2021
+---------------------------------------
diff --git a/tests/outputs/H2O_wire_z.sparc/H2O_wire.inpt b/tests/outputs/H2O_wire_z.sparc/H2O_wire.inpt
new file mode 100644
index 00000000..b914ee64
--- /dev/null
+++ b/tests/outputs/H2O_wire_z.sparc/H2O_wire.inpt
@@ -0,0 +1,21 @@
+# nprocs: 1
+# Test: H2O wire in z direction #
+LATVEC_SCALE: 14.00 14.00 6.00
+LATVEC:
+1 0 0
+0 1 0
+0 0 1
+MESH_SPACING: 0.1
+FD_ORDER: 12
+BC: D D P
+#BOUNDARY_CONDITION: 4
+EXCHANGE_CORRELATION: GGA_PBE
+ELEC_TEMP_TYPE: fd
+SMEARING: 0.003674932
+TOL_SCF: 1e-6
+
+#MIXING_VARIABLE: density
+PRINT_FORCES: 1
+PRINT_ATOMS: 1
+CALC_STRESS: 1
+
diff --git a/tests/outputs/H2O_wire_z.sparc/H2O_wire.ion b/tests/outputs/H2O_wire_z.sparc/H2O_wire.ion
new file mode 100644
index 00000000..25a57558
--- /dev/null
+++ b/tests/outputs/H2O_wire_z.sparc/H2O_wire.ion
@@ -0,0 +1,29 @@
+#=========================
+# format of ion file
+#=========================
+# ATOM_TYPE: <atom type name> 
+# PSEUDO_POT: <path/to/pseudopotential/file>
+# N_TYPE_ATOM: <num of atoms of this type>
+# COORD:
+# <xcoord> <ycoord> <zcoord>
+# ...
+# RELAX:
+# <xrelax> <yrelax> <zrelax>
+# ...
+
+# Reminder: when changing number of atoms, change the RELAX flags accordingly
+#           as well.
+
+ATOM_TYPE: H                 # atom type
+PSEUDO_POT: ../../../psps/01_H_1_1.0_1.0_pbe_v1.0.psp8
+N_TYPE_ATOM: 2               # number of atoms of this type
+COORD:                       # coordinates follows
+   7.369126503082081   7.000000000000000   4.430522676620148
+   7.369126503082081   7.000000000000000   1.569477323379852
+
+ATOM_TYPE: O               # atom type 
+PSEUDO_POT: ../../../psps/08_O_6_1.2_1.4_pbe_n_v1.0.psp8
+N_TYPE_ATOM: 1               # number of atoms of this type
+COORD:                       # coordinates follows
+   6.261746993835837   7.000000000000000   3.000000000000000
+
diff --git a/tests/outputs/H2O_wire_z.sparc/H2O_wire.out b/tests/outputs/H2O_wire_z.sparc/H2O_wire.out
new file mode 100644
index 00000000..88039a0e
--- /dev/null
+++ b/tests/outputs/H2O_wire_z.sparc/H2O_wire.out
@@ -0,0 +1,145 @@
+***************************************************************************
+*                       SPARC (version Feb 03, 2023)                      *
+*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *
+*           Distributed under GNU General Public License 3 (GPL)          *
+*                   Start time: Sun Feb  5 13:48:03 2023                  *
+***************************************************************************
+                           Input parameters                                
+***************************************************************************
+LATVEC_SCALE: 14 14 6 
+LATVEC:
+1 0 0 
+0 1 0 
+0 0 1 
+FD_GRID: 140 140 60
+FD_ORDER: 12
+BC: D D P
+KPOINT_GRID: 1 1 1
+KPOINT_SHIFT: 0 0 0
+SPIN_TYP: 0
+ELEC_TEMP_TYPE: Fermi-Dirac
+SMEARING: 0.003674932
+EXCHANGE_CORRELATION: GGA_PBE
+NSTATES: 9
+CHEB_DEGREE: 50
+CHEFSI_BOUND_FLAG: 0
+CALC_STRESS: 1
+MAXIT_SCF: 100
+MINIT_SCF: 2
+MAXIT_POISSON: 3000
+TOL_SCF: 1.00E-06
+POISSON_SOLVER: AAR
+TOL_POISSON: 1.00E-08
+TOL_LANCZOS: 1.00E-02
+TOL_PSEUDOCHARGE: 1.00E-09
+MIXING_VARIABLE: density
+MIXING_PRECOND: kerker
+TOL_PRECOND: 1.00E-05
+PRECOND_KERKER_KTF: 1
+PRECOND_KERKER_THRESH: 0.1
+MIXING_PARAMETER: 0.3
+MIXING_HISTORY: 7
+PULAY_FREQUENCY: 1
+PULAY_RESTART: 0
+REFERENCE_CUTOFF: 0.5
+RHO_TRIGGER: 4
+FIX_RAND: 0
+VERBOSITY: 1
+PRINT_FORCES: 1
+PRINT_ATOMS: 1
+PRINT_EIGEN: 0
+PRINT_DENSITY: 0
+PRINT_ENERGY_DENSITY: 0
+OUTPUT_FILE: H2O_wire/temp_run1/H2O_wire
+***************************************************************************
+                                Cell                                       
+***************************************************************************
+Lattice vectors (Bohr):
+14.000000000000000 0.000000000000000 0.000000000000000 
+0.000000000000000 14.000000000000000 0.000000000000000 
+0.000000000000000 0.000000000000000 6.000000000000000 
+Volume: 1.1760000000E+03 (Bohr^3)
+***************************************************************************
+                           Parallelization                                 
+***************************************************************************
+NP_SPIN_PARAL: 1
+NP_KPOINT_PARAL: 1
+NP_BAND_PARAL: 9
+NP_DOMAIN_PARAL: 2 5 1
+NP_DOMAIN_PHI_PARAL: 6 8 2
+EIG_SERIAL_MAXNS: 1500
+***************************************************************************
+                             Initialization                                
+***************************************************************************
+Number of processors               :  96
+Mesh spacing                       :  0.1 (Bohr)
+Number of symmetry adapted k-points:  1
+Output printed to                  :  H2O_wire/temp_run1/H2O_wire.out
+Total number of atom types         :  2
+Total number of atoms              :  3
+Total number of electrons          :  8
+Atom type 1  (valence electrons)   :  H 1
+Pseudopotential                    :  ../psps/01_H_1_1.0_1.0_pbe_v1.0.psp8
+Pseudocharge radii of atom type 1  :  3.60 3.60 3.60 (x, y, z dir)
+Number of atoms of type 1          :  2
+Atom type 2  (valence electrons)   :  O 6
+Pseudopotential                    :  ../psps/08_O_6_1.2_1.4_pbe_n_v1.0.psp8
+Pseudocharge radii of atom type 2  :  6.60 6.60 6.60 (x, y, z dir)
+Number of atoms of type 2          :  1
+Estimated total memory usage       :  746.27 MB
+Estimated memory per processor     :  7.77 MB
+===================================================================
+                    Self Consistent Field (SCF#1)                     
+===================================================================
+Iteration     Free Energy (Ha/atom)   SCF Error        Timing (sec)
+1            -5.8784017210E+00        1.129E-01        1.492
+2            -5.8982443200E+00        8.419E-02        0.761
+3            -5.8953947597E+00        4.527E-02        0.799
+4            -5.8952382397E+00        3.389E-02        0.829
+5            -5.8950032653E+00        6.974E-03        0.782
+6            -5.8949794718E+00        4.175E-03        0.797
+7            -5.8949968332E+00        1.602E-03        0.735
+8            -5.8950039595E+00        1.414E-03        0.796
+9            -5.8950207782E+00        5.333E-04        0.758
+10           -5.8950288210E+00        3.714E-04        0.763
+11           -5.8950357494E+00        1.305E-04        0.722
+12           -5.8950384443E+00        7.096E-05        0.671
+13           -5.8950395670E+00        3.247E-05        0.611
+14           -5.8950399363E+00        2.173E-05        0.598
+15           -5.8950399645E+00        1.134E-05        0.615
+16           -5.8950400046E+00        6.239E-06        0.541
+17           -5.8950400046E+00        3.638E-06        0.618
+18           -5.8950400152E+00        1.356E-06        0.543
+19           -5.8950400124E+00        8.640E-07        0.555
+Total number of SCF: 19    
+====================================================================
+                    Energy and force calculation                    
+====================================================================
+Free energy per atom               : -5.8950400124E+00 (Ha/atom)
+Total free energy                  : -1.7685120037E+01 (Ha)
+Band structure energy              : -4.2451123247E+00 (Ha)
+Exchange correlation energy        : -4.8893741610E+00 (Ha)
+Self and correction energy         : -2.6913658068E+01 (Ha)
+-Entropy*kb*T                      : -2.2710844514E-10 (Ha)
+Fermi level                        : -2.2199605614E-01 (Ha)
+RMS force                          :  1.9797939601E-02 (Ha/Bohr)
+Maximum force                      :  2.3630179638E-02 (Ha/Bohr)
+Time for force calculation         :  0.025 (sec)
+Maximum stress                     :  1.9995447051E-02 (Ha/Bohr)
+Maximum stress equiv. to periodic  :  3.0014610284E+00 (GPa)
+Time for stress calculation        :  0.044 (sec)
+***************************************************************************
+                               Timing info                                 
+***************************************************************************
+Total walltime                     :  15.440 sec
+___________________________________________________________________________
+
+***************************************************************************
+*             Material Physics & Mechanics Group, Georgia Tech            *
+*                       PI: Phanish Suryanarayana                         *
+*               List of contributors: See the documentation               *
+*         Citation: See README.md or the documentation for details        *
+*  Acknowledgements: U.S. DOE SC (DE-SC0019410), U.S. DOE NNSA (ASC)      *
+*      {Preliminary developments: U.S. NSF (1333500,1663244,1553212)}     *
+***************************************************************************
+                                                                           
diff --git a/tests/outputs/H2O_wire_z.sparc/H2O_wire.static b/tests/outputs/H2O_wire_z.sparc/H2O_wire.static
new file mode 100644
index 00000000..b6d882a3
--- /dev/null
+++ b/tests/outputs/H2O_wire_z.sparc/H2O_wire.static
@@ -0,0 +1,17 @@
+***************************************************************************
+                            Atom positions                                 
+***************************************************************************
+Fractional coordinates of H:
+      0.5263661788       0.5000000000       0.7384204461
+      0.5263661788       0.5000000000       0.2615795539
+Fractional coordinates of O:
+      0.4472676424       0.5000000000       0.5000000000
+Total free energy (Ha): -1.768512003719815E+01
+Atomic forces (Ha/Bohr):
+ -6.0667427045E-03   3.2520689190E-08  -2.2838126515E-02
+ -6.0667226226E-03   4.7399932144E-08   2.2838125847E-02
+  1.2133465327E-02  -7.9920621334E-08   6.6746163366E-10
+Stress (Ha/Bohr): 
+ -1.9995447051E-02 
+Stress equiv. to all periodic (GPa): 
+ -3.0014610284E+00 
diff --git a/tests/outputs/H2O_wire_z.sparc/output.sparc b/tests/outputs/H2O_wire_z.sparc/output.sparc
new file mode 100644
index 00000000..29ebc9d6
--- /dev/null
+++ b/tests/outputs/H2O_wire_z.sparc/output.sparc
@@ -0,0 +1,21 @@
+---------------------------------------
+Begin PBS Prologue Wed Oct 27 17:54:10 EDT 2021
+Job ID:     7040484.sched-hive.pace.gatech.edu
+User ID:    skumar416
+Job name:   testing_suite
+Queue:      hive-interact
+End PBS Prologue Wed Oct 27 17:54:10 EDT 2021
+---------------------------------------
+/storage/hive/scratch/5/skumar416/00merge_dev_sparc_10272021/sparc_public/tests_LHQ/H2O_wire/temp_run1
+---------------------------------------
+Begin PBS Epilogue Wed Oct 27 17:55:00 EDT 2021
+Job ID:     7040484.sched-hive.pace.gatech.edu
+User ID:    skumar416
+Job name:   testing_suite
+Resources:  nodes=1:ppn=1,mem=10gb,walltime=01:00:00,neednodes=1:ppn=1
+Rsrc Used:  cput=00:00:50,vmem=645340kb,walltime=00:00:50,mem=98232kb,energy_used=0
+Queue:      hive-interact
+Nodes:     
+atl1-1-01-017-3-l.pace.gatech.edu
+End PBS Epilogue Wed Oct 27 17:55:00 EDT 2021
+---------------------------------------
diff --git a/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.inpt b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.inpt
new file mode 100644
index 00000000..3c4d3593
--- /dev/null
+++ b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.inpt
@@ -0,0 +1,16 @@
+# nprocs: 24
+CELL: 30.349896582200003 0.299199300341885 3.135906672400000
+TWIST_ANGLE: 0.047705389796071
+FD_GRID: 120 42 25
+KPOINT_GRID: 1 3 4
+KPOINT_SHIFT: 0 0 0
+CHEB_DEGREE: 100
+BC: D C H 
+EXCHANGE_CORRELATION: GGA_PBE
+TOL_SCF: 1e-6
+SMEARING: 0.001
+ELEC_TEMP_TYPE: fd
+CALC_STRESS: 1
+PRINT_FORCES: 1
+MIXING_VARIABLE: potential
+MIXING_PRECOND: none
\ No newline at end of file
diff --git a/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.ion b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.ion
new file mode 100644
index 00000000..bd70e8ed
--- /dev/null
+++ b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.ion
@@ -0,0 +1,11 @@
+ATOM_TYPE: W
+N_TYPE_ATOM: 1
+PSEUDO_POT: ../../../psps/74_W_14_2.3_2.5_pbe_n_v1.0.psp8
+COORD:
+35.932050040279556    5.415857727989695    3.135799543546264
+ATOM_TYPE: Se
+N_TYPE_ATOM: 2
+PSEUDO_POT: ../../../psps/34_Se_16_1.9_2.1_pbe_n_v1.0.psp8
+COORD:    
+32.959304779880675    3.357462537482735    0.000032907355847
+39.278965212349100    3.880850013508314    0.000074221497889
\ No newline at end of file
diff --git a/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.out b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.out
new file mode 100644
index 00000000..0bc4c348
--- /dev/null
+++ b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.out
@@ -0,0 +1,132 @@
+***************************************************************************
+*                       SPARC (version Aug 01, 2023)                      *
+*   Copyright (c) 2020 Material Physics & Mechanics Group, Georgia Tech   *
+*           Distributed under GNU General Public License 3 (GPL)          *
+*                   Start time: Tue Sep  5 12:03:33 2023                  *
+***************************************************************************
+                           Input parameters                                
+***************************************************************************
+CELL: 30.3498965822 0.299199300341885 3.1359066724 
+TWIST_ANGLE: 0.047705 
+FD_GRID: 120 42 25
+FD_ORDER: 12
+BC: D C H
+KPOINT_GRID: 1 3 4
+KPOINT_SHIFT: 0 0 0
+SPIN_TYP: 0
+ELEC_TEMP_TYPE: Fermi-Dirac
+SMEARING: 0.001
+EXCHANGE_CORRELATION: GGA_PBE
+NSTATES: 32
+CHEB_DEGREE: 100
+CHEFSI_BOUND_FLAG: 0
+CALC_STRESS: 1
+MAXIT_SCF: 100
+MINIT_SCF: 2
+MAXIT_POISSON: 3000
+TOL_SCF: 1.00E-06
+POISSON_SOLVER: AAR
+TOL_POISSON: 1.00E-08
+TOL_LANCZOS: 1.00E-02
+TOL_PSEUDOCHARGE: 1.00E-09
+MIXING_VARIABLE: potential
+MIXING_PRECOND: none
+MIXING_PARAMETER: 0.3
+MIXING_HISTORY: 7
+PULAY_FREQUENCY: 1
+PULAY_RESTART: 0
+REFERENCE_CUTOFF: 0.5
+RHO_TRIGGER: 4
+NUM_CHEFSI: 1
+FIX_RAND: 0
+VERBOSITY: 1
+PRINT_FORCES: 1
+PRINT_ATOMS: 1
+PRINT_EIGEN: 0
+PRINT_DENSITY: 0
+PRINT_ENERGY_DENSITY: 0
+OUTPUT_FILE: WSe2_cyclix/temp_run/WSe2_cyclix
+***************************************************************************
+                                Cell                                       
+***************************************************************************
+Volume :  1.0336846642E+03 (Bohr^3)
+***************************************************************************
+                           Parallelization                                 
+***************************************************************************
+NP_SPIN_PARAL: 1
+NP_KPOINT_PARAL: 7
+NP_BAND_PARAL: 6
+NP_DOMAIN_PARAL: 1 1 1
+NP_DOMAIN_PHI_PARAL: 8 3 2
+EIG_SERIAL_MAXNS: 1500
+***************************************************************************
+                             Initialization                                
+***************************************************************************
+Number of processors               :  48
+Mesh spacing in x-direction        :  0.252916 (Bohr)
+Mesh spacing in y-direction        :  0.00712379 (Bohr)
+Mesh spacing in z-direction        :  0.125436 (Bohr)
+Output printed to                  :  WSe2_cyclix/temp_run/WSe2_cyclix.out
+Total number of atom types         :  2
+Total number of atoms              :  3
+Total number of electrons          :  46
+Atom type 1  (valence electrons)   :  W 14
+Pseudopotential                    :  ../psps/74_W_14_2.3_2.5_pbe_n_v1.0.psp8
+Atomic mass                        :  183.84
+Pseudocharge radii of atom type 1  :  7.33 0.52 7.15 (x, y, z dir)
+Number of atoms of type 1          :  1
+Atom type 2  (valence electrons)   :  Se 16
+Pseudopotential                    :  ../psps/34_Se_16_1.9_2.1_pbe_n_v1.0.psp8
+Atomic mass                        :  78.971
+Pseudocharge radii of atom type 2  :  7.33 0.53 7.15 (x, y, z dir)
+Number of atoms of type 2          :  2
+Estimated total memory usage       :  3.06 GB
+Estimated memory per processor     :  65.31 MB
+===================================================================
+                    Self Consistent Field (SCF#1)                     
+===================================================================
+Iteration     Free Energy (Ha/atom)   SCF Error        Timing (sec)
+1            -1.0183254549E+02        7.613E-02        27.229
+2            -1.0183429867E+02        5.101E-02        7.144
+3            -1.0183502413E+02        3.712E-02        7.149
+4            -1.0183546552E+02        1.235E-02        7.108
+5            -1.0183544523E+02        1.337E-02        9.091
+6            -1.0183551633E+02        2.416E-03        10.378
+7            -1.0183552075E+02        3.882E-04        7.066
+8            -1.0183552093E+02        8.927E-05        6.673
+9            -1.0183552093E+02        2.839E-05        6.699
+10           -1.0183552094E+02        8.981E-06        7.334
+11           -1.0183552092E+02        7.418E-06        7.366
+12           -1.0183552094E+02        6.178E-07        7.352
+Total number of SCF: 12    
+====================================================================
+                    Energy and force calculation                    
+====================================================================
+Free energy per atom               : -1.0183552094E+02 (Ha/atom)
+Total free energy                  : -3.0550656283E+02 (Ha)
+Band structure energy              : -6.1063578525E+01 (Ha)
+Exchange correlation energy        : -5.8895903012E+01 (Ha)
+Self and correction energy         : -2.9306005079E+02 (Ha)
+-Entropy*kb*T                      : -1.0332736669E-12 (Ha)
+Fermi level                        : -1.7486456952E-01 (Ha)
+RMS force                          :  6.4467289927E-03 (Ha/Bohr)
+Maximum force                      :  1.5322782852E-02 (Ha/Bohr)
+Time for force calculation         :  0.415 (sec)
+Maximum stress                     :  5.9636521405E-01 (Ha/Bohr**2)
+Maximum stress equiv. to periodic  :  5.7811301847E+02 (GPa)
+Time for stress calculation        :  0.443 (sec)
+***************************************************************************
+                               Timing info                                 
+***************************************************************************
+Total walltime                     :  115.440 sec
+___________________________________________________________________________
+
+***************************************************************************
+*             Material Physics & Mechanics Group, Georgia Tech            *
+*                       PI: Phanish Suryanarayana                         *
+*               List of contributors: See the documentation               *
+*         Citation: See README.md or the documentation for details        *
+*  Acknowledgements: U.S. DOE SC (DE-SC0019410), U.S. DOE NNSA (ASC)      *
+*      {Preliminary developments: U.S. NSF (1333500,1663244,1553212)}     *
+***************************************************************************
+                                                                           
diff --git a/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.static b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.static
new file mode 100644
index 00000000..86ad756c
--- /dev/null
+++ b/tests/outputs/WSe2_helix_static.sparc/WSe2_cyclix.static
@@ -0,0 +1,17 @@
+***************************************************************************
+                            Atom positions                                 
+***************************************************************************
+Fractional coordinates of W:
+      1.1972993167       0.0000145889       0.9999658380
+Fractional coordinates of Se:
+      1.0915974792       0.3392894457       0.0000104937
+      1.3005058466       0.3291420709       0.0000236683
+Total free energy (Ha): -3.055065628300148E+02
+Atomic forces (Ha/Bohr):
+ -1.9631231365E-04   2.0197957296E-03  -9.3342235956E-04
+ -1.4864255791E-03  -8.3883915724E-04   5.1817786556E-04
+ -1.5069571202E-02  -2.7428580192E-03   4.1524449400E-04
+Stress (Ha/Bohr): 
+  5.9636521405E-01 
+Stress equiv. to all periodic (GPa): 
+  2.5346951054E+00 
diff --git a/tests/test_000_ase_io.py b/tests/test_000_ase_io.py
index ac25eac8..39d25c2d 100644
--- a/tests/test_000_ase_io.py
+++ b/tests/test_000_ase_io.py
@@ -178,9 +178,31 @@ def fake_read_sparc(filename, *args, **kwargs):
     assert atoms.get_chemical_formula() == "Al"
 
 
-def test_ase_io_filetype():
-    """If hacked ase.io.formats.filetype correctly recognized sparc format"""
+def test_ase_io_filetype(fs):
+    """If hacked ase.io.formats.filetype correctly recognized sparc format
+
+    Due to the implementation of ase.io.formats, single file tests should be
+    done on non-empty files
+    """
     import sparc
     from ase.io.formats import filetype
 
+    fs.create_dir("test.sparc")
+    fs.create_file("test.sparc/test.ion")
+    fs.create_file("test.sparc/test.static")
+    fs.create_file("test.sparc/test.geopt")
+    fs.create_file("test.sparc/test.aimd")
+    with open("test.sparc/test.ion", "w") as fd:
+        fd.write("\n")
+    with open("test.sparc/test.static", "w") as fd:
+        fd.write("\n")
+    with open("test.sparc/test.geopt", "w") as fd:
+        fd.write("\n")
+    with open("test.sparc/test.aimd", "w") as fd:
+        fd.write("\n")
+
     assert filetype("test.sparc") == "sparc"
+    assert filetype("test.sparc/test.ion") == "ion"
+    assert filetype("test.sparc/test.static") == "static"
+    assert filetype("test.sparc/test.geopt") == "geopt"
+    assert filetype("test.sparc/test.aimd") == "aimd"
diff --git a/tests/test_bundle_to_calc.py b/tests/test_bundle_to_calc.py
index 2baf04a9..48b1a9e1 100644
--- a/tests/test_bundle_to_calc.py
+++ b/tests/test_bundle_to_calc.py
@@ -53,7 +53,5 @@ def test_multi_file_geopt_read():
     assert len(images) == 7
 
     last = read_sparc(bundle, index=-1, include_all_files=True)
-    assert np.isclose(
-        last.get_potential_energy(), -9.057488887961474 * Hartree, 1e-4
-    )
+    assert np.isclose(last.get_potential_energy(), -9.057488887961474 * Hartree, 1e-4)
     return
diff --git a/tests/test_calculator.py b/tests/test_calculator.py
index e783b48e..1e15f23f 100644
--- a/tests/test_calculator.py
+++ b/tests/test_calculator.py
@@ -36,6 +36,90 @@ def test_h_parameter():
         assert "FD_GRID: 25 25 25" in filecontent
 
 
+def test_xc_parameter():
+    from sparc.calculator import SPARC
+    from ase.build import bulk
+
+    atoms = bulk("Al", cubic=True)
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: GGA_PBE" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="pbe", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: GGA_PBE" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="lda", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: LDA_PZ" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="lda", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: LDA_PZ" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="rpbe", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: GGA_RPBE" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="pbesol", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: GGA_PBEsol" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="pbe0", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: PBE0" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="hf", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: HF" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="hse", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: HSE" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="hse03", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: HSE" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="vdw-df", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: vdWDF1" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="vdw-df2", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: vdWDF2" in filecontent
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        calc = SPARC(xc="scan", directory=tmpdir)
+        calc.write_input(atoms)
+        filecontent = open(Path(tmpdir) / "SPARC.inpt", "r").read()
+        assert "EXCHANGE_CORRELATION: SCAN" in filecontent
+
+
 def test_conflict_param():
     from sparc.calculator import SPARC
     from ase.build import bulk
@@ -48,9 +132,7 @@ def test_conflict_param():
             calc.write_input(atoms)
 
     with tempfile.TemporaryDirectory() as tmpdir:
-        calc = SPARC(
-            h=0.2, directory=tmpdir, FD_GRID=[25, 25, 25], MESH_SPACING=0.4
-        )
+        calc = SPARC(h=0.2, directory=tmpdir, FD_GRID=[25, 25, 25], MESH_SPACING=0.4)
         # FD_GRID and ECUT are conflict, but only detected during the writing
         with pytest.raises(Exception):
             calc.write_input(atoms)
diff --git a/tests/test_docparser.py b/tests/test_docparser.py
index 3749a985..b943d90a 100644
--- a/tests/test_docparser.py
+++ b/tests/test_docparser.py
@@ -118,8 +118,8 @@ def test_json():
     import json
 
     sp = SPARCDocParser(test_doc_dir)
-    json_string = sp.to_json()
-    loaded = json.loads(json_string)
+    # json_string = sp.to_json()
+    loaded = sp.to_dict()
     assert all(
         [
             p in loaded.keys()
@@ -137,7 +137,7 @@ def test_json():
 def test_class_load():
     from sparc.docparser import SPARCDocParser
 
-    sp = SPARCDocParser.from_directory(test_doc_dir)
+    sp = SPARCDocParser.json_from_directory(test_doc_dir)
 
 
 def test_text2value():
@@ -155,18 +155,14 @@ def test_text2value():
     assert np.isclose(
         text2value(" 1 2 3 ", desired_type="integer array"), np.array([1, 2, 3])
     ).all()
-    assert isinstance(
-        text2value(" 1 2 3 ", desired_type="integer array")[0], int
-    )
+    assert isinstance(text2value(" 1 2 3 ", desired_type="integer array")[0], int)
 
     assert np.isclose(
         text2value(" 1. 2. 3. ", desired_type="double array"),
         np.array([1, 2, 3]),
     ).all()
 
-    assert isinstance(
-        text2value(" 1 2 3 ", desired_type="double array")[0], float
-    )
+    assert isinstance(text2value(" 1 2 3 ", desired_type="double array")[0], float)
     assert np.array(
         text2value(" 1 2 3\n 4 5 6", desired_type="double array")
     ).shape == (2, 3)
@@ -179,6 +175,4 @@ def test_text2value():
 def test_docparser_main():
     import subprocess
 
-    subprocess.run(
-        ["python", "-m", "sparc.docparser", f"{test_doc_dir.as_posix()}"]
-    )
+    subprocess.run(["python", "-m", "sparc.docparser", f"{test_doc_dir.as_posix()}"])
diff --git a/tests/test_inpt_parser.py b/tests/test_inpt_parser.py
index 6665062c..6522d325 100644
--- a/tests/test_inpt_parser.py
+++ b/tests/test_inpt_parser.py
@@ -69,8 +69,7 @@ def test_write_inpt():
             "PRINT_DENSITY",
         ):
             assert (
-                data_dict["inpt"]["params"][key]
-                == new_data_dict["inpt"]["params"][key]
+                data_dict["inpt"]["params"][key] == new_data_dict["inpt"]["params"][key]
             )
         for key in ("LATVEC", "LATVEC_SCALE"):
             assert np.isclose(
@@ -84,7 +83,7 @@ def test_cell_conversion():
     from ase.units import Bohr, Angstrom
 
     # 0. invalid
-    # 1. invalid
+    # 1. valid, equivalent to LATVEC = diag(1.5, 1.5, 1.5)
     data_dict = {
         "inpt": {
             "params": {
@@ -92,14 +91,12 @@ def test_cell_conversion():
             }
         }
     }
-    with pytest.raises(KeyError):
-        _inpt_cell_to_ase_cell(data_dict)
+    cell = _inpt_cell_to_ase_cell(data_dict)
+    assert np.isclose(cell, np.eye(3) * 1.5 * Bohr).all()
 
     # 1. invalid
     data_dict = {
-        "inpt": {
-            "params": {"CELL": [1.5, 1.5, 1.5], "LATVEC_SCALE": [1.5, 1.5, 1.5]}
-        }
+        "inpt": {"params": {"CELL": [1.5, 1.5, 1.5], "LATVEC_SCALE": [1.5, 1.5, 1.5]}}
     }
     with pytest.raises(ValueError):
         _inpt_cell_to_ase_cell(data_dict)
@@ -135,3 +132,14 @@ def test_cell_conversion():
     cell = _inpt_cell_to_ase_cell(data_dict)
 
     assert np.isclose(cell, np.diag([8, 8, 9]) * Bohr).all()
+
+    # 6. valid, equivalent to LATVEC = diag(1.5, 1.5, 1.5)
+    data_dict = {
+        "inpt": {
+            "params": {
+                "LATVEC_SCALE": [1.5, 1.5, 1.5],
+            }
+        }
+    }
+    cell = _inpt_cell_to_ase_cell(data_dict)
+    assert np.isclose(cell, np.eye(3) * 1.5 * Bohr).all()
diff --git a/tests/test_ion_parser.py b/tests/test_ion_parser.py
index 759c0d62..1eb324a9 100644
--- a/tests/test_ion_parser.py
+++ b/tests/test_ion_parser.py
@@ -132,9 +132,7 @@ def test_read_write_ion_w_sort(fs):
         fd.write(ion_content)
 
     data_dict = _read_ion("test.ion")
-    assert all(
-        ["ASE-SORT" not in line for line in data_dict["ion"]["comments"]]
-    )
+    assert all(["ASE-SORT" not in line for line in data_dict["ion"]["comments"]])
     assert tuple(data_dict["ion"]["sorting"]["sort"]) == (4, 0, 1, 2, 3)
     assert tuple(data_dict["ion"]["sorting"]["resort"]) == (1, 2, 3, 4, 0)
 
diff --git a/tests/test_parse_atoms.py b/tests/test_parse_atoms.py
index c683bc6b..737c2edb 100644
--- a/tests/test_parse_atoms.py
+++ b/tests/test_parse_atoms.py
@@ -1,5 +1,9 @@
 import pytest
 import numpy as np
+from pathlib import Path
+
+curdir = Path(__file__).parent
+test_output_dir = curdir / "outputs"
 
 
 def test_atoms2dict():
@@ -26,9 +30,7 @@ def test_atoms2dict():
     mol1.set_initial_charges([0.1] * 8)
     with pytest.warns(UserWarning, match="initial charges"):
         # Charge cannot be written
-        adict = atoms_to_dict(
-            mol1, sort=False, direct=True, comments=["Ethane"]
-        )
+        adict = atoms_to_dict(mol1, sort=False, direct=True, comments=["Ethane"])
 
     # Spin
     mol2 = mol.copy()
@@ -219,3 +221,25 @@ def test_relax_from_constraint():
     with pytest.warns(UserWarning, match="only support freezing entire"):
         # Freezing diagonal is not possible
         assert relax_from_constraint(FixedPlane(0, [1, 1, 1])) == {}
+
+
+def test_atoms_pbc_conversion():
+    from sparc.sparc_parsers.atoms import atoms_to_dict
+    from ase.build import molecule, bulk, mx2
+
+    h2 = molecule("H2")
+    sparc_dict = atoms_to_dict(h2)
+
+    assert sparc_dict["inpt"]["params"]["BC"] == "D D D"
+    h2.pbc = True
+
+    sparc_dict = atoms_to_dict(h2)
+    assert sparc_dict["inpt"]["params"]["BC"] == "P P P"
+
+    al = bulk("Al", cubic=True)
+    sparc_dict = atoms_to_dict(al)
+    assert sparc_dict["inpt"]["params"]["BC"] == "P P P"
+
+    mos2 = mx2("MoS2")
+    sparc_dict = atoms_to_dict(mos2)
+    assert sparc_dict["inpt"]["params"]["BC"] == "P P D"
diff --git a/tests/test_parser_utils.py b/tests/test_parser_utils.py
index 9cdf72f1..395ad29f 100644
--- a/tests/test_parser_utils.py
+++ b/tests/test_parser_utils.py
@@ -118,6 +118,4 @@ def test_make_reverse_mapping():
         np.random.shuffle(sort)
         reverse = make_reverse_mapping(sort)
         rere = make_reverse_mapping(reverse)
-        assert np.isclose(
-            sort, rere
-        ).all(), f"Reverse error! {sort}, {reverse}, {rere}"
+        assert np.isclose(sort, rere).all(), f"Reverse error! {sort}, {reverse}, {rere}"
diff --git a/tests/test_psp.py b/tests/test_psp.py
index 6137563c..8a591aa7 100644
--- a/tests/test_psp.py
+++ b/tests/test_psp.py
@@ -82,10 +82,7 @@ def test_pseudo_infer():
     with pytest.raises(NoMatchingPseudopotential):
         infer_pseudo_path("As", ".")
 
-    assert (
-        infer_pseudo_path("As", psp_dir).name
-        == "33_As_15_1.8_2.1_pbe_n_v1.0.psp8"
-    )
+    assert infer_pseudo_path("As", psp_dir).name == "33_As_15_1.8_2.1_pbe_n_v1.0.psp8"
 
     with pytest.raises(NoMatchingPseudopotential):
         infer_pseudo_path("Hf", psp_dir)
diff --git a/tests/test_read_all_examples.py b/tests/test_read_all_examples.py
new file mode 100644
index 00000000..6aeb07ce
--- /dev/null
+++ b/tests/test_read_all_examples.py
@@ -0,0 +1,200 @@
+"""Test reading of all SPARC calculation examples
+
+The test needs to be combined with the downloadable test outputs from SPARC-X's repo
+and will only be activated when the environment variable $SPARC_TESTS_DIR is set
+
+The ref
+"""
+import pytest
+import numpy as np
+from pathlib import Path
+import os
+import tempfile
+import shutil
+
+skipped_names = [
+    "Si2_domain_paral",
+    "Si2_kpt_paral",
+    "SiH4",
+    "SiH4_quick",
+    "H2O_sheet_quick",
+    "H2O_sheet",
+    "CdS_bandstruct",
+]
+
+selected_quick_tests = [
+    "AlSi_orthogonal_quick_scf/standard",
+    "AlSi_primitive_quick_relax/standard",
+    "Cu_FCC/standard",
+    "BaTiO3_quick/standard",
+    "H2O_wire_quick/standard",
+]
+
+
+def test_read_all_tests():
+    """Search all .inpt files within the tests dir."""
+
+    from sparc.io import read_sparc
+
+    skipped_names = []
+    tests_dir = os.environ.get("SPARC_TESTS_DIR", "")
+    print(f"Current test dir is {tests_dir}")
+    if len(tests_dir) == 0:
+        pytest.skip(allow_module_level=True)
+
+    tests_dir = Path(tests_dir)
+    failed_counts = 0
+    for inpt_file in tests_dir.glob("**/*.inpt"):
+        workdir = inpt_file.parent
+        parent_name = inpt_file.parents[1].name
+        if parent_name in skipped_names:
+            continue
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmpdir = Path(tmpdir)
+            for ext in [".ion", ".inpt"]:
+                shutil.copy(workdir / f"{parent_name}{ext}", tmpdir)
+            for ext in [".refout", ".refstatic", ".refgeopt", ".refaimd"]:
+                origin_file = workdir / f"{parent_name}{ext}"
+                new_ext = ext.replace("ref", "")
+                new_file = tmpdir / f"{parent_name}{new_ext}"
+                if origin_file.is_file():
+                    shutil.copy(origin_file, new_file)
+            try:
+                read_sparc(tmpdir)
+                # print("Passed: ", parent_name, workdir)
+            except Exception as e:
+                print("Failed: ", parent_name, workdir)
+                print("\t: Error is ", e)
+                failed_counts += 1
+    if failed_counts > 0:
+        raise RuntimeError("More than 1 test in output read test failed")
+
+
+def test_write_all_inputs():
+    """Search all .inpt files within the tests dir."""
+
+    from sparc.io import read_sparc, read_ion, write_ion
+    from sparc.sparc_parsers.inpt import _read_inpt
+
+    # Skipped tests are to avoid unwanted keywords
+    tests_dir = os.environ.get("SPARC_TESTS_DIR", "")
+    failed_counts = 0
+    print(f"Current test dir is {tests_dir}")
+    if len(tests_dir) == 0:
+        pytest.skip(allow_module_level=True)
+
+    tests_dir = Path(tests_dir)
+    for inpt_file in tests_dir.glob("**/*.inpt"):
+        workdir = inpt_file.parent
+        parent_name = inpt_file.parents[1].name
+        ion_file = inpt_file.with_suffix(".ion")
+        if parent_name in skipped_names:
+            continue
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmpdir = Path(tmpdir)
+            origin_atoms = read_ion(ion_file)
+            origin_inpt_dict = _read_inpt(inpt_file)
+            for key in ["CELL", "LATVEC_SCALE", "LATVEC", "BC"]:
+                origin_inpt_dict["inpt"]["params"].pop(key, None)
+            # Re-write the ion and inpt files
+            try:
+                write_ion(
+                    tmpdir / "test.ion",
+                    origin_atoms,
+                    **origin_inpt_dict["inpt"]["params"],
+                )
+                new_atoms = read_ion(tmpdir / "test.ion")
+                new_inpt_dict = _read_inpt(tmpdir / "test.inpt")
+                assert np.all(origin_atoms.pbc == new_atoms.pbc)
+                for key in origin_inpt_dict["inpt"]["params"].keys():
+                    origin_val = origin_inpt_dict["inpt"]["params"][key]
+                    new_val = new_inpt_dict["inpt"]["params"][key]
+                    if isinstance(origin_val, (int, bool)):
+                        assert origin_val == new_val
+                    elif isinstance(origin_val, float):
+                        assert np.isclose(origin_val, new_val, 1e-6)
+                    elif isinstance(origin_val, str):
+                        assert origin_val == new_val
+                    # Vector types can be list compared
+                    elif isinstance(origin_val, (list, np.ndarray)):
+                        assert np.all(origin_val == new_val)
+
+            except Exception as e:
+                print("Failed: ", parent_name, workdir)
+                print("\t: Error is ", e)
+                failed_counts += 1
+    if failed_counts > 0:
+        raise RuntimeError("More than 1 test in inpt write test failed")
+
+
+def test_quick_examples():
+    """Perform quick tests on selected ref examples"""
+    from sparc.calculator import SPARC
+    from ase.build import molecule
+    from pathlib import Path
+    from sparc.io import read_sparc, read_ion, write_ion
+    from sparc.sparc_parsers.inpt import _read_inpt
+
+    dummy_calc = SPARC()
+    try:
+        cmd = dummy_calc._make_command()
+    except EnvironmentError:
+        print("Skip test since no sparc command found")
+        pytest.skip()
+
+    # Skipped tests are to avoid unwanted keywords
+    tests_dir = os.environ.get("SPARC_TESTS_DIR", "")
+    failed_counts = 0
+    print(f"Current test dir is {tests_dir}")
+    if len(tests_dir) == 0:
+        pytest.skip(allow_module_level=True)
+
+    tests_dir = Path(tests_dir)
+    for test_name in selected_quick_tests:
+        bundle = tests_dir / test_name
+        parent_name = test_name.split("/")[0]
+
+        atoms = read_sparc(bundle, index=0)
+        inpt_file = bundle / f"{parent_name}.inpt"
+        params = _read_inpt(inpt_file)["inpt"]["params"]
+        for key in ["CELL", "LATVEC_SCALE", "LATVEC", "BC"]:
+            params.pop(key, None)
+
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmpdir = Path(tmpdir)
+            os.makedirs(tmpdir, exist_ok=True)
+            for ext in [".ion", ".inpt"]:
+                shutil.copy(bundle / f"{parent_name}{ext}", tmpdir)
+            for ext in [".refout", ".refstatic", ".refgeopt", ".refaimd"]:
+                origin_file = bundle / f"{parent_name}{ext}"
+                new_ext = ext.replace("ref", "")
+                new_file = tmpdir / f"{parent_name}{new_ext}"
+                if origin_file.is_file():
+                    shutil.copy(origin_file, new_file)
+            old_atoms = read_sparc(tmpdir, index=-1, include_all_files=True)
+
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmpdir = Path(tmpdir)
+            calc = SPARC(directory=tmpdir, label=parent_name, **params)
+            calc.calculate(atoms=atoms)
+            new_atoms = read_sparc(tmpdir, index=-1, include_all_files=True)
+
+        print("old atoms", old_atoms)
+        print("new atoms", new_atoms)
+        assert len(old_atoms) == len(new_atoms)
+        assert np.all(old_atoms.pbc) == np.all(new_atoms.pbc)
+        if "energy" in old_atoms.calc.results:
+            assert np.isclose(
+                old_atoms.get_potential_energy(),
+                new_atoms.get_potential_energy(),
+                rtol=1.0e-6,
+                atol=1.0e-3,
+            )
+        if "forces" in old_atoms.calc.results:
+            assert np.isclose(
+                old_atoms.get_forces(), new_atoms.get_forces(), rtol=1.0e-3, atol=1.0e-2
+            ).all()
+        if "stress" in old_atoms.calc.results:
+            assert np.isclose(
+                old_atoms.get_stress(), new_atoms.get_stress(), rtol=1.0e-3, atol=1.0e-2
+            ).all()
diff --git a/tests/test_read_sparc.py b/tests/test_read_sparc.py
index ee89d49c..a217b4a6 100644
--- a/tests/test_read_sparc.py
+++ b/tests/test_read_sparc.py
@@ -16,3 +16,16 @@ def test_read_sparc_all():
             results = read_sparc(bundle)
         else:
             results = read_sparc(bundle, include_all_files=True)
+
+
+def test_atoms_read_pbc():
+    from sparc.io import read_sparc
+    from sparc.sparc_parsers.atoms import atoms_to_dict
+
+    # Case 1: H2O sheet
+    water_sheet = read_sparc(test_output_dir / "H2O_sheet_yz.sparc")
+    assert all(water_sheet.pbc == [False, True, True])
+
+    # Case 2: H2O wire
+    water_wire = read_sparc(test_output_dir / "H2O_wire_z.sparc")
+    assert all(water_wire.pbc == [False, False, True])
diff --git a/tests/test_sparc_bundle.py b/tests/test_sparc_bundle.py
index 546f91a3..51ee995e 100644
--- a/tests/test_sparc_bundle.py
+++ b/tests/test_sparc_bundle.py
@@ -48,9 +48,7 @@ def test_default_psp(monkeypatch):
     def _fake_psp_check(directory):
         return True
 
-    monkeypatch.setattr(
-        sparc_io_bundle, "is_psp_download_complete", _fake_psp_check
-    )
+    monkeypatch.setattr(sparc_io_bundle, "is_psp_download_complete", _fake_psp_check)
 
     from sparc.io import SparcBundle
     from sparc.common import psp_dir as default_psp_dir
@@ -69,9 +67,7 @@ def test_bundle_label():
     sb.label == "Cu_FCC"
     assert sb._indir(".ion").name == "Cu_FCC.ion"
 
-    sb = SparcBundle(
-        directory=test_output_dir / "Cu_FCC.sparc", label="Something"
-    )
+    sb = SparcBundle(directory=test_output_dir / "Cu_FCC.sparc", label="Something")
     sb.label == "Something"
 
     assert sb._indir(ext=".ion").name == "Something.ion"
@@ -121,9 +117,7 @@ def test_read_ion_inpt():
         * Bohr,
     ).all()
 
-    sb = SparcBundle(
-        directory=test_output_dir / "AlSi_primitive_quick_relax.sparc"
-    )
+    sb = SparcBundle(directory=test_output_dir / "AlSi_primitive_quick_relax.sparc")
     atoms = sb._read_ion_and_inpt()
     assert atoms.get_chemical_formula() == "AlSi"
 
@@ -132,9 +126,7 @@ def test_read_ion_inpt():
     assert atoms.get_chemical_formula() == "Fe2"
     assert tuple(atoms.get_initial_magnetic_moments()) == (1.0, 1.0)
 
-    sb = SparcBundle(
-        directory=test_output_dir / "TiO2_orthogonal_quick_md.sparc"
-    )
+    sb = SparcBundle(directory=test_output_dir / "TiO2_orthogonal_quick_md.sparc")
     atoms = sb._read_ion_and_inpt()
     assert atoms.get_chemical_formula() == "O4Ti2"
 
diff --git a/tests/test_static_parser.py b/tests/test_static_parser.py
index e6be7c1b..86107caa 100644
--- a/tests/test_static_parser.py
+++ b/tests/test_static_parser.py
@@ -39,13 +39,9 @@ def test_static_parser(fs):
     static_dict = data_dict["static"]
     assert "atoms" in static_dict
     assert tuple(static_dict["atoms"]["symbols"]) == ("Fe", "Fe")
-    assert np.isclose(
-        static_dict["free energy"], -2.283157353113279e02 * Hartree
-    )
+    assert np.isclose(static_dict["free energy"], -2.283157353113279e02 * Hartree)
     assert static_dict["forces"].shape == (2, 3)
-    assert np.isclose(
-        static_dict["forces"][0, 0], 8.0738249305e-01 * Hartree / Bohr
-    )
+    assert np.isclose(static_dict["forces"][0, 0], 8.0738249305e-01 * Hartree / Bohr)
     assert static_dict["stress"].shape == (6,)