193 read the s2 operator expectation value and spin contamination for adf calculations

src/tcutility/results/adf.py

@@ -179,6 +179,9 @@ def get_properties(info: Result) -> Result:
             - **vibrations.character (str)** – Character of the molecule based on the number of imaginary vibrational modes. Can be "minimum" or "transition state".
             - **vdd.charges (nparray[float] (1D))** - 1D array of Voronoi Deformation Denisty (VDD) charges in [electrons], being the difference between the final (SCF) and initial VDD charges.
             - **vdd.charges.{symmetry label} (nparray[float] (1D))** - 1D array of Voronoi Deformation Denisty (VDD) charges in [electrons] per irrep.
+            - **s2** - expectation value of the :math:`S^2` operator.
+            - **s2_expected** - ideal expectation value of the :math:`S^2` operator. For restricted calculations this should always equal ``s2``.
+            - **spin_contamination** - the amount of spin-contamination observed in this calculation. It is equal to (s2 - s2_expected) / (s2_expected). Ideally this value should be below 0.1.
     """
 
     assert info.engine == "adf", f"This function reads ADF data, not {info.engine} data"
@@ -218,6 +221,23 @@ def get_properties(info: Result) -> Result:
     except KeyError:
         pass
 
+    # read spin-squared operator info
+    # the total spin
+    S = info.adf.spin_polarization * 1 / 2
+    ret.s2_expected = S * (S + 1)
+    # this is the real expectation value
+    if ('Properties', 'S2calc') in reader_adf:
+        ret.s2 = reader_adf.read('Properties', 'S2calc')
+    else:
+        ret.s2 = 0
+
+    # calculate the contamination
+    # if S is 0 then we will get a divide by zero error, but spin-contamination should be 0
+    if S != 0:
+        ret.spin_contamination = (ret.s2 - ret.s2_expected) / (ret.s2_expected)
+    else:
+        ret.spin_contamination = 0
+
     return ret
 
 

test/test_adf_properties.py

@@ -9,8 +9,29 @@ def ethane_res() -> results.Result:
     return results.read(pl.Path(__file__).parent / "fixtures" / "ethane_adf")
 
 
+@pytest.fixture()
+def unrestricted_res() -> results.Result:
+    return results.read(pl.Path(__file__).parent / "fixtures" / "radical_addition_ts/geometry")
+
+
 def test_energy_terms(ethane_res):
     print(ethane_res.keys())
     assert ethane_res.properties.energy.gibbs == pytest.approx(-522.38, abs=1e-2)
     assert ethane_res.properties.energy.enthalpy == pytest.approx(-507.64, abs=1e-2)
     assert ethane_res.properties.energy.nuclear_internal == pytest.approx(29.09, abs=1e-2)
+
+
+def test_s2_restricted(ethane_res):
+    assert ethane_res.properties.s2 == 0
+    assert ethane_res.properties.s2_expected == 0
+    assert ethane_res.properties.spin_contamination == 0
+
+
+def test_s2_spinpol1(unrestricted_res):
+    assert round(unrestricted_res.properties.s2, 5) == 0.75845
+    assert unrestricted_res.properties.s2_expected == 0.75
+    assert round(unrestricted_res.properties.spin_contamination, 4) == 0.0113
+
+
+if __name__ == '__main__':
+    pytest.main()