Phonon density of states from velocity autocorrelation function

Lammps Users,
I am attempting to calculate the phonon density of states from the
velocity autocorrelation function, but I obtain some strange results.
After running the lammps script below, I obtain an fft of the VACF output.
The phonon section for Ni is reproduced, however, I find is a very long
1/f-type tail at high frequencies.

Does anyone have any experience with this technique? Has anyone run into
the long-tail phenomena or have a suggestion on how to avoid it? I have
tried the obvious: increased system size, run for longer times, used
smaller time steps, averaging multiple VACFs, and employed zero padding.
Nothing I have tried removes or changes the tail. I have also tried
multiple Ni potentials, but the tail remains.

Thanks in advance,
Chris O'Brien

### Phi=0.000000
variable t equal 400.0
variable alat equal
3.52296665615187+2.98018854145148e-5*$t+1.82740327672739e-08*t^2 units metal boundary p p p atom\_style atomic \#\#\# set up simulation domain lattice fcc {alat}
region box block 0.0 10.0 0.0 10.0 0.0 10.0 units lattice
create_box 1 box
create_atoms 1 region box
# EAM potential -- Ni
pair_style eam/alloy
pair_coeff * * /usr/local/lammps-pot/ni1.set Ni
mass * 58.700
#re-neighbor
neighbor 1.5 bin
neigh_modify delay 5
#Computes
compute PE all pe/atom
compute xPE all pe/atom
compute sumpe all reduce sum c_xPE
compute xke all ke
variable xetotal equal c_xke+c_sumpe
compute gbtemp all temp
compute CNA all cna/atom 3.033
compute CENTRO all centro/atom fcc
dump 3 all custom 5000 dump.phi_0.0 id type x y z c_CNA
c_CENTRO c_PE
dump_modify 3 sort id
thermo 1000
#thermo_style custom step pxx pyy pzz lx lz c_gbtemp v_xetotal
velocity all create $t 19277
#fix 7 all nve
#fix 8 all langevin $t $t 0.1 19277
#fix 8 all npt temp $t $t 0.01 iso 1.0 1.0 0.1
fix 8 all nvt temp $t $t 0.01
timestep 0.0001
# equilibrate
run 100000
unfix 8
fix 7 all nve

variable i loop 10
        label ensembleloop
        #store initial velocity for each instance of vacf calculation
        fix v0save all store/state 10000 vx vy vz #save
velocitities every 10000 steps
        variable v0v0 atom
f_v0save[1]*f_v0save[1]+f_v0save[2]*f_v0save[2]+f_v0save[3]*f_v0save[3]
        compute v0v0ave all reduce ave v_v0v0
        compute velcorr all vacf
        variable normvelcorr equal c_velcorr[4]/c_v0v0ave
        thermo 100
        thermo_style custom step pxx pyy pzz lx lz c_gbtemp v_xetotal
v_normvelcorr c_v0v0ave
        fix corrout all ave/time 1 1 1 v_normvelcorr file
normvelcorr.$i
        run 10000
        uncompute velcorr
        uncompute v0v0ave
        next i
        jump SELF ensembleloop

Lammps Users,
I am attempting to calculate the phonon density of states from the
velocity autocorrelation function, but I obtain some strange results.
After running the lammps script below, I obtain an fft of the VACF output.
The phonon section for Ni is reproduced, however, I find is a very long
1/f-type tail at high frequencies.

how did you do the FFT?

did you "symmetrize" the VACF data before the transform?

Lammps Users,
I am attempting to calculate the phonon density of states from the
velocity autocorrelation function, but I obtain some strange results.
After running the lammps script below, I obtain an fft of the VACF
output.
The phonon section for Ni is reproduced, however, I find is a very long
1/f-type tail at high frequencies.

how did you do the FFT?

did you "symmetrize" the VACF data before the transform?

I performed the FFT using the numpy¹s fft function ³rfft² to perform the
FFT. This function is meant for real data only.

Unfortunately, I forgot to symmetrize the VACF. The lack of symmetry seems
to have caused the long 1/f-type tail.

You could also compare your results to what the USER-PHONON
package and its fix phonon command can do.

Steve