ASE: Shape mismatch when using FrechetCellFilter filter

Question

I want to use a FrechetCellFilter into a MD simulation with ASE to make sure both the atomic forces and unit cell are simultaneously minimised.

I encounter a shape mismatch error: ValueError: shape mismatch: value array of shape (11,3) could not be broadcast to indexing result of shape (8,3)

I perform first a QuasiNewton optimisation which goes through correctly and then a NPT relaxation, where the code complains about the filter( len(ecf)=11 ) being broadcasted to the atoms (len(atoms)=8).

Why is the FrechetCellFilter working for the optimization and not for the NPT run? Please find below the xyz input file, the python script and the output produced.

Cheers

"""Demonstrate use of FrechetCellFilter in molecular dynamics exemple."""

import ase
import numpy as np
from ase.md.npt import NPT
from asap3 import Atoms, units
from xtb.ase.calculator import XTB
from ase.io.trajectory import Trajectory
from ase.filters import FrechetCellFilter
from ase.optimize import QuasiNewton

input_file = 'H2_box.xyz'
cell_size = 2

calculator = XTB(method="GFN1-xTB")

atoms = ase.io.read(input_file)
atoms.rattle(stdev=0.01)
atoms.set_cell(np.array([cell_size, cell_size, cell_size]))
atoms.set_pbc(True)
atoms.center(vacuum=1)
atoms.set_calculator(calculator)

ecf = FrechetCellFilter(atoms)

print('len(atoms)) :', len(atoms))
print('len(ecf))   :', len(ecf))

qn = QuasiNewton(ecf)
traj = Trajectory('full.traj', 'w', atoms)
qn.attach(traj, interval=10)
qn.run(fmax=0.5)

print('len(atoms)) :', len(atoms))
print('len(ecf))   :', len(ecf))


print('---- NPT begins ----')
dyn = NPT(
    ecf,
    timestep=0.1 * units.fs,
    temperature_K=300,
    externalstress=0.000101325 * units.GPa,
    trajectory=traj,
    logfile='NPT.log'
    )

dyn.attach(traj.write, interval=10)
dyn.run(50)

Output:

len(atoms)) : 8
len(ecf))   : 11
                Step[ FC]     Time          Energy          fmax
BFGSLineSearch:    0[  0] 09:39:12     -108.889007        6.1586
BFGSLineSearch:    1[  2] 09:39:13     -109.473179        1.1762
BFGSLineSearch:    2[  3] 09:39:13     -109.623021        2.2571
BFGSLineSearch:    3[  5] 09:39:13     -110.469132        2.7407
BFGSLineSearch:    4[  7] 09:39:14     -110.635861        2.8123
BFGSLineSearch:    5[  8] 09:39:14     -110.798278        1.1939
BFGSLineSearch:    6[ 10] 09:39:15     -110.852617        1.2630
BFGSLineSearch:    7[ 12] 09:39:15     -110.907418        1.1981
BFGSLineSearch:    8[ 13] 09:39:16     -110.980643        1.0792
BFGSLineSearch:    9[ 14] 09:39:16     -111.035826        0.1971
len(atoms)) : 8
len(ecf))   : 11
---- NPT begins ----
Traceback (most recent call last):
  File "/Users/marcodigennaro/WORK/TME/PGEL/ASE_tutorials/FrechetCellFilter.py", line 40, in <module>
    dyn = NPT(
          ^^^^
  File "/Users/marcodigennaro/WORK/external_packages/ase/ase/md/npt.py", line 149, in __init__
    MolecularDynamics.__init__(self, atoms, timestep, trajectory,
  File "/Users/marcodigennaro/WORK/external_packages/ase/ase/md/md.py", line 119, in __init__
    self.atoms.set_momenta(np.zeros([len(self.atoms), 3]))
  File "/Users/marcodigennaro/WORK/external_packages/ase/ase/filters.py", line 94, in set_momenta
    mom[self.index] = momenta
    ~~~^^^^^^^^^^^^
ValueError: shape mismatch: value array of shape (11,3) could not be broadcast to indexing result of shape (8,3)

Input file: H2_box.xyz

           8
  Built with Packmol
  H            0.013589       -0.001149        1.023230
  H            0.000832        0.011254        0.323456
  H            0.051896        1.974564        0.714629
  H            0.370690        1.998602        0.091899
  H            1.944549        1.995178        2.000676
  H            1.963944        2.001188        1.300971
  H            1.933873        0.445614        0.043857
  H            1.974646        0.077936        0.638122

Answer 1

I think your NPT object expects a plain Atoms object, and not the Filter object that you get when doing ecf = FrechetCellFilter(atoms).

ecf.atoms will yield the optimized Atoms object without filter, that should work.

Also note that your external stress is 1 eV/AA**3 = 160 GPa.

Edit: Note that none of the NPT classes in ASE implements a "proper" NPT sampling. Just as a word of caution.