Adding Abinit magmoms from netCDF files to Structure.site_properties (#…

…3936)

* Added magmoms reading from nc files with Abinit.

* Added tests.

* Adding GSR files to tests setup.

* Fixing test for collinear magmom.

* Fixing test for noncollinear magmom.

* Revert backward incompatible change introduced in #2a3608f and #565a8b4.

* Added a test for PAW pseudopotentials with abinit.

* Fix pseudo test.

* Compressed large netCDF files for abinit.

* Compressed pseudo file for abinit.

src/pymatgen/io/abinit/netcdf.py

@@ -15,6 +15,7 @@
 from pymatgen.core.structure import Structure
 from pymatgen.core.units import ArrayWithUnit
 from pymatgen.core.xcfunc import XcFunc
+from pymatgen.electronic_structure.core import Magmom
 
 if TYPE_CHECKING:
     from typing_extensions import Self
@@ -315,6 +316,23 @@ def structure_from_ncdata(ncdata, site_properties=None, cls=Structure):
     # type_atom[:natom] --> index Between 1 and number of atom species
     type_atom = ncdata.read_value("atom_species")
 
+    try:
+        intgden = ncdata.read_value("intgden")
+        nspden = intgden.shape[1]
+    except NetcdfReaderError:
+        intgden = None
+        nspden = None
+
+    if intgden is not None:
+        if nspden == 2:
+            magmoms = Magmom(intgden[:, 1] - intgden[:, 0])
+        elif nspden == 4:
+            magmoms = [Magmom([intg_at[1], intg_at[2], intg_at[3]]) for intg_at in intgden]
+        else:
+            magmoms = None
+    else:
+        magmoms = None
+
     # Fortran to C index and float --> int conversion.
     species = n_atom * [None]
     for atom in range(n_atom):
@@ -326,6 +344,9 @@ def structure_from_ncdata(ncdata, site_properties=None, cls=Structure):
         for prop in site_properties:
             dct[prop] = ncdata.read_value(prop)
 
+    if magmoms is not None and "magmom" not in dct:
+        dct["magmom"] = magmoms
+
     structure = cls(lattice, species, red_coords, site_properties=dct)
 
     # Quick and dirty hack.

src/pymatgen/io/abinit/pseudos.py

@@ -1078,8 +1078,8 @@ def read_ppdesc(self, filename):
         # Assume file with the abinit header.
         lines = _read_nlines(filename, 80)
 
-        for lineno, line in enumerate(lines, start=1):
-            if lineno == 3:
+        for lineno, line in enumerate(lines):
+            if lineno == 2:
                 try:
                     tokens = line.split()
                     pspcod, _pspxc = map(int, tokens[:2])
@@ -1095,7 +1095,7 @@ def read_ppdesc(self, filename):
 
                 if pspcod == 7:
                     # PAW -> need to know the format pspfmt
-                    tokens = line.split()
+                    tokens = lines[lineno + 1].split()
                     pspfmt, _creatorID = tokens[:2]
 
                     ppdesc = ppdesc._replace(format=pspfmt)

tests/io/abinit/test_netcdf.py

@@ -1,7 +1,11 @@
 from __future__ import annotations
 
+import os
+import tarfile
+
 import numpy as np
 import pytest
+from monty.tempfile import ScratchDir
 from numpy.testing import assert_allclose, assert_array_equal
 from pymatgen.core.structure import Structure
 from pymatgen.io.abinit import EtsfReader
@@ -72,12 +76,34 @@ def test_read_si2(self):
             # Initialize pymatgen structure from GSR.
             structure = data.read_structure()
             assert isinstance(structure, Structure)
+            assert "magmom" not in structure.site_properties
 
             # Read ixc.
             # TODO: Upgrade GSR file.
             # xc = data.read_abinit_xcfunc()
             # assert xc == "LDA"
 
+    @pytest.mark.skipif(netCDF4 is None, reason="Requires Netcdf4")
+    def test_read_fe(self):
+        with ScratchDir(".") as tmp_dir:
+            with tarfile.open(f"{TEST_DIR}/Fe_magmoms_collinear_GSR.tar.xz", mode="r:xz") as t:
+                t.extractall(tmp_dir)
+                ref_magmom_collinear = [-0.5069359730980665]
+                path = os.path.join(tmp_dir, "Fe_magmoms_collinear_GSR.nc")
+
+                with EtsfReader(path) as data:
+                    structure = data.read_structure()
+                    assert structure.site_properties["magmom"] == ref_magmom_collinear
+
+            with tarfile.open(f"{TEST_DIR}/Fe_magmoms_noncollinear_GSR.tar.xz", mode="r:xz") as t:
+                t.extractall(tmp_dir)
+                ref_magmom_noncollinear = [[0.357939487, 0.357939487, 0]]
+                path = os.path.join(tmp_dir, "Fe_magmoms_noncollinear_GSR.nc")
+
+                with EtsfReader(path) as data:
+                    structure = data.read_structure()
+                    assert structure.site_properties["magmom"] == ref_magmom_noncollinear
+
 
 class TestAbinitHeader(PymatgenTest):
     def test_api(self):

tests/io/abinit/test_pseudos.py

@@ -1,9 +1,11 @@
 from __future__ import annotations
 
 import os.path
+import tarfile
 from collections import defaultdict
 
 import pytest
+from monty.tempfile import ScratchDir
 from pymatgen.io.abinit.pseudos import Pseudo, PseudoTable
 from pymatgen.util.testing import TEST_FILES_DIR, PymatgenTest
 from pytest import approx
@@ -17,6 +19,7 @@ def setUp(self):
         nc_pseudo_fnames["Si"] = [f"{TEST_DIR}/{file}" for file in ("14si.pspnc", "14si.4.hgh", "14-Si.LDA.fhi")]
 
         self.nc_pseudos = defaultdict(list)
+        self.paw_pseudos = defaultdict(list)
 
         for symbol, file_names in nc_pseudo_fnames.items():
             for file_name in file_names:
@@ -90,6 +93,35 @@ def test_nc_pseudos(self):
         # Test pickle
         self.serialize_with_pickle(table, test_eq=False)
 
+    def test_paw_pseudos(self):
+        """Test 28ni.paw."""
+        file_name = f"{TEST_DIR}/28ni.paw.tar.xz"
+        symbol = "Ni"
+        with ScratchDir(".") as tmp_dir, tarfile.open(file_name, mode="r:xz") as t:
+            t.extractall(tmp_dir)
+            path = os.path.join(tmp_dir, "28ni.paw")
+            pseudo = Pseudo.from_file(path)
+
+            assert repr(pseudo)
+            assert str(pseudo)
+            assert not pseudo.isnc
+            assert pseudo.ispaw
+            assert pseudo.Z == 28
+            assert pseudo.symbol == symbol
+            assert pseudo.Z_val == 18
+            assert pseudo.paw_radius >= 0.0
+
+            assert pseudo.l_max == 2
+            assert pseudo.l_local == 0
+            assert pseudo.supports_soc
+            assert pseudo.md5 is not None
+
+            # Test pickle
+            self.serialize_with_pickle(pseudo)
+
+            # Test MSONable
+            self.assert_msonable(pseudo)
+
     def test_pawxml_pseudos(self):
         """Test O.GGA_PBE-JTH-paw.xml."""
         oxygen = Pseudo.from_file(f"{TEST_DIR}/O.GGA_PBE-JTH-paw.xml")