Bug risk: If a calculator re-orders atom indices, run_calc will raise an error if a geom_file is provided

[Andrew-S-Rosen]

Details about the quacc environment

• quacc version: 0.6.0
• Python version: 3.9+

What is the issue?

In the quacc.runners.ase.run_calc function, there is a section that updates the Atoms object based on a provided geometry file as a way to ensure the output is updated appropriately when an internal relaxation is carried out. While some calculators (e.g. Vasp at the time of writing) do this automatically in ASE, others do not. So, the method ensures that all calculators update Atoms with the new positions, cell, and properties.

However, there is a bug risk here. If ASE rewrites the order of the atoms when writing the input file or if the calculator itself reorders atoms during the simulation, then the input Atoms object's atom ordering will no longer be the same ordering as the output geometry file. This will result in a mismatch, and a ValueError (in most cases) will be raised.

This is not an issue for any of the current recipes, but it could be an issue for future ones. For instance, CASTEP reorders atoms, the D3 calculator reorders atoms, and VASP can too. For VASP, it is not currently an issue because we are relying on the fact that the calculator always returns the updated Atoms object, but there is no guarantee that this behavior will stay the same in the future. In fact, it most certainly won't if/when the new GenericFileIOCalculator is implemented.

It is not immediately clear to me the best mechanism to handle this in a general way. Right now, the logic is as follows: take the Atoms object resulting from Atoms.get_potential_energy() (which has the attached calculator and its properties), read in a geometry file to get the new Atoms object, and then patch the positions and cell in the original Atoms object to match the updated structure. This will break if there has been reordering. One could switch this around and read in the new Atoms object and then patch that with the calculator (and .info) from the original Atoms object, but that doesn't solve anything because then the atom-wise calculator properties (e.g. forces, magmoms) would be mismatched on reordering.

It is possible that there may be an ASE-certified mechanism to do this in the future (see here), in which case we can defer to that instead of trying to make atoms.get_potential_energy() consistent.