Remove hessian conversion on normal mode calculation (#243)

* Ensure no unit conversion in hessian

* Use a copy for value and valuearray

* Raise exception on value `to` with inplace

* Add missed changelog updates

Author: t-young31

autode/hessians.py

@@ -234,9 +234,9 @@ def _mass_weighted(self) -> np.ndarray:
             axis=np.newaxis,
         )
 
-        return Hessian(
-            H / np.sqrt(np.outer(mass_array, mass_array)), units="J m^-2 kg^-1"
-        )
+        return np.array(
+            H / np.sqrt(np.outer(mass_array, mass_array))
+        )  # J Å^-2 kg^-1
 
     @cached_property
     def _proj_mass_weighted(self) -> np.ndarray:

autode/units.py

@@ -211,15 +211,15 @@ def energy_unit_from_name(name: str) -> "autode.units.Unit":
 
 
 ha_per_ang = CompositeUnit(
-    ha, per=[ang], aliases=["ha Å-1", "ha Å^-1", "ha/ang"]
+    ha, per=[ang], aliases=["ha / Å", "ha Å-1", "ha Å^-1", "ha/ang"]
 )
 
 ha_per_a0 = CompositeUnit(
-    ha, per=[a0], aliases=["ha a0-1", "ha a0^-1", "ha/bohr"]
+    ha, per=[a0], aliases=["ha / a0", "ha a0-1", "ha a0^-1", "ha/bohr"]
 )
 
 ev_per_ang = CompositeUnit(
-    ev, per=[ang], aliases=["ha a0-1", "ev Å^-1", "ev/ang"]
+    ev, per=[ang], aliases=["ev / Å", "ev Å^-1", "ev/ang"]
 )
 
 kcalmol_per_ang = CompositeUnit(

autode/values.py

@@ -37,9 +37,12 @@
 )
 
 
-def _to(value: Union["Value", "ValueArray"], units: Union[Unit, str]):
+def _to(
+    value: Union["Value", "ValueArray"], units: Union[Unit, str], inplace: bool
+) -> Any:
     """
-    Convert a value or value array to a new unit and return a copy
+    Convert a value or value array to a new unit and return a copy if
+    inplace=False
 
     ---------------------------------------------------------------------------
     Arguments:
@@ -65,23 +68,26 @@ def _to(value: Union["Value", "ValueArray"], units: Union[Unit, str]):
             f"No viable unit conversion from {value.units} -> {units}"
         )
 
-    # Convert to the base unit, then to the new units
-    c = float(units.conversion / value.units.conversion)
-
-    if isinstance(value, Value):
-        return value.__class__(float(value) * c, units=units)
-
-    elif isinstance(value, ValueArray):
-        value[:] = np.array(value, copy=True) * c
-        value.units = units
-        return value
-
-    else:
+    if not (isinstance(value, Value) or isinstance(value, ValueArray)):
         raise ValueError(
             f"Cannot convert {value} to new units. Must be one of"
             f" Value of ValueArray"
         )
 
+    if isinstance(value, Value) and inplace:
+        raise ValueError(
+            "Cannot modify a value inplace as floats are immutable"
+        )
+
+    # Convert to the base unit, then to the new units
+    c = float(units.conversion / value.units.conversion)
+
+    new_value = value if inplace else value.copy()
+    new_value *= c
+    new_value.units = units
+
+    return None if inplace else new_value
+
 
 def _units_init(value, units: Union[Unit, str, None]):
     """Initialise the units of this value
@@ -171,6 +177,11 @@ def _other_same_units(self, other):
 
         return other.to(self.units)
 
+    def _like_self_from_float(self, value: float) -> "Value":
+        new_value = self.__class__(value, units=self.units)
+        new_value.__dict__.update(self.__dict__)
+        return new_value
+
     def __eq__(self, other) -> bool:
         """Equality of two values, which may be in different units"""
 
@@ -210,8 +221,8 @@ def __add__(self, other) -> "Value":
         if isinstance(other, np.ndarray):
             return other + float(self)
 
-        return self.__class__(
-            float(self) + self._other_same_units(other), units=self.units
+        return self._like_self_from_float(
+            float(self) + self._other_same_units(other)
         )
 
     def __mul__(self, other) -> Union[float, "Value"]:
@@ -225,8 +236,8 @@ def __mul__(self, other) -> Union[float, "Value"]:
             )
             return float(self) * self._other_same_units(other)
 
-        return self.__class__(
-            float(self) * self._other_same_units(other), units=self.units
+        return self._like_self_from_float(
+            float(self) * self._other_same_units(other)
         )
 
     def __rmul__(self, other) -> Union[float, "Value"]:
@@ -240,17 +251,11 @@ def __sub__(self, other) -> "Value":
 
     def __floordiv__(self, other) -> Union[float, "Value"]:
         x = float(self) // self._other_same_units(other)
-        if isinstance(other, Value):
-            return x
-
-        return self.__class__(x, units=self.units)
+        return x if isinstance(other, Value) else self._like_self_from_float(x)
 
     def __truediv__(self, other) -> Union[float, "Value"]:
         x = float(self) / self._other_same_units(other)
-        if isinstance(other, Value):
-            return x
-
-        return self.__class__(x, units=self.units)
+        return x if isinstance(other, Value) else self._like_self_from_float(x)
 
     def __abs__(self) -> "Value":
         """Absolute value"""
@@ -269,7 +274,7 @@ def to(self, units):
         Raises:
             (TypeError):
         """
-        return _to(self, units)
+        return _to(self, units, inplace=False)
 
 
 class Energy(Value):
@@ -652,7 +657,21 @@ def to(self, units) -> Any:
         Raises:
             (TypeError):
         """
-        return _to(self, units)
+        return _to(self, units, inplace=False)
+
+    def to_(self, units) -> None:
+        """
+        Convert this array into a set of new units, inplace. This will not copy
+        the array
+
+        -----------------------------------------------------------------------
+        Returns:
+            (None)
+
+        Raises:
+            (TypeError):
+        """
+        _to(self, units, inplace=True)
 
     def __array_finalize__(self, obj):
         """See https://numpy.org/doc/stable/user/basics.subclassing.html"""

doc/changelog.rst

@@ -8,19 +8,21 @@ Changelog
 
 Usability improvements/Changes
 ******************************
-*
+- :code:`autode.value.ValueArray.to()` now defaults to copying the object rather than inplace modification
 
 
 Functionality improvements
 **************************
--
+- Adds a :code:`to_` method to :code:`autode.value.ValueArray` for explicit inplace modification of the array
 
 
 Bug Fixes
 *********
 - Fixes :code:`ERROR` logging level being ignored from environment variable :code:`AUTODE_LOG_LEVEL`
 - Fixes :code:`autode.values.Value` instances generating items with units on division, and throw a warning if multiplying
 - Fixes the printing of cartesian constraints in XTB input files, meaning they are no longer ignored
+- Fixes :code:`Hessian` instances changing units when normal modes are calculated
+- Fixes an incorrect alias for :code:`ev_per_ang`
 
 
 1.3.4
@@ -32,7 +34,6 @@ Feature additions.
 Usability improvements/Changes
 ******************************
 * Throw useful exception for invalid :code:`ade.Config.ts_template_folder_path`
-* Adds the reaction temperature to the unique reaction hash
 
 
 Functionality improvements
@@ -44,7 +45,7 @@ Functionality improvements
 
 Bug Fixes
 *********
--
+- Fixes calculation :code:`clean_up()` failing with a null filename
 
 
 1.3.3

tests/test_hessian.py

@@ -13,7 +13,7 @@
 from autode.species import Molecule
 from autode.values import Frequency
 from autode.geom import calc_rmsd
-from autode.units import wavenumber
+from autode.units import wavenumber, ha_per_ang_sq
 from autode.exceptions import CalculationException
 from autode.wrappers.keywords import pbe0
 from autode.transition_states.base import displaced_species_along_mode
@@ -243,6 +243,7 @@ def assert_correct_co2_frequencies(hessian, expected=(666, 1415, 2517)):
     """Ensure the projected frequencies of CO2 are roughly right"""
     nu_1, nu_2, nu_3 = expected
 
+    print(hessian.frequencies_proj)
     assert sum(freq == 0.0 for freq in hessian.frequencies_proj) == 5
 
     # Should have a degenerate bending mode for CO2 with ν = 666 cm-1
@@ -419,6 +420,7 @@ def test_hessian_modes():
 
     h2o = Molecule("H2O_hess_orca.xyz")
     h2o.hessian = h2o_hessian_arr
+    assert h2o.hessian.units == ha_per_ang_sq
 
     # The structure is a minimum, thus there should be no imaginary frequencies
     assert h2o.imaginary_frequencies is None
@@ -441,6 +443,9 @@ def test_hessian_modes():
             vib_mode, h2o.hessian.normal_modes[6 + i], atol=0.1
         ) or np.allclose(vib_mode, -h2o.hessian.normal_modes[6 + i], atol=0.1)
 
