From 8a43eab19f86c8fedf11b86fd932bcd55d984d2e Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Tue, 2 Jun 2026 20:43:55 +0800
Subject: [PATCH 01/13] Fix Molden GTO normalization and coordinate conversion

Title: Improve ABACUS Molden output for wavefunction analysis

Summary:
This PR fixes several Molden conversion issues in tools/molden/molden.py while keeping the default workflow unchanged as much as possible.

Changes:
- Correct the primitive Gaussian coefficient convention when writing Molden GTO data. The NAO-to-GTO fit uses unnormalized radial primitives, while Molden readers usually normalize primitive Gaussian functions internally.
- Fix Cartesian_angstrom coordinate conversion. Coordinates in Angstrom are now converted to Bohr for the Molden [Atoms] AU section by dividing by 0.529177210903.
- Add optional multi-start NAO-to-GTO fitting. A single -r value keeps the old single-start behavior; comma-separated -r values enable multi-start fitting and keep the fit with the lowest nonlinear error.
- Add optional Molden [Nval] output via --write-nval. The values are read from UPF z_valence. This option is disabled by default.

Notes:
- The changes are limited to the Molden converter.

Validation:
- Ran the existing molden.py unit tests successfully.
- Checked that default output does not contain [Nval].
- Checked that --write-nval writes C/O/H valence charges for the PhenolDimer test case.
- Checked that Cartesian_angstrom coordinates are written at the correct Bohr scale.
---
 tools/molden/molden.py | 94 ++++++++++++++++++++++++++++++++++++++----
 1 file changed, 85 insertions(+), 9 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index 10745c04992..2c4fb2352ab 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -211,6 +211,16 @@ def _build_gto(a, c, l, r):
         import numpy as np
         g = c * np.exp(-a * r**2) * r**l
         return g
+
+    def _gto_radial_norm(a, l):
+        """Normalization factor of r^l exp(-a r^2) with normalized spherical harmonics.
+
+        The NAO fit uses unnormalized radial primitives, while Molden readers
+        usually normalize primitive Gaussian functions internally. Therefore the
+        fitted coefficient should be divided by this factor when writing Molden.
+        """
+        import math
+        return math.sqrt(2.0 * (2.0 * a)**(l + 1.5) / math.gamma(l + 1.5))
     
     def __str__(self) -> str:
         """print the GTOs in the Gaussian format. Different CGTO are printed as different section."""
@@ -247,6 +257,7 @@ def molden_all(self) -> str:
                     out += f"{spectra[l]:>25s}{ngto:>8d}{'1.00':>8s}\n"
                     for ig in range(ngto):
                         a, c = NumericalRadial[l][ic][ig]
+                        c /= GTORadials._gto_radial_norm(a, l)
                         out += f"{a:>62.3f} {c:>12.3f}\n"
             out += "\n"
         return out
@@ -261,11 +272,12 @@ def molden(self, it, iat) -> str:
                     out += f"{spectra[l]:>25s}{ngto:>8d}{'1.00':>8s}\n"
                     for ig in range(ngto):
                         a, c = NumericalRadial[l][ic][ig]
+                        c /= GTORadials._gto_radial_norm(a, l)
                         out += f"{a:>62.3f} {c:>12.3f}\n"
         out += "\n"
         return out
 
-def fit_radial_with_gto(nao, ngto, l, r, rel_r=2):
+def _fit_radial_with_gto_result(nao, ngto, l, r, rel_r=2):
     """fit one radial function mapped on grid with GTOs
     
     Args:
@@ -334,16 +346,35 @@ def gto_guess(nao, ngto, l, r, rel_r=2):
     norm_nao = simpson(nao**2 * r**2, x=r)
     norm_gto = simpson(out[0][l][0]**2 * r**2, x=r)
     factor = np.sqrt(norm_nao / norm_gto)
-    print(f"NAO2GTO: Renormalize the CGTO from NAO2GTO method with factor {factor:.4f}")
     c *= factor # renormalize the coefficients to make the norm of GTO equals to that of NAO
 
+    return a, c, err, factor
+
+def _print_fit_radial_with_gto_result(a, c, err, factor, ngto, l, best_rel_r=None):
+    print(f"NAO2GTO: Renormalize the CGTO from NAO2GTO method with factor {factor:.4f}")
+    best_rel_r_line = "" if best_rel_r is None else f"         Best rel_r: {best_rel_r}\n"
+
     print(f"""NAO2GTO: Angular momentum {l}, with {ngto} superposition to fit numerical atomic orbitals on given grid, 
-         Nonlinear fitting error: {err:.4e}
+{best_rel_r_line}         Nonlinear fitting error: {err:.4e}
          Exponential and contraction coefficients of primitive GTOs in a.u.:
 {"a":>10} {"c":>10}\n---------------------""")
     for i in range(ngto):
         print(f"{a[i]:10.6f} {c[i]:10.6f}")
     print(f"\nNAO2GTO: The fitted GTOs are saved in the CGTO instance.")
+
+def fit_radial_with_gto(nao, ngto, l, r, rel_r=2):
+    a, c, err, factor = _fit_radial_with_gto_result(nao, ngto, l, r, rel_r)
+    _print_fit_radial_with_gto_result(a, c, err, factor, ngto, l)
+    return a, c
+
+def fit_radial_with_gto_multistart(nao, ngto, l, r, rel_rs):
+    """Run independent rel_r starts and keep the fit with the lowest error."""
+    results = []
+    for rel_r in rel_rs:
+        a, c, err, factor = _fit_radial_with_gto_result(nao, ngto, l, r, rel_r)
+        results.append((err, rel_r, a, c, factor))
+    err, best_rel_r, a, c, factor = min(results, key=lambda item: item[0])
+    _print_fit_radial_with_gto_result(a, c, err, factor, ngto, l, best_rel_r)
     return a, c
 
 def read_nao(fpath):
@@ -399,6 +430,15 @@ def read_nao(fpath):
 
     return {'elem': elem, 'ecut': ecut, 'rcut': rcut, 'nr': nr, 'dr': dr, 'chi': chi}
 
+def parse_rel_r(rel_r):
+    """Return a single rel_r value, or a list for optional multi-start fitting."""
+    if isinstance(rel_r, str):
+        rel_r = rel_r.strip()
+        if "," in rel_r:
+            return [float(item) for item in rel_r.split(",") if item.strip()]
+        return float(rel_r)
+    return rel_r
+
 def convert_nao_to_gto(fnao, fgto = None, ngto: int = 7, rel_r: float = 2):
     """convert the numerical atomic orbitals to GTOs. Each chi (or say the zeta function)
     corresponds to a CGTO (contracted GTO), and the GTOs are fitted to the radial functions.
@@ -408,6 +448,7 @@ def convert_nao_to_gto(fnao, fgto = None, ngto: int = 7, rel_r: float = 2):
     import numpy as np
     import os
 
+    rel_r = parse_rel_r(rel_r)
     gto = GTORadials()
     # read the numerical atomic orbitals
     nao = read_nao(fnao)
@@ -418,7 +459,10 @@ def convert_nao_to_gto(fnao, fgto = None, ngto: int = 7, rel_r: float = 2):
     for l in range(lmax+1):
         nchi = len(nao["chi"][l])
         for i in range(nchi):
-            a, c = fit_radial_with_gto(nao["chi"][l][i], ngto, l, rgrid, rel_r)
+            if isinstance(rel_r, list):
+                a, c = fit_radial_with_gto_multistart(nao["chi"][l][i], ngto, l, rgrid, rel_r)
+            else:
+                a, c = fit_radial_with_gto(nao["chi"][l][i], ngto, l, rgrid, rel_r)
             gto.register_cgto(a, c, l, symbol, 'a')
     
     # draw the fitted GTOs
@@ -791,6 +835,34 @@ def write_molden_atoms(labels, kinds, labels_kinds_map, coords):
         out += f"{elem:<2s}{i+1:>8d}{ptable(elem):>8d}{coords[i][0]:>15.6f}{coords[i][1]:>15.6f}{coords[i][2]:>15.6f}\n"
     return out
 
+def read_upf_z_valence(fupf):
+    """Read z_valence from a UPF pseudopotential file."""
+    import re
+    with open(fupf, "r") as file:
+        data = file.read()
+    m = re.search(r'\bz_valence\s*=\s*[\'"]([^\'"]+)[\'"]', data)
+    if not m:
+        raise ValueError(f"Cannot find z_valence in UPF file {fupf}")
+    zval = float(m.group(1))
+    if abs(zval - round(zval)) < 1e-8:
+        return int(round(zval))
+    return zval
+
+def write_molden_nval(kinds, nvals):
+    """Write Molden [Nval] effective valence charges per element."""
+    out = "[Nval]\n"
+    for elem, nval in zip(kinds, nvals):
+        out += f"{elem} {nval:g}\n"
+    return out
+
+def write_molden_nval_from_stru(stru, pseudo_dir):
+    """Build [Nval] from STRU species and their UPF z_valence values."""
+    import os
+    kinds = [spec['symbol'] for spec in stru['species']]
+    fpps = [os.path.abspath(os.path.join(pseudo_dir, spec['pp_file'])) for spec in stru["species"]]
+    nvals = [read_upf_z_valence(fpp) for fpp in fpps]
+    return write_molden_nval(kinds, nvals)
+
 def read_abacus_input(finput):
     """Read the ABACUS input file and return the key-value pairs
     
@@ -889,7 +961,7 @@ def indexing_mo(total_gto: GTORadials, labels: list):
                 i += 2*l+1
     return out
 
-def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
+def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_nval=False):
     """Entrance function: generate molden file by reading the outdir of ABACUS, for only LCAO 
     calculation.
     
