4
$\begingroup$

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
$\endgroup$

1 Answer 1

5
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thank you @flonko. This solves the mismatch issue, but the volume and stress are not changing as I'd expect. $\endgroup$ Oct 24, 2023 at 12:32
  • $\begingroup$ In the new version of the question, I show how neither NVT nor NPT change volume or pressure. Do you know why? $\endgroup$ Oct 24, 2023 at 12:47
  • $\begingroup$ @MarcoDiGennaro please ask a new question. I rolled back your edit. $\endgroup$ Oct 24, 2023 at 13:48
  • $\begingroup$ @Nike Dattani, fine by me. I think we could also modify the title otherwise it's misleading. How about "shape mismatch in ASE filter"? $\endgroup$ Oct 24, 2023 at 13:57
  • $\begingroup$ @MarcoDiGennaro great idea! You can go ahead with that title change! $\endgroup$ Oct 24, 2023 at 15:19

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .