Velocities has a wrong shape

[fekad]

It looks like the velocities ('velo' in ase's arrays dictionary ) stored in a transposed way.

from ase import Atoms, units
from ase.io import write

from ase.lattice.cubic import Diamond

from ase.md.langevin import Langevin
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

from quippy.potential import Potential
# Build Si lattice
size=(1,2,2)
lattice = Diamond('Si', latticeconstant=5.44, size=size)
atoms = Atoms(lattice)

# attach tight binding calculator
qm_pot = Potential('TB DFTB', param_filename='/opt/quip/doc/Tutorials/params.xml')
atoms.set_calculator(qm_pot)

T = 1000.0 # Temperature [Kelvin]
timestep = 1.0 * units.fs

MaxwellBoltzmannDistribution(atoms, 2.0 * T * units.kB)
dynamics = Langevin(atoms, timestep, T * units.kB, 0.002)
dynamics.run(steps=1)

print(atoms.arrays['velo'].shape, atoms.arrays['positions'].shape)

(3, 32) (32, 3)

write('test.xyz', atoms)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-aa5bc91af6c4> in <module>
----> 1 write('test.xyz', atoms)

/opt/conda/lib/python3.7/site-packages/ase/io/formats.py in write(filename, images, format, parallel, append, **kwargs)
    403 
    404     _write(filename, fd, format, io, images, parallel=parallel, append=append,
--> 405            **kwargs)
    406 
    407 

/opt/conda/lib/python3.7/site-packages/ase/io/formats.py in _write(filename, fd, format, io, images, parallel, append, **kwargs)
    436                 mode = mode.replace('w', 'a')
    437             fd = open_with_compression(filename, mode)
--> 438         io.write(fd, images, **kwargs)
    439         if open_new:
    440             fd.close()

/opt/conda/lib/python3.7/site-packages/ase/io/extxyz.py in write_xyz(fileobj, images, comment, columns, write_info, write_results, plain, vec_cell, append, tolerant)
    862             else:
    863                 for c in range(ncol):
--> 864                     data[column + str(c)] = value[:, c]
    865 
    866         nat = natoms

ValueError: could not broadcast input array from shape (3) into shape (32)

[jameskermode]

This is related to #144, so a similar fix of listing velo as a property not to copy back to ASE would work. But perhaps we need a more general fix which tranposes any 3xN arrays. I’ll take a look.

[stenczelt]

I have taken a look now and can fix it easily.

If velocities were copied, they need to be converted to ase units or actually to momenta. There are back and forth converters implemented for this in quippy.convert but I think the potential object should not be doing this, but the dynamics, which is indeed dong it.

[jameskermode]

Thanks @stenczelt - can you try for a general fix, since the same issue cam up with force in #144

[stenczelt]

arrays not in _skip_keys are now transposed and velocities are not updated by the potential object, should use atoms.get_velocities() to obtain velocities, which is in ase units, can be converted to QUIP units with quippy.convert.velocities_ase_to_quip()

[jameskermode]

Thanks Tamas, looks good!

[fekad]

Thank you both for the quick fix! It works like a charm now.