+    # Hessian units should be retained
+    assert h2o.hessian.units == ha_per_ang_sq
+
 
 @testutils.work_in_zipped_dir(os.path.join(here, "data", "hessians.zip"))
 def test_proj_modes():
@@ -595,6 +600,8 @@ def test_nwchem_hessian_co2():
         keywords=ade.HessianKeywords(),
     )
     calc.set_output_filename("CO2_hess_nwchem.out")
+    print(co2.hessian)
+    print(co2.hessian._mass_weighted)
     assert_correct_co2_frequencies(
         hessian=co2.hessian, expected=(659.76, 1406.83, 2495.73)
     )

tests/test_value.py

@@ -3,6 +3,7 @@
 from autode.constants import Constants
 from autode.units import ha, kjmol, kcalmol, ev, ang, a0, nm, pm, m, rad, deg
 from autode.values import (
+    _to,
     Value,
     Distance,
     MWDistance,
@@ -278,3 +279,25 @@ def test_div_mul_generate_floats():
 
     # Note: this behaviour is not ideal. But it is better than having the wrong units
     assert isinstance(e * e, float)
+
+
+def test_operations_maintain_other_attrs():
+
+    e = Energy(1, estimated=True, units="eV")
+    assert e.is_estimated and e.units == ev
+
+    e *= 2
+    assert e.is_estimated and e.units == ev
+
+    e /= 2
+    assert e.is_estimated and e.units == ev
+
+    a = e * 2
+    assert a.is_estimated and e.units == ev
+
+
+def test_inplace_value_modification_raises():
+
+    e = Energy(1, units="Ha")
+    with pytest.raises(ValueError):  # floats are immutable
+        _to(e, units="eV", inplace=True)

tests/test_values.py

@@ -106,4 +106,22 @@ class InvalidValue(float):
 def test_to_unsupported():
 
     with pytest.raises(ValueError):
-        _ = _to(InvalidValue(), Unit())
+        _ = _to(InvalidValue(), Unit(), inplace=True)
+
+
+def test_inplace_modification():
+
+    x = Gradient([[1.0, 1.0, 1.0]], units="Ha / Å")
+    return_value = x.to_("eV / Å")
+    assert return_value is None
+
+    assert not np.allclose(x, np.ones(shape=(1, 3)))
+
+
+def test_copy_conversion():
+
+    x = Gradient([[1.0, 1.0, 1.0]], units="Ha / Å")
+    y = x.to("eV / Å")
+
+    assert not np.allclose(x, y)
+    assert np.allclose(x, np.ones(shape=(1, 3)))