Fix nelectrons not updating when replacing species in Molecule (#3269)

* make property Molecule.nelectrons dynamic

* add test to check nelectrons is updated when replacing species

Author: janosh

pymatgen/core/structure.py

@@ -2902,22 +2902,16 @@ def __init__(
             raise StructureError("Molecule contains sites that are less than 0.01 Angstrom apart!")
 
         self._charge = charge
-        nelectrons = 0.0
-        for site in sites:
-            for sp, amt in site.species.items():
-                if not isinstance(sp, DummySpecies):
-                    nelectrons += sp.Z * amt
-        nelectrons -= charge
-        self._nelectrons = nelectrons
+        n_electrons = self.nelectrons
         if spin_multiplicity:
-            if charge_spin_check and (nelectrons + spin_multiplicity) % 2 != 1:
+            if charge_spin_check and (n_electrons + spin_multiplicity) % 2 != 1:
                 raise ValueError(
                     f"Charge of {self._charge} and spin multiplicity of {spin_multiplicity} "
                     "is not possible for this molecule!"
                 )
             self._spin_multiplicity = spin_multiplicity
         else:
-            self._spin_multiplicity = 1 if nelectrons % 2 == 0 else 2
+            self._spin_multiplicity = 1 if n_electrons % 2 == 0 else 2
 
     @property
     def charge(self) -> float:
@@ -2932,7 +2926,13 @@ def spin_multiplicity(self) -> float:
     @property
     def nelectrons(self) -> float:
         """Number of electrons in the molecule."""
-        return self._nelectrons
+        n_electrons = 0.0
+        for site in self:
+            for sp, amt in site.species.items():
+                if not isinstance(sp, DummySpecies):
+                    n_electrons += sp.Z * amt
+        n_electrons -= self.charge
+        return n_electrons
 
     @property
     def center_of_mass(self) -> np.ndarray:

tests/core/test_structure.py

@@ -1779,6 +1779,10 @@ def test_prop(self):
         mol = Molecule(["C", "H", "H", "H", "H"], self.coords, charge=1)
         assert mol.spin_multiplicity == 2
         assert mol.nelectrons == 9
+        # https://github.com/materialsproject/pymatgen/issues/3265
+        # replace species and ensure nelectrons is updated
+        mol[0] = "N"
+        assert mol.nelectrons == 10
 
         # Triplet O2
         mol = IMolecule(["O"] * 2, [[0, 0, 0], [0, 0, 1.2]], spin_multiplicity=3)