@@ -921,6 +993,8 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
     _ = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
     stru = read_abacus_stru(kv.get("stru_file", "STRU"))
     out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
+    if write_nval:
+        out += write_molden_nval_from_stru(stru, kv.get("pseudo_dir", "./"))
     
     ####################
     # write the atoms  #
@@ -941,7 +1015,7 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
         coords = np.dot(coords, vec)
     elif stru['coord_type'].startswith("Cartesian_angstrom"):
         # including *_center_xy, *_center_z, *_center_xyz, ... cases
-        coords *= 0.529177249
+        coords /= 0.529177210903
     else:
         raise NotImplementedError(f"Unknown coordinate type {stru['coord_type']}")
     out += write_molden_atoms(labels, kinds, labels_kinds_map, coords.tolist())
@@ -1318,6 +1392,7 @@ def _argparse():
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
     -g, --ngto: the number of GTOs to fit ABACUS NAOs. The default is 7.
     -r, --rel_r: the relative cutoff radius for the GTOs. The default is 2.
+    --write-nval: write the Molden [Nval] section from UPF z_valence.
     -o, --output: the output Molden file name. The default is ABACUS.molden.
     """
     import argparse
@@ -1330,7 +1405,8 @@ def _argparse():
     parser.add_argument("-f", "--folder", type=str, help="the folder of the ABACUS calculation")
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
-    parser.add_argument("-r", "--rel_r", type=int, default=2, help="the relative cutoff radius for the GTOs")
+    parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("--write-nval", action="store_true", help="write the Molden [Nval] section from UPF z_valence")
     parser.add_argument("-o", "--output", type=str, default="ABACUS.molden", help="the output Molden file name")
     args = parser.parse_args()
     return args
@@ -1338,7 +1414,7 @@ def _argparse():
 if __name__ == "__main__":
     #unittest.main(exit=False)
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.write_nval)
     print(" ".join("*"*10).center(80, " "))
     print(f"""MOLDEN: Generated Molden file {args.output} from ABACUS calculation in folder {args.folder}.
 WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore
@@ -1352,4 +1428,4 @@ def _argparse():
 Qin X, Shang H, Xiang H, et al. HONPAS: A linear scaling open-source solution for large system simulations[J].
 International Journal of Quantum Chemistry, 2015, 115(10): 647-655.
 """
-    print(citation, flush=True)
\ No newline at end of file
+    print(citation, flush=True)

From 0bb08be02c0a367af774934aefadbb0c40b8100a Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Tue, 2 Jun 2026 23:52:37 +0800
Subject: [PATCH 02/13] Show default values in molden.py CLI help

---
 tools/molden/molden.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index 2c4fb2352ab..9908c6e7ce7 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -1390,13 +1390,16 @@ def _argparse():
     
     -f, --folder: the folder of the ABACUS calculation, in which the STRU, INPUT, KPT, and OUT* folders are located.
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
-    -g, --ngto: the number of GTOs to fit ABACUS NAOs. The default is 7.
-    -r, --rel_r: the relative cutoff radius for the GTOs. The default is 2.
+    -g, --ngto: the number of GTOs to fit ABACUS NAOs.
+    -r, --rel_r: the relative cutoff radius for the GTOs.
     --write-nval: write the Molden [Nval] section from UPF z_valence.
-    -o, --output: the output Molden file name. The default is ABACUS.molden.
+    -o, --output: the output Molden file name.
     """
     import argparse
-    parser = argparse.ArgumentParser(description="Generate Molden file from ABACUS LCAO calculation via NAO2GTO method")
+    parser = argparse.ArgumentParser(
+        description="Generate Molden file from ABACUS LCAO calculation via NAO2GTO method",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
     welcome = """WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore
 the total number of electrons may not be conserved. Always use after a re-normalization operation.
 Once meet any problem, please submit an issue at: https://github.com/deepmodeling/abacus-develop/issues

From 15941879df734a12af72aafa1773c07bfae4a35c Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Wed, 3 Jun 2026 14:05:20 +0800
Subject: [PATCH 03/13] Refactor Molden writing functions and update CLI

Refactor Molden file writing functions to handle effective nuclear charges instead of valence charges. Update command line arguments and main function to remove Nval section.

 - Write [Pseudo] metadata by default in tools/molden/molden.py
 - Remove the [Nval] command-line option from the core Molden generator
---
 tools/molden/molden.py | 31 +++++++++++++++++--------------
 1 file changed, 17 insertions(+), 14 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index 9908c6e7ce7..81056e1118e 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -1,3 +1,9 @@
+'''
+This script is implemented referring to Molden file standard:
+https://www.theochem.ru.nl/molden/molden_format.html
+https://zorkzou.github.io/Molden2AIM/readme.html
+'''
+
 ####################
 # Impl. of NAO2GTO #
 ####################
@@ -848,20 +854,20 @@ def read_upf_z_valence(fupf):
         return int(round(zval))
     return zval
 
-def write_molden_nval(kinds, nvals):
-    """Write Molden [Nval] effective valence charges per element."""
-    out = "[Nval]\n"
-    for elem, nval in zip(kinds, nvals):
-        out += f"{elem} {nval:g}\n"
+def write_molden_pseudo(kinds, labels_kinds_map, nvals):
+    """Write Molden [Pseudo] effective nuclear charges per atom."""
+    out = "[Pseudo]\n"
+    for i, itype in enumerate(labels_kinds_map):
+        out += f"{kinds[itype]} {i + 1} {nvals[itype]:g}\n"
     return out
 
-def write_molden_nval_from_stru(stru, pseudo_dir):
-    """Build [Nval] from STRU species and their UPF z_valence values."""
+def write_molden_pseudo_from_stru(stru, pseudo_dir, labels_kinds_map):
+    """Build [Pseudo] from STRU species and their UPF z_valence values."""
     import os
     kinds = [spec['symbol'] for spec in stru['species']]
     fpps = [os.path.abspath(os.path.join(pseudo_dir, spec['pp_file'])) for spec in stru["species"]]
     nvals = [read_upf_z_valence(fpp) for fpp in fpps]
-    return write_molden_nval(kinds, nvals)
+    return write_molden_pseudo(kinds, labels_kinds_map, nvals)
 
 def read_abacus_input(finput):
     """Read the ABACUS input file and return the key-value pairs
@@ -961,7 +967,7 @@ def indexing_mo(total_gto: GTORadials, labels: list):
                 i += 2*l+1
     return out
 
-def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_nval=False):
+def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
     """Entrance function: generate molden file by reading the outdir of ABACUS, for only LCAO 
     calculation.
     
@@ -993,8 +999,6 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     _ = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
     stru = read_abacus_stru(kv.get("stru_file", "STRU"))
     out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
-    if write_nval:
-        out += write_molden_nval_from_stru(stru, kv.get("pseudo_dir", "./"))
     
     ####################
     # write the atoms  #
@@ -1019,6 +1023,7 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     else:
         raise NotImplementedError(f"Unknown coordinate type {stru['coord_type']}")
     out += write_molden_atoms(labels, kinds, labels_kinds_map, coords.tolist())
+    out += write_molden_pseudo_from_stru(stru, kv.get("pseudo_dir", "./"), labels_kinds_map)
 
     ####################
     # write the basis  #
@@ -1392,7 +1397,6 @@ def _argparse():
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
     -g, --ngto: the number of GTOs to fit ABACUS NAOs.
     -r, --rel_r: the relative cutoff radius for the GTOs.
-    --write-nval: write the Molden [Nval] section from UPF z_valence.
     -o, --output: the output Molden file name.
     """
     import argparse
@@ -1409,7 +1413,6 @@ def _argparse():
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
-    parser.add_argument("--write-nval", action="store_true", help="write the Molden [Nval] section from UPF z_valence")
     parser.add_argument("-o", "--output", type=str, default="ABACUS.molden", help="the output Molden file name")
     args = parser.parse_args()
     return args
