BUG: issue with running NPT MD

[fraricci]

Hi there.

I've been playing around with some MD, in particular with NPT ensemble
and I think I found two issue/bugs.

1. The "Schur" decomposition to transform the lattice into a upper triangular form does not work properly. The output cell is completely changed from the original

si_structure = Structure(
    lattice=[[0, 2.73, 2.73], [2.73, 0, 2.73], [2.73, 2.73, 0]],
    species=["Si", "Si"],
    coords=[[0, 0, 0], [0.25, 0.25, 0.25]],
)


atoms = si_structure.to_ase_atoms()
schur_decomp, _ = schur(atoms.get_cell(complete=True), output="complex")
atoms.set_cell(schur_decomp.real, scale_atoms=True)
st_schur = AseAtomsAdaptor.get_structure(
    atoms,
    cls=Structure)
print(si_structure.lattice.abc)
print(st_schur.lattice.abc)

: (3.8608030252785492, 3.8608030252785492, 3.8608030252785492)
: (2.73, 5.460000000000001, 2.73)

An alternative approach that works is to use the ASE cell.standard_form() that uses a QR decomposition:

atoms = si_structure.to_ase_atoms()
print(atoms.cell.lengths())
print(atoms.cell[:])

atoms.set_cell(atoms.cell.standard_form("upper")[0])
print(atoms.cell[:])
print(atoms.cell.lengths())

: [3.86080303 3.86080303 3.86080303]
: [[0.   2.73 2.73]
:  [2.73 0.   2.73]
:  [2.73 2.73 0.  ]]
: [[ 3.15233247  1.11451783  1.93040151]
:  [-0.          3.3435535   1.93040151]
:  [-0.          0.          3.86080303]]
: [3.86080303 3.86080303 3.86080303]

2. the "berendsen" dynamcs is listed among the possible choices, but the default ase_md_kwargs include the externalstress which is only for NPT ase function, but not for NPTBerendsen, which actually require the compressibilty instead. See

   atomate2/src/atomate2/ase/md.py

   Line 244 in 3672ac4

   self.ase_md_kwargs["externalstress"] = self.p_schedule[0] * 1e3 * units.bar

Please, let me know if this makes sense, and I could provide a fix in case.

Thanks!

[JaGeo]

Thanks @fraricci ! I would be happy about a fix but do not have enough knowledge about the modules to confirm the bugs.

@esoteric-ephemera could you maybe have a look? Or refer it to someone who can?

[esoteric-ephemera]

Hey @fraricci probably best to use the QR decomposition from ASE - the direct coordinates aren't being updated by ASE after the decomposition, so the structural symmetry is broken

For the default MD kwargs, some logic around not over-setting what the user requests is also needed. I can implement this if you'd like, you're also welcome to submit a PR to get this started!

[fraricci]

Looks good to me, thanks!