The characteristic length is sqrt(hbar/(m*w)). I think that the previous code had gotten the units mixed up? In any case, the last line in the code below should return True, which it will with this PR. Essentially the force constant matrix is the hessian matrix for the total energy, so if we displace by the characteristic length via E=1/2 <u|FC|u>, the energy should change by the characteristic energy (which is hw/2).
from Inelastica import Phonons
import numpy as np; from pathlib import Path
dyns=np.load("../eph_rse.npy") # or however you get the dynamic atoms
fdfs = [str(p) for p in Path().glob("atm*/RUN.fdf")] # or wherever your fcruns are
D = Phonons.DynamicalMatrix(fdfs ,DynamicAtoms=list(dyns+1))
D.SetMasses()
D.ComputePhononModes(D.mean)
fcmat = D.mean.real[:,:,dyns,:].reshape(72, 72) # I know I have 24*3=72 'springs'
fcmat += fcmat.T
fcmat /= 2
assert np.allclose(D.hw[D.hw > 0], np.diag(D.Ucl.dot(fcmat).dot(D.Ucl.T)).real[D.hw > 0])
The characteristic length is
sqrt(hbar/(m*w))
. I think that the previous code had gotten the units mixed up? In any case, the last line in the code below should return True, which it will with this PR. Essentially the force constant matrix is the hessian matrix for the total energy, so if we displace by the characteristic length viaE=1/2 <u|FC|u>
, the energy should change by the characteristic energy (which ishw/2
).