@@ -1417,7 +1420,7 @@ def _argparse():
 if __name__ == "__main__":
     #unittest.main(exit=False)
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.write_nval)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output)
     print(" ".join("*"*10).center(80, " "))
     print(f"""MOLDEN: Generated Molden file {args.output} from ABACUS calculation in folder {args.folder}.
 WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore

From c5e0dc8b78216df1f14489984a98c6df1a47663b Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Wed, 3 Jun 2026 14:09:02 +0800
Subject: [PATCH 04/13] Write [Nval] metadata by default in molden_nval.py

This change modifies the molden_nval.py script to write [Nval] metadata by default instead of [Pseudo] when generating Molden files for Multiwfn.
---
 .../molden/Multiwfn_interface/molden_nval.py  | 59 +++++++++++++++++++
 1 file changed, 59 insertions(+)
 create mode 100644 interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py

diff --git a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py
new file mode 100644
index 00000000000..ac94157c45c
--- /dev/null
+++ b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py	
@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+"""Generate ABACUS Molden files with [Nval] metadata for Multiwfn."""
+
+from pathlib import Path
+import argparse
+import importlib.util
+import os
+
+
+def load_molden_module():
+    """Load tools/molden/molden.py without requiring it to be installed."""
+    molden_path = Path(__file__).resolve().parents[2] / "tools" / "molden" / "molden.py"
+    spec = importlib.util.spec_from_file_location("abacus_molden", molden_path)
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+def write_molden_nval_from_stru(molden, stru, pseudo_dir, labels_kinds_map=None):
+    """Build a Molden [Nval] section from STRU species and UPF z_valence."""
+    kinds = [spec["symbol"] for spec in stru["species"]]
+    fpps = [os.path.abspath(os.path.join(pseudo_dir, spec["pp_file"])) for spec in stru["species"]]
+    nvals = [molden.read_upf_z_valence(fpp) for fpp in fpps]
+    out = "[Nval]\n"
+    for elem, nval in zip(kinds, nvals):
+        out += f"{elem} {nval:g}\n"
+    return out
+
+
+def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden"):
+    """Call tools/molden/molden.py and write [Nval] instead of [Pseudo]."""
+    molden = load_molden_module()
+    old_writer = molden.write_molden_pseudo_from_stru
+    molden.write_molden_pseudo_from_stru = (
+        lambda stru, pseudo_dir, labels_kinds_map:
+        write_molden_nval_from_stru(molden, stru, pseudo_dir, labels_kinds_map)
+    )
+    try:
+        return molden.moldengen(folder, ndigits, ngto, rel_r, fmolden)
+    finally:
+        molden.write_molden_pseudo_from_stru = old_writer
+
+
+def _argparse():
+    parser = argparse.ArgumentParser(
+        description="Generate Molden file from ABACUS LCAO calculation for Multiwfn",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
+    parser.add_argument("-f", "--folder", type=str, help="the folder of the ABACUS calculation")
+    parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
+    parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
+    parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("-o", "--output", type=str, default="ABACUS_Multiwfn.molden", help="the output Molden file name")
+    return parser.parse_args()
+
+
+if __name__ == "__main__":
+    args = _argparse()
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output)

From 46f4b28ebfabf86350125965f902cc3161337e2d Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Wed, 3 Jun 2026 16:28:05 +0800
Subject: [PATCH 05/13] Add a Multiwfn interface wrapper with [Nval] and [Cell]
 enabled by default

---
 .../{molden_nval.py => molden.py}             | 56 ++++++++++---------
 1 file changed, 30 insertions(+), 26 deletions(-)
 rename interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/{molden_nval.py => molden.py} (51%)

diff --git a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py
similarity index 51%
rename from interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py
rename to interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py
index ac94157c45c..ef7709e4457 100644
--- a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden_nval.py	
+++ b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py	
@@ -1,11 +1,18 @@
 #!/usr/bin/env python3
-"""Generate ABACUS Molden files with [Nval] metadata for Multiwfn."""
+"""Generate ABACUS Molden files with [Nval] metadata for Multiwfn.
+
+Keep this script in the Multiwfn_interface directory so it can locate
+tools/molden/molden.py by relative path. Users may add this directory to PATH
+and run chmod +x on this file to use it as a command.
+"""
 
 from pathlib import Path
 import argparse
 import importlib.util
 import os
 
+DEFAULT_FOLDER = "$(pwd)"
+
 
 def load_molden_module():
     """Load tools/molden/molden.py without requiring it to be installed."""
@@ -16,44 +23,41 @@ def load_molden_module():
     return module
 
 
-def write_molden_nval_from_stru(molden, stru, pseudo_dir, labels_kinds_map=None):
-    """Build a Molden [Nval] section from STRU species and UPF z_valence."""
-    kinds = [spec["symbol"] for spec in stru["species"]]
-    fpps = [os.path.abspath(os.path.join(pseudo_dir, spec["pp_file"])) for spec in stru["species"]]
-    nvals = [molden.read_upf_z_valence(fpp) for fpp in fpps]
-    out = "[Nval]\n"
-    for elem, nval in zip(kinds, nvals):
-        out += f"{elem} {nval:g}\n"
-    return out
-
-
-def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden"):
-    """Call tools/molden/molden.py and write [Nval] instead of [Pseudo]."""
+def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden", write_cell=True, with_nval=True):
+    """Call tools/molden/molden.py with Multiwfn-oriented [Nval] metadata."""
     molden = load_molden_module()
-    old_writer = molden.write_molden_pseudo_from_stru
-    molden.write_molden_pseudo_from_stru = (
-        lambda stru, pseudo_dir, labels_kinds_map:
-        write_molden_nval_from_stru(molden, stru, pseudo_dir, labels_kinds_map)
-    )
-    try:
-        return molden.moldengen(folder, ndigits, ngto, rel_r, fmolden)
-    finally:
-        molden.write_molden_pseudo_from_stru = old_writer
+    pseudo = "nval" if with_nval else "none"
+    return molden.moldengen(folder, ndigits, ngto, rel_r, fmolden, write_cell, pseudo)
 
 
 def _argparse():
+    def str_to_bool(value):
+        if isinstance(value, bool):
+            return value
+        value = value.lower()
+        if value in ("true", "t", "yes", "y", "1"):
+            return True
+        if value in ("false", "f", "no", "n", "0"):
+            return False
+        raise argparse.ArgumentTypeError("expected true or false")
+
     parser = argparse.ArgumentParser(
         description="Generate Molden file from ABACUS LCAO calculation for Multiwfn",
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
     )
-    parser.add_argument("-f", "--folder", type=str, help="the folder of the ABACUS calculation")
+    parser.add_argument("-f", "--folder", type=str, default=DEFAULT_FOLDER, help="the folder of the ABACUS calculation")
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("--with-cell", type=str_to_bool, default=True, help="whether to write the Molden [Cell] section")
+    parser.add_argument("--with-Nval", dest="with_nval", type=str_to_bool, default=True, help="whether to write the Molden [Nval] section")
     parser.add_argument("-o", "--output", type=str, default="ABACUS_Multiwfn.molden", help="the output Molden file name")
-    return parser.parse_args()
+    args = parser.parse_args()
+    if args.folder == DEFAULT_FOLDER:
+        args.folder = os.getcwd()
+    return args
 
 
 if __name__ == "__main__":
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval)

From cf2e3935b4b7772950fc5dbcde614fd9e4b14c6c Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Wed, 3 Jun 2026 16:43:49 +0800
Subject: [PATCH 06/13] Add option to include [Pseudo] metadata in Molden

---
 .../molden/Multiwfn_interface/molden.py       | 46 +++++++++++++++++--
 1 file changed, 42 insertions(+), 4 deletions(-)

diff --git a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py
index ef7709e4457..bccca0ee3e7 100644
--- a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py	
+++ b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py	
@@ -23,11 +23,48 @@ def load_molden_module():
     return module
 
 
-def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden", write_cell=True, with_nval=True):
+def write_molden_nval(kinds, zvals):
+    """Write Multiwfn-oriented Molden [Nval] values per element."""
+    out = "[Nval]\n"
+    for elem, zval in zip(kinds, zvals):
+        out += f"{elem} {zval:g}\n"
+    return out
+
+
+def make_molden_nval_section(molden, folder):
+    """Build [Nval] from UPF z_valence without adding this logic to tools/molden."""
+    folder = os.path.abspath(folder)
+    kv = molden.read_abacus_input(os.path.join(folder, "INPUT"))
+    stru = molden.read_abacus_stru(os.path.join(folder, kv.get("stru_file", "STRU")))
+    pseudo_dir = kv.get("pseudo_dir", "./")
+    if not os.path.isabs(pseudo_dir):
+        pseudo_dir = os.path.join(folder, pseudo_dir)
+    kinds = [spec["symbol"] for spec in stru["species"]]
+    fpps = [os.path.abspath(os.path.join(pseudo_dir, spec["pp_file"])) for spec in stru["species"]]
+    zvals = [molden.read_upf_z_valence(fpp) for fpp in fpps]
+    return write_molden_nval(kinds, zvals)
+
+
+def insert_molden_nval(text, nval):
+    """Insert [Nval] after [Atoms] and before [GTO]."""
+    marker = "[GTO]\n"
+    index = text.find(marker)
+    if index < 0:
+        raise ValueError("Cannot insert [Nval]: [GTO] section was not found")
+    return text[:index] + nval + text[index:]
+
+
+def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden", write_cell=True, with_nval=True, with_pseudo=False):
     """Call tools/molden/molden.py with Multiwfn-oriented [Nval] metadata."""
     molden = load_molden_module()
-    pseudo = "nval" if with_nval else "none"
-    return molden.moldengen(folder, ndigits, ngto, rel_r, fmolden, write_cell, pseudo)
+    pseudo = "pseudo" if with_pseudo else "none"
+    out = molden.moldengen(folder, ndigits, ngto, rel_r, fmolden, write_cell, pseudo)
+    if with_nval:
+        out = insert_molden_nval(out, make_molden_nval_section(molden, folder))
+        with open(fmolden, "w") as file:
+            file.write(out)
+    return out
+
 
 
 def _argparse():
@@ -51,6 +88,7 @@ def str_to_bool(value):
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
     parser.add_argument("--with-cell", type=str_to_bool, default=True, help="whether to write the Molden [Cell] section")
     parser.add_argument("--with-Nval", dest="with_nval", type=str_to_bool, default=True, help="whether to write the Molden [Nval] section")
+    parser.add_argument("--with-pseudo", type=str_to_bool, default=False, help="whether to write the Molden [Pseudo] section")
     parser.add_argument("-o", "--output", type=str, default="ABACUS_Multiwfn.molden", help="the output Molden file name")
     args = parser.parse_args()
     if args.folder == DEFAULT_FOLDER:
@@ -60,4 +98,4 @@ def str_to_bool(value):
 
 if __name__ == "__main__":
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval, args.with_pseudo)

From ac15e75df94e4440e4f0f210975da499e160ec91 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Wed, 3 Jun 2026 16:47:39 +0800
Subject: [PATCH 07/13] Refactor Molden functions for pseudopotential handling

---
 tools/molden/molden.py | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index 81056e1118e..50ef058cf26 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -854,20 +854,24 @@ def read_upf_z_valence(fupf):
         return int(round(zval))
     return zval
 
-def write_molden_pseudo(kinds, labels_kinds_map, nvals):
+def write_molden_pseudo(kinds, labels_kinds_map, zvals):
     """Write Molden [Pseudo] effective nuclear charges per atom."""
     out = "[Pseudo]\n"
     for i, itype in enumerate(labels_kinds_map):
-        out += f"{kinds[itype]} {i + 1} {nvals[itype]:g}\n"
+        out += f"{kinds[itype]} {i + 1} {zvals[itype]:g}\n"
     return out
 
-def write_molden_pseudo_from_stru(stru, pseudo_dir, labels_kinds_map):
-    """Build [Pseudo] from STRU species and their UPF z_valence values."""
+def write_molden_pseudopotential_from_stru(stru, pseudo_dir, labels_kinds_map, pseudo="pseudo"):
+    """Build optional pseudopotential-related sections from UPF z_valence values."""
     import os
+    if pseudo == "none":
+        return ""
     kinds = [spec['symbol'] for spec in stru['species']]
     fpps = [os.path.abspath(os.path.join(pseudo_dir, spec['pp_file'])) for spec in stru["species"]]
-    nvals = [read_upf_z_valence(fpp) for fpp in fpps]
-    return write_molden_pseudo(kinds, labels_kinds_map, nvals)
+    zvals = [read_upf_z_valence(fpp) for fpp in fpps]
+    if pseudo == "pseudo":
+        return write_molden_pseudo(kinds, labels_kinds_map, zvals)
+    raise ValueError(f"Unknown Molden pseudopotential mode {pseudo}")
 
 def read_abacus_input(finput):
     """Read the ABACUS input file and return the key-value pairs
@@ -967,7 +971,7 @@ def indexing_mo(total_gto: GTORadials, labels: list):
                 i += 2*l+1
     return out
 
-def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
+def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_cell=False, pseudo="none"):
     """Entrance function: generate molden file by reading the outdir of ABACUS, for only LCAO 
     calculation.
     
@@ -998,7 +1002,8 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
     kv = read_abacus_input("INPUT")
     _ = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
     stru = read_abacus_stru(kv.get("stru_file", "STRU"))
-    out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
+    if write_cell:
+        out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
     
     ####################
     # write the atoms  #
@@ -1023,7 +1028,7 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden"):
     else:
         raise NotImplementedError(f"Unknown coordinate type {stru['coord_type']}")
     out += write_molden_atoms(labels, kinds, labels_kinds_map, coords.tolist())
-    out += write_molden_pseudo_from_stru(stru, kv.get("pseudo_dir", "./"), labels_kinds_map)
+    out += write_molden_pseudopotential_from_stru(stru, kv.get("pseudo_dir", "./"), labels_kinds_map, pseudo)
 
     ####################
     # write the basis  #
@@ -1397,6 +1402,8 @@ def _argparse():
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
     -g, --ngto: the number of GTOs to fit ABACUS NAOs.
     -r, --rel_r: the relative cutoff radius for the GTOs.
+    --with-cell: write the Molden [Cell] section.
+    --with-pseudo: write the Molden [Pseudo] section from UPF z_valence.
     -o, --output: the output Molden file name.
     """
     import argparse
