Improve phonon code.

matcalc/phonon.py

@@ -23,7 +23,7 @@ def __init__(
         calculator: Calculator,
         atom_disp=0.015,
         supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)),
-        fmax=0.1,
+        fmax=0.01,
         relax_structure=True,
         t_step=10,
         t_max=1000,
@@ -34,7 +34,7 @@ def __init__(
             calculator: ASE Calculator to use.
             atom_disp: Atomic displacement
             supercell_matrix: Supercell matrix to use. Defaults to 2x2x2 supercell.
-            fmax: Max forces.
+            fmax: Max forces. This criterion is more stringent than for simple relaxation.
             relax_structure: Whether to first relax the structure. Set to False if structures provided are pre-relaxed
                 with the same calculator.
             t_step: Temperature step.

matcalc/relaxation.py

@@ -104,7 +104,7 @@ def __init__(
             interval (int): the step interval for saving the trajectories.
         """
         self.calculator = calculator
-        self.opt_class: Optimizer = OPTIMIZERS[optimizer] if isinstance(optimizer, str) else optimizer
+        self.optimizer: Optimizer = OPTIMIZERS[optimizer] if isinstance(optimizer, str) else optimizer
         self.fmax = fmax
         self.interval = interval
         self.steps = steps
@@ -133,14 +133,15 @@ def calc(self, structure) -> dict:
         atoms.set_calculator(self.calculator)
         stream = io.StringIO()
         with contextlib.redirect_stdout(stream):
-            obs = TrajectoryObserver(atoms)
             atoms = ExpCellFilter(atoms)
-            optimizer = self.opt_class(atoms)
-            optimizer.attach(obs, interval=self.interval)
+            optimizer = self.optimizer(atoms)
+            if self.traj_file is not None:
+                obs = TrajectoryObserver(atoms)
+                optimizer.attach(obs, interval=self.interval)
             optimizer.run(fmax=self.fmax, steps=self.steps)
-            obs()
-        if self.traj_file is not None:
-            obs.save(self.traj_file)
+            if self.traj_file is not None:
+                obs()
+                obs.save(self.traj_file)
         atoms = atoms.atoms
 
         final_structure = ase_adaptor.get_structure(atoms)

tests/test_phonon.py

@@ -6,17 +6,19 @@
 from matcalc.phonon import PhononCalc
 
 
-def test_PhononCalc(Li2O, M3GNetUPCalc):
+def test_PhononCalc(Li2O, LiFePO4, M3GNetUPCalc):
     """Tests for PhononCalc class"""
     calculator = M3GNetUPCalc
-    pcalc = PhononCalc(calculator, supercell_matrix=((3, 0, 0), (0, 3, 0), (0, 0, 3)), fmax=0.001, t_step=2, t_max=1500)
+    # Note that the fmax is probably too high. This is for testing purposes only.
+    pcalc = PhononCalc(calculator, supercell_matrix=((2, 0, 0), (0, 2, 0), (0, 0, 2)), fmax=0.1, t_step=50, t_max=1000)
     results = pcalc.calc(Li2O)
 
     # Test values at 100 K
-    assert results["heat_capacity"][100] == pytest.approx(46.686590135216676, abs=0.01)
-    assert results["entropy"][100] == pytest.approx(29.347199723052057, abs=0.01)
-    assert results["free_energy"][100] == pytest.approx(16.273094775993236, abs=0.01)
+    ind = results["temperatures"].tolist().index(300)
+    assert results["heat_capacity"][ind] == pytest.approx(59.918928933451305, abs=0.01)
+    assert results["entropy"][ind] == pytest.approx(51.9081928335805, abs=0.01)
+    assert results["free_energy"][ind] == pytest.approx(11.892105644441045, abs=0.01)
 
-    results = list(pcalc.calc_many([Li2O] * 2))
+    results = list(pcalc.calc_many([Li2O, LiFePO4]))
     assert len(results) == 2
-    assert results[-1]["heat_capacity"][100] == pytest.approx(46.68654747038025, abs=0.01)
+    assert results[-1]["heat_capacity"][ind] == pytest.approx(550.6419940551511, abs=0.01)