# Units of get_vibrations().get_energies_and_modes()

Hello,

I’m using ASE for the diagonalization of the hessian matrix. My code is:

import os
from ase.units import Ry
from ase.vibrations import Vibrations
from ase.io.trajectory import Trajectory
from ase.calculators.siesta import Siesta
from ase.calculators.siesta.parameters import Species

from basis import basisSet
from arguments import singlePointArguments

traj = Trajectory('./PBEstructures.traj')

oxygenIndex = [atom.index for atom in traj[0] if atom.symbol == 'O'][0]
hydrogen1Index = [atom.index for atom in traj[0] if atom.symbol == 'H'][0]
hydrogen2Index = [atom.index for atom in traj[0] if atom.symbol == 'H'][1]

for j in range(len(traj)):
system = traj[j]

for k in range(0, oxygenIndex):
if k <= 31:
system[k].tag = 1
elif k >= 32:
system[k].tag = 2

siesta = Siesta(
label=f'structure-{j}',
mesh_cutoff=700*Ry,
kpts=[1, 1, 1],
xc="PBE",
species=[
Species(symbol="H", basis_set=basisSet("H"), pseudopotential="H.psf"),
Species(symbol="O", basis_set=basisSet("O"), pseudopotential="O.psf"),
Species(symbol='Au', basis_set=basisSet('Au2'), pseudopotential='Au.gga.psf', tag=1),
Species(symbol='Au', basis_set=basisSet('Au3'), pseudopotential='Au.gga.psf', tag=1)
],
fdf_arguments=singlePointArguments(),
)

system.calc = siesta

if not os.path.exists(f'structure-{j}'):
os.makedirs(f'structure-{j}')

os.chdir(f'image-{j}')

vib = Vibrations(
atoms=system,
indices=[ oxygenIndex, hydrogen1Index, hydrogen2Index ],
name=f'structure-{j}',
delta=0.02,
nfree=2
)

vib.run()
print(vib.get_vibrations().get_energies_and_modes())
vib.summary()

os.chdir('../')

system = None
siesta = None


PBEStructures.traj is a file containing a box with 64 gold atoms and a water a molecular on top of it.

My goal here is to understand the units of measure (mainly of the eigenvectors) used when I print vib.get_vibrations().get_energies_and_modes(). The output I receive is:

(array([0.        +0.02172963j, 0.        +0.00898803j,
0.0047047 +0.j        , 0.00515926+0.j        ,
0.01071083+0.j        , 0.01819307+0.j        ,
0.20329078+0.j        , 0.46236624+0.j        ,
0.4736801 +0.j        ]), array([[[-4.42204141e-02,  8.21599997e-04, -8.20717964e-04],
[ 1.88182925e-01,  4.34003246e-02,  6.77329557e-01],
[ 1.87967471e-01, -2.82722309e-02, -6.54830052e-01]],

[[-1.30904884e-03, -1.17128050e-01,  8.54004482e-02],
[-1.48444828e-02,  1.62819355e-01,  5.50239486e-01],
[-2.37987090e-02,  1.16354127e-01,  5.61187346e-01]],

[[ 2.30667250e-01,  1.59479522e-02,  1.83724813e-02],
[ 2.57152279e-01,  2.91414458e-03,  6.41077081e-02],
[ 2.57173126e-01,  2.47372183e-02, -3.41282499e-02]],

[[-2.12276397e-02,  2.00084623e-01,  1.18940029e-01],
[-3.12744446e-02,  2.17942020e-01,  1.27328885e-01],
[-3.23472198e-02,  1.96342333e-01,  1.43769534e-01]],

[[ 9.45816964e-03,  5.59756452e-02, -1.97466996e-01],
[-8.46615478e-03,  3.37616114e-01,  2.52215042e-01],
[-1.63751228e-02,  2.77428558e-01,  2.59051203e-01]],

[[ 4.77412315e-02, -1.63801442e-03,  4.83822863e-03],
[-3.29851589e-01,  5.74999313e-01,  4.96722335e-02],
[-3.28998731e-01, -6.24454755e-01, -1.24080348e-01]],

[[ 1.21337881e-05,  5.97117696e-02, -3.33461626e-02],
[ 4.09683326e-01, -4.74611830e-01,  2.56556369e-01],
[-4.09675251e-01, -4.74696304e-01,  2.56387055e-01]],

[[-2.03379881e-04, -4.23625906e-02,  2.44361192e-02],
[ 5.74482847e-01,  3.37161265e-01, -1.90362467e-01],
[-5.71214032e-01,  3.34886317e-01, -1.89000171e-01]],

[[-6.72311778e-02,  1.32489072e-04, -6.62638940e-05],
[ 5.31582749e-01,  3.64169723e-01, -2.05093709e-01],
[ 5.35119980e-01, -3.66289334e-01,  2.06244020e-01]]]))
---------------------
#    meV     cm^-1
---------------------
0   21.7i    175.3i
1    9.0i     72.5i
2    4.7      37.9
3    5.2      41.6
4   10.7      86.4
5   18.2     146.7
6  203.3    1639.7
7  462.4    3729.2
8  473.7    3820.5
---------------------
Zero-point energy: 0.589 eV


Thank you in advance,

Regards,

João Rheinheimer

Hi João,

the energies should be in eV as this is the “standard” ASE unit for energy; that looks consistent with the summary table in meV.

For the eigenvectors let’s look at the implementation at ase/vibrations/data.py · 2838969fc26938f94813a454ffdbd910604ae458 · ase / ase · GitLab - it is only a few lines long!
The raw eigenvectors from np.linalg.eigh are already normalised, then they are multiplied by the square roots of the atomic masses to obtain the “modes”.

This makes things simple for iter_animated_mode(): the eigenvectors are scaled by \sqrt{k_B T / \epsilon} to obtain physical-looking displacements. I’m not sure this is enough to get a “correct” distance scaling though? For visualisation purposes a bit of exaggeration is generally useful…

If you want to get back to the normalised eigenvectors you can multiply by the square roots of the masses again.