@@ -1413,6 +1420,8 @@ def _argparse():
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("--with-cell", action="store_true", help="write the Molden [Cell] section")
+    parser.add_argument("--with-pseudo", action="store_true", help="write the Molden [Pseudo] section from UPF z_valence")
     parser.add_argument("-o", "--output", type=str, default="ABACUS.molden", help="the output Molden file name")
     args = parser.parse_args()
     return args
@@ -1420,7 +1429,8 @@ def _argparse():
 if __name__ == "__main__":
     #unittest.main(exit=False)
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output)
+    pseudo = "pseudo" if args.with_pseudo else "none"
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, pseudo)
     print(" ".join("*"*10).center(80, " "))
     print(f"""MOLDEN: Generated Molden file {args.output} from ABACUS calculation in folder {args.folder}.
 WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore

From 0619f19e6ea02555d3e870b59a97916aed99a427 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 02:08:33 +0800
Subject: [PATCH 08/13] Enhance ABACUS lowf file handling and moldengen
 function

Updated regex pattern for reading ABACUS lowf files to support new naming conventions. Added functions to handle multiple k points and modified the moldengen function to accept ikpoint parameter.
---
 tools/molden/molden.py | 245 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 207 insertions(+), 38 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index 50ef058cf26..e9662b0f709 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -506,14 +506,14 @@ def write_molden_gto(total_gto: GTORadials, labels: list):
         out += total_gto.molden(l, iat + 1) # iat starts from 0, molden starts from 1
     return out
 
-def read_abacus_lowf(flowf, pat=r'^WFC_NAO_(K\d+|GAMMA(1|2)).txt$'):
-    """Read the ABACUS WFC_NAO_(K|GAMMA) file, which contains the coefficients of MOs.
+def read_abacus_lowf(flowf, pat=r'^(WFC_NAO_(K\d+|GAMMA\d+)(_ION\d+)?|wf((k\d+)?(s\d+)?|(s\d+)?(k\d+)?)(g\d+)?_nao)\.txt$'):
+    """Read the ABACUS LCAO wavefunction file, which contains the coefficients of MOs.
     Please note, only Gamma-only calculation is supported. For ABACUS version earlier
     than 3.7.1, the `flowf` file is something like LOWF_K_. Manually change the file
     name and run is okay.
     
     Args:
-        flowf (str): the file name of the ABACUS WFC_NAO_(K|GAMMA) file
+        flowf (str): the file name of the ABACUS LCAO wavefunction file
         pat (str): the pattern to match the file name
     
     Returns:
@@ -555,7 +555,8 @@ def read_abacus_lowf(flowf, pat=r'^WFC_NAO_(K\d+|GAMMA(1|2)).txt$'):
         i += 1
 
     # check if the data we collected has the correct number of elements
-    if "WFC_NAO_K" in flowf: # multi-k case, the coef will be complex
+    basename = os.path.basename(flowf)
+    if "WFC_NAO_K" in basename or re.match(r"^wf(k\d+|s\d+k\d+)", basename): # multi-k case, the coef will be complex
         data = [complex(data[i], data[i+1]) for i in range(0, len(data), 2)]
         data = [d.real for d in data]
     data = np.array(data).reshape(nbands, nlocal)
@@ -566,6 +567,132 @@ def read_abacus_lowf(flowf, pat=r'^WFC_NAO_(K\d+|GAMMA(1|2)).txt$'):
 
     return nbands, nlocal, occ, band, ener, data
 
+def append_molden_filename_suffix(fmolden, suffix):
+    """Append a suffix before the Molden file extension."""
+    import os
+
+    root, ext = os.path.splitext(fmolden)
+    if ext == "":
+        ext = ".molden"
+    return f"{root}_{suffix}{ext}"
+
+def find_abacus_lowf_files(nspin, ikpoint=1, nkpoints=1):
+    """Find ABACUS LCAO wavefunction files in OUT.*.
+
+    The converter keeps compatibility with older WFC_NAO_* names and newer
+    wf*_nao.txt names. The nspin=4 spinor/SOC format is not mapped to Molden
+    because the current writer only emits restricted/Alpha/Beta MOs.
+    """
+    import glob
+    import os
+    import re
+
+    if nspin == 4:
+        raise NotImplementedError("nspin=4/SOC spinor wavefunctions are not supported by this Molden writer")
+    if nspin not in (1, 2):
+        raise NotImplementedError(f"Unsupported nspin={nspin}")
+    if ikpoint < 1:
+        raise ValueError("ikpoint must be a positive integer when selecting one k point")
+    if ikpoint > nkpoints:
+        raise ValueError(f"ikpoint={ikpoint} is out of range for {nkpoints} k point(s)")
+
+    def newest(pattern):
+        def key(path):
+            nums = [int(x) for x in re.findall(r"\d+", os.path.basename(path))]
+            return nums, path
+        matches = sorted(glob.glob(pattern), key=key)
+        return matches[-1] if matches else None
+
+    def old_k_file(spin):
+        """Select an old WFC_NAO_K* file for a spin channel and k point.
+
+        Older ABACUS files flatten spin and k into WFC_NAO_KN.txt. For nspin=2,
+        files are ordered as all k points of spin 1 followed by all k points of
+        spin 2, so the requested file index is ik + (spin - 1) * nk.
+        """
+        files_by_k = {}
+        for path in glob.glob("WFC_NAO_K*.txt"):
+            match = re.match(r"^WFC_NAO_K(\d+)(?:_ION(\d+))?\.txt$", os.path.basename(path))
+            if match is None:
+                continue
+            ik = int(match.group(1))
+            ion = int(match.group(2) or 0)
+            if ik not in files_by_k or (ion, path) > files_by_k[ik][0]:
+                files_by_k[ik] = ((ion, path), path)
+        target = ikpoint if nspin == 1 else ikpoint + (spin - 1) * nkpoints
+        return files_by_k.get(target, (None, None))[1]
+
+    if nspin == 1:
+        pattern_groups = [[
+            "OLD_WFC_NAO_K_SPIN1",
+            f"wfk{ikpoint}_nao.txt",
+            f"WFC/wfk{ikpoint}g*_nao.txt",
+        ]]
+        if ikpoint == 1:
+            pattern_groups[0] = [
+                "WFC_NAO_GAMMA1.txt",
+                "WFC_NAO_GAMMA1_ION*.txt",
+                "wf_nao.txt",
+                "WFC/wfg*_nao.txt",
+            ] + pattern_groups[0]
+    else:
+        pattern_groups = [
+            [
+                "OLD_WFC_NAO_K_SPIN1",
+                f"wfs1k{ikpoint}_nao.txt",
+                f"wfk{ikpoint}s1_nao.txt",
+                f"WFC/wfs1k{ikpoint}g*_nao.txt",
+                f"WFC/wfk{ikpoint}s1g*_nao.txt",
+            ],
+            [
+                "OLD_WFC_NAO_K_SPIN2",
+                f"wfs2k{ikpoint}_nao.txt",
+                f"wfk{ikpoint}s2_nao.txt",
+                f"WFC/wfs2k{ikpoint}g*_nao.txt",
+                f"WFC/wfk{ikpoint}s2g*_nao.txt",
+            ],
+        ]
+        if ikpoint == 1:
+            pattern_groups[0] = [
+                "WFC_NAO_GAMMA1.txt",
+                "WFC_NAO_GAMMA1_ION*.txt",
+                "wfs1_nao.txt",
+                "WFC/wfs1g*_nao.txt",
+            ] + pattern_groups[0]
+            pattern_groups[1] = [
+                "WFC_NAO_GAMMA2.txt",
+                "WFC_NAO_GAMMA2_ION*.txt",
+                "wfs2_nao.txt",
+                "WFC/wfs2g*_nao.txt",
+            ] + pattern_groups[1]
+
+    files = []
+    for patterns in pattern_groups:
+        found = None
+        for pattern in patterns:
+            if pattern == "OLD_WFC_NAO_K_SPIN1":
+                found = old_k_file(1)
+            elif pattern == "OLD_WFC_NAO_K_SPIN2":
+                found = old_k_file(2)
+            else:
+                found = newest(pattern)
+            if found is not None:
+                break
+        if found is None:
+            raise FileNotFoundError(f"Cannot find ABACUS LCAO wavefunction file. Tried: {', '.join(patterns)}")
+        files.append(found)
+    return files
+
+def find_abacus_lowf_file_groups(nspin, ikpoint=1, nkpoints=1):
+    """Find wavefunction files grouped by selected k point."""
+    if ikpoint == -1:
+        ikpoints = range(1, nkpoints + 1)
+    elif ikpoint > 0:
+        ikpoints = [ikpoint]
+    else:
+        raise ValueError("ikpoint must be a positive integer, or -1 to convert all k points")
+    return [(ik, find_abacus_lowf_files(nspin, ik, nkpoints)) for ik in ikpoints]
+
 def write_molden_mo(coefs, mo_map, occ, ener, spin=0, ndigits=3):
     """Write Molden [MO] section from the coefficients, occupations and energies.
 
@@ -899,31 +1026,34 @@ def read_abacus_input(finput):
     return kv
 
 def read_abacus_kpt(fkpt):
-    """the way to organize information of KPT file of ABACUS still has some degree of freedom.
-    However, the only one wanted should have content like the following:
+    """Read the number of k points from an ABACUS KPT file.
+
+    The automatic mesh format is expected to look like:
     K_POINTS
     0
     Gamma
     1 1 1 0 0 0
 
-    in which the "Gamma" can be replaced by "MP" but the number of kpoints should be 1.
-    In the future the multiple kpoints is planned to be supported in a relatively naive way that
-    simply combining all MOs at different kpoints together but not for now. The occupation of MO
-    at different kpoints is already multiplied by the weight of kpoints, is it expected?
-
-    This function is not really read kpoints and return something, instead, it is for assert
-    the number of kpoints is 1.
+    In this script, multiple k points are converted one k point at a time. The
+    k-point number is used only to select the corresponding ABACUS
+    wavefunction file; different k points are not merged into one Molden file.
     """
     with open(fkpt, 'r') as file:
         lines = file.readlines()
-    lines = [line.strip() for line in lines]
-    if lines[0] == "K_POINTS":
-        if lines[1] == "0":
-            if lines[2] in ["Gamma", "MP"]:
-                if lines[3] == "1 1 1 0 0 0":
-                    return
-    raise ValueError(f"Invalid KPT file {fkpt}. Presently only 1 kpoint calculation \
-(implicit or explicit) Gamma-only calculation is supported.")
+    lines = [line.split("#", 1)[0].strip() for line in lines]
+    lines = [line for line in lines if line]
+    if len(lines) < 2 or lines[0] != "K_POINTS":
+        raise ValueError(f"Invalid KPT file {fkpt}.")
+    if lines[1] == "0":
+        if len(lines) < 4 or lines[2] not in ["Gamma", "MP"]:
+            raise ValueError(f"Invalid KPT file {fkpt}. Only Gamma/MP automatic mesh is supported.")
+        mesh = [int(x) for x in lines[3].split()[:3]]
+        nkpoints = mesh[0] * mesh[1] * mesh[2]
+    else:
+        nkpoints = int(lines[1])
+    if nkpoints < 1:
+        raise ValueError(f"Invalid KPT file {fkpt}. Number of k points should be positive.")
+    return nkpoints
 
 def CondonShortleyPhase(index):
     """Imposing the Condon-Shortley phase on the MO index. 
@@ -971,7 +1101,7 @@ def indexing_mo(total_gto: GTORadials, labels: list):
                 i += 2*l+1
     return out
 
-def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_cell=False, pseudo="none"):
+def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_cell=False, pseudo="none", ikpoint=1, return_files=False, output_is_default=False):
     """Entrance function: generate molden file by reading the outdir of ABACUS, for only LCAO 
     calculation.
     
@@ -980,15 +1110,23 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
         ndigits (int): the number of digits to be printed for the coefficients
         ngto (int): the number of GTOs to be fitted to the numerical atomic orbitals
         fmolden (str): the file name of the molden file
+        ikpoint (int): selected k point index; -1 converts all k points separately
     
     Returns:
         str: the content of the molden file"""
     import os
     import numpy as np
 
+    folder = os.path.abspath(folder)
     files = os.listdir(folder)
     assert ("STRU" in files) and ("INPUT" in files) and ("KPT" in files),\
         "STRU, INPUT, KPT files are required in the folder."
+    input_file = os.path.join(folder, "INPUT")
+    print(f"MOLDEN: Processing ABACUS calculation directory {folder}")
+    kv = read_abacus_input(input_file)
+    nspin = int(kv.get("nspin", 1))
+    if nspin == 4:
+        raise NotImplementedError("nspin=4/SOC spinor wavefunctions are not supported by this Molden writer")
     cwd = os.path.abspath(os.getcwd())
     os.chdir(folder)
     ####################
@@ -999,8 +1137,7 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     ####################
     # write the cell   #
     ####################
-    kv = read_abacus_input("INPUT")
-    _ = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
+    nkpoints = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
     stru = read_abacus_stru(kv.get("stru_file", "STRU"))
     if write_cell:
         out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
@@ -1049,20 +1186,38 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     ####################
     suffix = kv.get("suffix", "ABACUS")
     out_dir = ".".join(["OUT", suffix])
+    print(f"MOLDEN: Parsing ABACUS output directory {os.path.abspath(out_dir)}")
     os.chdir(out_dir)
-    nspin = int(kv.get("nspin", 1))
-    mo_files = "WFC_NAO_GAMMA1.txt"
-    mo_files = "WFC_NAO_K1.txt" if mo_files not in os.listdir() else mo_files
-    mo_files = [mo_files] if nspin == 1 else [mo_files, mo_files.replace("1.txt", "2.txt")]
-    for ispin, mo_file in enumerate(mo_files):
-        nbands, nlocal, occ, band, ener, data_np = read_abacus_lowf(mo_file)
-        coefs = data_np.tolist()
-        out += write_molden_mo(coefs, mo_map, occ, ener, spin=ispin, ndigits=ndigits)
-
+    mo_file_groups = find_abacus_lowf_file_groups(nspin, ikpoint, nkpoints)
+    default_fmolden = output_is_default or fmolden in (None, "ABACUS.molden")
+    if fmolden is None:
+        fmolden = "ABACUS.molden"
+    outputs = []
+    for selected_ikpoint, mo_files in mo_file_groups:
+        out_ik = out
+        print(f"MOLDEN: Converting k point {selected_ikpoint} of {nkpoints}")
+        for ispin, mo_file in enumerate(mo_files):
+            print(f"MOLDEN: Parsing wavefunction file {os.path.abspath(mo_file)}")
+            nbands, nlocal, occ, band, ener, data_np = read_abacus_lowf(mo_file)
+            coefs = data_np.tolist()
+            out_ik += write_molden_mo(coefs, mo_map, occ, ener, spin=ispin, ndigits=ndigits)
+        fmolden_ik = fmolden
+        if len(mo_file_groups) > 1 or (nkpoints > 1 and default_fmolden):
+            fmolden_ik = append_molden_filename_suffix(fmolden, f"k{selected_ikpoint}")
+        outputs.append((fmolden_ik, out_ik))
     os.chdir(cwd)
-    with open(fmolden, "w") as file:
-        file.write(out)
-    return out
+    written = []
+    for fmolden_ik, out_ik in outputs:
+        fmolden_ik = os.path.abspath(fmolden_ik)
+        print(f"MOLDEN: Writing Molden file {fmolden_ik}")
+        with open(fmolden_ik, "w") as file:
+            file.write(out_ik)
+        written.append((fmolden_ik, out_ik))
+    if return_files:
+        return written
+    if len(written) == 1:
+        return written[0][1]
+    return [out_ik for _, out_ik in written]
 
 import unittest
 class TestNAO2GTO(unittest.TestCase):
@@ -1402,6 +1557,7 @@ def _argparse():
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
     -g, --ngto: the number of GTOs to fit ABACUS NAOs.
     -r, --rel_r: the relative cutoff radius for the GTOs.
+    --ikpoint: the k point to convert. Use -1 to convert all k points into separate Molden files.
     --with-cell: write the Molden [Cell] section.
     --with-pseudo: write the Molden [Pseudo] section from UPF z_valence.
     -o, --output: the output Molden file name.
@@ -1420,6 +1576,7 @@ def _argparse():
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("--ikpoint", type=int, default=1, help="the k point to convert; use -1 to convert all k points into separate Molden files")
     parser.add_argument("--with-cell", action="store_true", help="write the Molden [Cell] section")
     parser.add_argument("--with-pseudo", action="store_true", help="write the Molden [Pseudo] section from UPF z_valence")
     parser.add_argument("-o", "--output", type=str, default="ABACUS.molden", help="the output Molden file name")
@@ -1430,9 +1587,21 @@ def _argparse():
     #unittest.main(exit=False)
     args = _argparse()
     pseudo = "pseudo" if args.with_pseudo else "none"
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, pseudo)
+    generated = moldengen(
+        args.folder,
+        args.ndigits,
+        args.ngto,
+        args.rel_r,
+        args.output,
+        args.with_cell,
+        pseudo,
+        args.ikpoint,
+        return_files=True,
+        output_is_default=args.output == "ABACUS.molden",
+    )
+    generated_files = ", ".join(path for path, _ in generated)
     print(" ".join("*"*10).center(80, " "))
-    print(f"""MOLDEN: Generated Molden file {args.output} from ABACUS calculation in folder {args.folder}.
+    print(f"""MOLDEN: Generated Molden file(s) {generated_files} from ABACUS calculation in folder {args.folder}.
 WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore
 the total number of electrons may not be conserved. Always use after a re-normalization operation.""")
     citation = """If you use this script in your research, please cite the following paper:\n

From 070a9f737fb283c295c4b2dafdda6f7d8f81dc67 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 02:10:25 +0800
Subject: [PATCH 09/13] Molden wavefunction filename compatibility

---
 .../molden/Multiwfn_interface/molden.py       | 33 +++++++++++++++----
 1 file changed, 26 insertions(+), 7 deletions(-)

diff --git a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py
index bccca0ee3e7..708eab87190 100644
--- a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py	
+++ b/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py	
@@ -54,16 +54,34 @@ def insert_molden_nval(text, nval):
     return text[:index] + nval + text[index:]
 
 
-def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS.molden", write_cell=True, with_nval=True, with_pseudo=False):
+def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS_Multiwfn.molden", write_cell=True, with_nval=True, with_pseudo=False, ikpoint=1):
     """Call tools/molden/molden.py with Multiwfn-oriented [Nval] metadata."""
     molden = load_molden_module()
     pseudo = "pseudo" if with_pseudo else "none"
-    out = molden.moldengen(folder, ndigits, ngto, rel_r, fmolden, write_cell, pseudo)
+    outputs = molden.moldengen(
+        folder,
+        ndigits,
+        ngto,
+        rel_r,
+        fmolden,
+        write_cell,
+        pseudo,
+        ikpoint,
+        return_files=True,
+        output_is_default=fmolden == "ABACUS_Multiwfn.molden",
+    )
     if with_nval:
-        out = insert_molden_nval(out, make_molden_nval_section(molden, folder))
-        with open(fmolden, "w") as file:
-            file.write(out)
-    return out
+        nval = make_molden_nval_section(molden, folder)
+        outputs_with_nval = []
+        for output_file, out in outputs:
+            out = insert_molden_nval(out, nval)
+            with open(output_file, "w") as file:
+                file.write(out)
+            outputs_with_nval.append((output_file, out))
+        outputs = outputs_with_nval
+    if len(outputs) == 1:
+        return outputs[0][1]
+    return [out for _, out in outputs]
 
 
 
@@ -86,6 +104,7 @@ def str_to_bool(value):
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
+    parser.add_argument("--ikpoint", type=int, default=1, help="the k point to convert; use -1 to convert all k points into separate Molden files")
     parser.add_argument("--with-cell", type=str_to_bool, default=True, help="whether to write the Molden [Cell] section")
     parser.add_argument("--with-Nval", dest="with_nval", type=str_to_bool, default=True, help="whether to write the Molden [Nval] section")
     parser.add_argument("--with-pseudo", type=str_to_bool, default=False, help="whether to write the Molden [Pseudo] section")
@@ -98,4 +117,4 @@ def str_to_bool(value):
 
 if __name__ == "__main__":
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval, args.with_pseudo)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval, args.with_pseudo, args.ikpoint)

From a0b0e253fa314b2210597d4012d063284fd111b6 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 03:08:00 +0800
Subject: [PATCH 10/13] Add default metadata writing in molden.py

---
 .../molden => }/Multiwfn_interface/molden.py                      | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename interfaces/{ - Write [Pseudo] metadata by default in tools/molden => }/Multiwfn_interface/molden.py (100%)

diff --git a/interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py b/interfaces/Multiwfn_interface/molden.py
similarity index 100%
rename from interfaces/ - Write [Pseudo] metadata by default in tools/molden/Multiwfn_interface/molden.py
rename to interfaces/Multiwfn_interface/molden.py

From aab9ee0a3cb08a47a0823771b5f463ee12f02672 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 03:08:59 +0800
Subject: [PATCH 11/13] Update tools/molden/molden.py

Co-authored-by: SY Wang <uwsy1059@qq.com>
---
 tools/molden/molden.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index e9662b0f709..ea49642de33 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 '''
 This script is implemented referring to Molden file standard:
 https://www.theochem.ru.nl/molden/molden_format.html

From f82508e4ec34cbe691360ea42307c757028a9b54 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 04:21:08 +0800
Subject: [PATCH 12/13] Refactor moldengen to remove ikpoint parameter

Removed ikpoint parameter from moldengen function and adjusted related logic.
---
 interfaces/Multiwfn_interface/molden.py | 33 ++++++-------------------
 1 file changed, 7 insertions(+), 26 deletions(-)

diff --git a/interfaces/Multiwfn_interface/molden.py b/interfaces/Multiwfn_interface/molden.py
index 708eab87190..7149898e90e 100644
--- a/interfaces/Multiwfn_interface/molden.py
+++ b/interfaces/Multiwfn_interface/molden.py
@@ -54,34 +54,16 @@ def insert_molden_nval(text, nval):
     return text[:index] + nval + text[index:]
 
 
-def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS_Multiwfn.molden", write_cell=True, with_nval=True, with_pseudo=False, ikpoint=1):
+def moldengen(folder, ndigits=3, ngto=7, rel_r="2", fmolden="ABACUS_Multiwfn.molden", write_cell=True, with_nval=True, with_pseudo=False):
     """Call tools/molden/molden.py with Multiwfn-oriented [Nval] metadata."""
     molden = load_molden_module()
     pseudo = "pseudo" if with_pseudo else "none"
-    outputs = molden.moldengen(
-        folder,
-        ndigits,
-        ngto,
-        rel_r,
-        fmolden,
-        write_cell,
-        pseudo,
-        ikpoint,
-        return_files=True,
-        output_is_default=fmolden == "ABACUS_Multiwfn.molden",
-    )
+    out = molden.moldengen(folder, ndigits, ngto, rel_r, fmolden, write_cell, pseudo)
     if with_nval:
-        nval = make_molden_nval_section(molden, folder)
-        outputs_with_nval = []
-        for output_file, out in outputs:
-            out = insert_molden_nval(out, nval)
-            with open(output_file, "w") as file:
-                file.write(out)
-            outputs_with_nval.append((output_file, out))
-        outputs = outputs_with_nval
-    if len(outputs) == 1:
-        return outputs[0][1]
-    return [out for _, out in outputs]
+        out = insert_molden_nval(out, make_molden_nval_section(molden, folder))
+        with open(fmolden, "w") as file:
+            file.write(out)
+    return out
 
 
 
@@ -104,7 +86,6 @@ def str_to_bool(value):
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
-    parser.add_argument("--ikpoint", type=int, default=1, help="the k point to convert; use -1 to convert all k points into separate Molden files")
     parser.add_argument("--with-cell", type=str_to_bool, default=True, help="whether to write the Molden [Cell] section")
     parser.add_argument("--with-Nval", dest="with_nval", type=str_to_bool, default=True, help="whether to write the Molden [Nval] section")
     parser.add_argument("--with-pseudo", type=str_to_bool, default=False, help="whether to write the Molden [Pseudo] section")
@@ -117,4 +98,4 @@ def str_to_bool(value):
 
 if __name__ == "__main__":
     args = _argparse()
-    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval, args.with_pseudo, args.ikpoint)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, args.with_nval, args.with_pseudo)

From 47dc44a178193fcc6c20e8668515ad3c27cfa030 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Thu, 4 Jun 2026 04:28:04 +0800
Subject: [PATCH 13/13] Refactor Molden file handling functions

Refactor find_abacus_lowf_files function to remove ikpoint parameter and update related logic. Remove unused append_molden_filename_suffix and find_abacus_lowf_file_groups functions.
---
 tools/molden/molden.py | 195 ++++++++++++++---------------------------
 1 file changed, 67 insertions(+), 128 deletions(-)

diff --git a/tools/molden/molden.py b/tools/molden/molden.py
index ea49642de33..c789f201e2d 100644
--- a/tools/molden/molden.py
+++ b/tools/molden/molden.py
@@ -568,21 +568,13 @@ def read_abacus_lowf(flowf, pat=r'^(WFC_NAO_(K\d+|GAMMA\d+)(_ION\d+)?|wf((k\d+)?
 
     return nbands, nlocal, occ, band, ener, data
 
-def append_molden_filename_suffix(fmolden, suffix):
-    """Append a suffix before the Molden file extension."""
-    import os
-
-    root, ext = os.path.splitext(fmolden)
-    if ext == "":
-        ext = ".molden"
-    return f"{root}_{suffix}{ext}"
-
-def find_abacus_lowf_files(nspin, ikpoint=1, nkpoints=1):
+def find_abacus_lowf_files(nspin):
     """Find ABACUS LCAO wavefunction files in OUT.*.
 
     The converter keeps compatibility with older WFC_NAO_* names and newer
     wf*_nao.txt names. The nspin=4 spinor/SOC format is not mapped to Molden
-    because the current writer only emits restricted/Alpha/Beta MOs.
+    because the current writer only emits restricted/Alpha/Beta MOs. K-style
+    file selection uses the first k point when several k-point files are found.
     """
     import glob
     import os
@@ -592,10 +584,6 @@ def find_abacus_lowf_files(nspin, ikpoint=1, nkpoints=1):
         raise NotImplementedError("nspin=4/SOC spinor wavefunctions are not supported by this Molden writer")
     if nspin not in (1, 2):
         raise NotImplementedError(f"Unsupported nspin={nspin}")
-    if ikpoint < 1:
-        raise ValueError("ikpoint must be a positive integer when selecting one k point")
-    if ikpoint > nkpoints:
-        raise ValueError(f"ikpoint={ikpoint} is out of range for {nkpoints} k point(s)")
 
     def newest(pattern):
         def key(path):
@@ -605,12 +593,7 @@ def key(path):
         return matches[-1] if matches else None
 
     def old_k_file(spin):
-        """Select an old WFC_NAO_K* file for a spin channel and k point.
-
-        Older ABACUS files flatten spin and k into WFC_NAO_KN.txt. For nspin=2,
-        files are ordered as all k points of spin 1 followed by all k points of
-        spin 2, so the requested file index is ik + (spin - 1) * nk.
-        """
+        """Select the first-k old WFC_NAO_K* file for a spin channel."""
         files_by_k = {}
         for path in glob.glob("WFC_NAO_K*.txt"):
             match = re.match(r"^WFC_NAO_K(\d+)(?:_ION(\d+))?\.txt$", os.path.basename(path))
@@ -620,52 +603,53 @@ def old_k_file(spin):
             ion = int(match.group(2) or 0)
             if ik not in files_by_k or (ion, path) > files_by_k[ik][0]:
                 files_by_k[ik] = ((ion, path), path)
-        target = ikpoint if nspin == 1 else ikpoint + (spin - 1) * nkpoints
+        k_indices = sorted(files_by_k)
+        if not k_indices:
+            return None
+        if nspin == 1:
+            target = k_indices[0]
+        else:
+            nk = len(k_indices) // 2
+            if nk == 0:
+                return None
+            target = k_indices[0] if spin == 1 else k_indices[nk]
         return files_by_k.get(target, (None, None))[1]
 
     if nspin == 1:
         pattern_groups = [[
+            "WFC_NAO_GAMMA1.txt",
+            "WFC_NAO_GAMMA1_ION*.txt",
             "OLD_WFC_NAO_K_SPIN1",
-            f"wfk{ikpoint}_nao.txt",
-            f"WFC/wfk{ikpoint}g*_nao.txt",
+            "wf_nao.txt",
+            "wfk1_nao.txt",
+            "WFC/wfg*_nao.txt",
+            "WFC/wfk1g*_nao.txt",
         ]]
-        if ikpoint == 1:
-            pattern_groups[0] = [
-                "WFC_NAO_GAMMA1.txt",
-                "WFC_NAO_GAMMA1_ION*.txt",
-                "wf_nao.txt",
-                "WFC/wfg*_nao.txt",
-            ] + pattern_groups[0]
     else:
         pattern_groups = [
             [
-                "OLD_WFC_NAO_K_SPIN1",
-                f"wfs1k{ikpoint}_nao.txt",
-                f"wfk{ikpoint}s1_nao.txt",
-                f"WFC/wfs1k{ikpoint}g*_nao.txt",
-                f"WFC/wfk{ikpoint}s1g*_nao.txt",
-            ],
-            [
-                "OLD_WFC_NAO_K_SPIN2",
-                f"wfs2k{ikpoint}_nao.txt",
-                f"wfk{ikpoint}s2_nao.txt",
-                f"WFC/wfs2k{ikpoint}g*_nao.txt",
-                f"WFC/wfk{ikpoint}s2g*_nao.txt",
-            ],
-        ]
-        if ikpoint == 1:
-            pattern_groups[0] = [
                 "WFC_NAO_GAMMA1.txt",
                 "WFC_NAO_GAMMA1_ION*.txt",
+                "OLD_WFC_NAO_K_SPIN1",
                 "wfs1_nao.txt",
+                "wfs1k1_nao.txt",
+                "wfk1s1_nao.txt",
                 "WFC/wfs1g*_nao.txt",
-            ] + pattern_groups[0]
-            pattern_groups[1] = [
+                "WFC/wfs1k1g*_nao.txt",
+                "WFC/wfk1s1g*_nao.txt",
+            ],
+            [
                 "WFC_NAO_GAMMA2.txt",
                 "WFC_NAO_GAMMA2_ION*.txt",
+                "OLD_WFC_NAO_K_SPIN2",
                 "wfs2_nao.txt",
+                "wfs2k1_nao.txt",
+                "wfk1s2_nao.txt",
                 "WFC/wfs2g*_nao.txt",
-            ] + pattern_groups[1]
+                "WFC/wfs2k1g*_nao.txt",
+                "WFC/wfk1s2g*_nao.txt",
+            ],
+        ]
 
     files = []
     for patterns in pattern_groups:
@@ -684,16 +668,6 @@ def old_k_file(spin):
         files.append(found)
     return files
 
-def find_abacus_lowf_file_groups(nspin, ikpoint=1, nkpoints=1):
-    """Find wavefunction files grouped by selected k point."""
-    if ikpoint == -1:
-        ikpoints = range(1, nkpoints + 1)
-    elif ikpoint > 0:
-        ikpoints = [ikpoint]
-    else:
-        raise ValueError("ikpoint must be a positive integer, or -1 to convert all k points")
-    return [(ik, find_abacus_lowf_files(nspin, ik, nkpoints)) for ik in ikpoints]
-
 def write_molden_mo(coefs, mo_map, occ, ener, spin=0, ndigits=3):
     """Write Molden [MO] section from the coefficients, occupations and energies.
 
@@ -1027,34 +1001,31 @@ def read_abacus_input(finput):
     return kv
 
 def read_abacus_kpt(fkpt):
-    """Read the number of k points from an ABACUS KPT file.
-
-    The automatic mesh format is expected to look like:
+    """the way to organize information of KPT file of ABACUS still has some degree of freedom.
+    However, the only one wanted should have content like the following:
     K_POINTS
     0
     Gamma
     1 1 1 0 0 0
 
-    In this script, multiple k points are converted one k point at a time. The
-    k-point number is used only to select the corresponding ABACUS
-    wavefunction file; different k points are not merged into one Molden file.
+    in which the "Gamma" can be replaced by "MP" but the number of kpoints should be 1.
+    In the future the multiple kpoints is planned to be supported in a relatively naive way that
+    simply combining all MOs at different kpoints together but not for now. The occupation of MO
+    at different kpoints is already multiplied by the weight of kpoints, is it expected?
+
+    This function is not really read kpoints and return something, instead, it is for assert
+    the number of kpoints is 1.
     """
     with open(fkpt, 'r') as file:
         lines = file.readlines()
-    lines = [line.split("#", 1)[0].strip() for line in lines]
-    lines = [line for line in lines if line]
-    if len(lines) < 2 or lines[0] != "K_POINTS":
-        raise ValueError(f"Invalid KPT file {fkpt}.")
-    if lines[1] == "0":
-        if len(lines) < 4 or lines[2] not in ["Gamma", "MP"]:
-            raise ValueError(f"Invalid KPT file {fkpt}. Only Gamma/MP automatic mesh is supported.")
-        mesh = [int(x) for x in lines[3].split()[:3]]
-        nkpoints = mesh[0] * mesh[1] * mesh[2]
-    else:
-        nkpoints = int(lines[1])
-    if nkpoints < 1:
-        raise ValueError(f"Invalid KPT file {fkpt}. Number of k points should be positive.")
-    return nkpoints
+    lines = [line.strip() for line in lines]
+    if lines[0] == "K_POINTS":
+        if lines[1] == "0":
+            if lines[2] in ["Gamma", "MP"]:
+                if lines[3] == "1 1 1 0 0 0":
+                    return
+    raise ValueError(f"Invalid KPT file {fkpt}. Presently only 1 kpoint calculation \
+(implicit or explicit) Gamma-only calculation is supported.")
 
 def CondonShortleyPhase(index):
     """Imposing the Condon-Shortley phase on the MO index. 
@@ -1102,7 +1073,7 @@ def indexing_mo(total_gto: GTORadials, labels: list):
                 i += 2*l+1
     return out
 
-def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_cell=False, pseudo="none", ikpoint=1, return_files=False, output_is_default=False):
+def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden", write_cell=False, pseudo="none"):
     """Entrance function: generate molden file by reading the outdir of ABACUS, for only LCAO 
     calculation.
     
@@ -1111,7 +1082,6 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
         ndigits (int): the number of digits to be printed for the coefficients
         ngto (int): the number of GTOs to be fitted to the numerical atomic orbitals
         fmolden (str): the file name of the molden file
-        ikpoint (int): selected k point index; -1 converts all k points separately
     
     Returns:
         str: the content of the molden file"""
@@ -1138,7 +1108,7 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     ####################
     # write the cell   #
     ####################
-    nkpoints = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
+    _ = read_abacus_kpt(kv.get("kpoint_file", "KPT"))
     stru = read_abacus_stru(kv.get("stru_file", "STRU"))
     if write_cell:
         out += write_molden_cell(stru['lat']['const'], stru['lat']['vec'])
@@ -1189,36 +1159,19 @@ def moldengen(folder: str, ndigits=3, ngto=7, rel_r=2, fmolden="ABACUS.molden",
     out_dir = ".".join(["OUT", suffix])
     print(f"MOLDEN: Parsing ABACUS output directory {os.path.abspath(out_dir)}")
     os.chdir(out_dir)
-    mo_file_groups = find_abacus_lowf_file_groups(nspin, ikpoint, nkpoints)
-    default_fmolden = output_is_default or fmolden in (None, "ABACUS.molden")
-    if fmolden is None:
-        fmolden = "ABACUS.molden"
-    outputs = []
-    for selected_ikpoint, mo_files in mo_file_groups:
-        out_ik = out
-        print(f"MOLDEN: Converting k point {selected_ikpoint} of {nkpoints}")
-        for ispin, mo_file in enumerate(mo_files):
-            print(f"MOLDEN: Parsing wavefunction file {os.path.abspath(mo_file)}")
-            nbands, nlocal, occ, band, ener, data_np = read_abacus_lowf(mo_file)
-            coefs = data_np.tolist()
-            out_ik += write_molden_mo(coefs, mo_map, occ, ener, spin=ispin, ndigits=ndigits)
-        fmolden_ik = fmolden
-        if len(mo_file_groups) > 1 or (nkpoints > 1 and default_fmolden):
-            fmolden_ik = append_molden_filename_suffix(fmolden, f"k{selected_ikpoint}")
-        outputs.append((fmolden_ik, out_ik))
+    mo_files = find_abacus_lowf_files(nspin)
+    for ispin, mo_file in enumerate(mo_files):
+        print(f"MOLDEN: Parsing wavefunction file {os.path.abspath(mo_file)}")
+        nbands, nlocal, occ, band, ener, data_np = read_abacus_lowf(mo_file)
+        coefs = data_np.tolist()
+        out += write_molden_mo(coefs, mo_map, occ, ener, spin=ispin, ndigits=ndigits)
+
     os.chdir(cwd)
-    written = []
-    for fmolden_ik, out_ik in outputs:
-        fmolden_ik = os.path.abspath(fmolden_ik)
-        print(f"MOLDEN: Writing Molden file {fmolden_ik}")
-        with open(fmolden_ik, "w") as file:
-            file.write(out_ik)
-        written.append((fmolden_ik, out_ik))
-    if return_files:
-        return written
-    if len(written) == 1:
-        return written[0][1]
-    return [out_ik for _, out_ik in written]
+    fmolden = os.path.abspath(fmolden)
+    print(f"MOLDEN: Writing Molden file {fmolden}")
+    with open(fmolden, "w") as file:
+        file.write(out)
+    return out
 
 import unittest
 class TestNAO2GTO(unittest.TestCase):
@@ -1558,7 +1511,6 @@ def _argparse():
     -n, --ndigits: the number of digits for the MO coefficients. For MO coefficients smaller than 10^-n, they will be set to 0.
     -g, --ngto: the number of GTOs to fit ABACUS NAOs.
     -r, --rel_r: the relative cutoff radius for the GTOs.
-    --ikpoint: the k point to convert. Use -1 to convert all k points into separate Molden files.
     --with-cell: write the Molden [Cell] section.
     --with-pseudo: write the Molden [Pseudo] section from UPF z_valence.
     -o, --output: the output Molden file name.
@@ -1577,7 +1529,6 @@ def _argparse():
     parser.add_argument("-n", "--ndigits", type=int, default=3, help="the number of digits for the MO coefficients")
     parser.add_argument("-g", "--ngto", type=int, default=7, help="the number of GTOs to fit ABACUS NAOs")
     parser.add_argument("-r", "--rel_r", type=str, default="2", help="the relative cutoff radius for the GTOs; comma-separated values enable multi-start fitting")
-    parser.add_argument("--ikpoint", type=int, default=1, help="the k point to convert; use -1 to convert all k points into separate Molden files")
     parser.add_argument("--with-cell", action="store_true", help="write the Molden [Cell] section")
     parser.add_argument("--with-pseudo", action="store_true", help="write the Molden [Pseudo] section from UPF z_valence")
     parser.add_argument("-o", "--output", type=str, default="ABACUS.molden", help="the output Molden file name")
@@ -1588,21 +1539,9 @@ def _argparse():
     #unittest.main(exit=False)
     args = _argparse()
     pseudo = "pseudo" if args.with_pseudo else "none"
-    generated = moldengen(
-        args.folder,
-        args.ndigits,
-        args.ngto,
-        args.rel_r,
-        args.output,
-        args.with_cell,
-        pseudo,
-        args.ikpoint,
-        return_files=True,
-        output_is_default=args.output == "ABACUS.molden",
-    )
-    generated_files = ", ".join(path for path, _ in generated)
+    moldengen(args.folder, args.ndigits, args.ngto, args.rel_r, args.output, args.with_cell, pseudo)
     print(" ".join("*"*10).center(80, " "))
-    print(f"""MOLDEN: Generated Molden file(s) {generated_files} from ABACUS calculation in folder {args.folder}.
+    print(f"""MOLDEN: Generated Molden file {args.output} from ABACUS calculation in folder {args.folder}.
 WARNING: use at your own risk because the NAO2GTO will not always conserve the shape of radial function, therefore
 the total number of electrons may not be conserved. Always use after a re-normalization operation.""")
     citation = """If you use this script in your research, please cite the following paper:\n