VDOS results do not follow logical trends

Hi there,

I am trying to get phonon vibrational density of states (VDOS) from LAMMPS dump file. I have tried several approaches, either using LAMMPS commands to get the VACF and convert it into the VDOS with python codes, or using the available codes in web for this purpose. From all of these methods, I am getting a VDOS distribution which does not show zero intensity at zero frequencies.

For getting the trajectory files, I am using NVE with 0.1 time steps and for 200000 run times.

#  This is the control script for LAMMPS.

echo			both

#-------------------------------------------------------------------------------
# Stage 2.1: Initialize LAMMPS run for 3-d periodic
#-------------------------------------------------------------------------------

units			real
boundary		p p p
atom_style		full
#restart 200        restart

pair_style		lj/class2/coul/long 10 15
pair_modify             mix sixthpower
pair_modify		tail yes
bond_style		class2
angle_style		class2
dihedral_style		class2
improper_style		class2
special_bonds           lj 0.0 0.0 1.0 coul 0.0 0.0 1.0

box                     tilt large
read_data		structure.dat
#read_restart	 	structure.dat

include			parameters.dat

neighbor		2.0 bin
neigh_modify		delay 0 every 1 check yes
kspace_style		pppm 0.0001

#-------------------------------------------------------------------------------
# Stage 2.2: Minimization
#-------------------------------------------------------------------------------

#min_style		cg
#min_modify		dmax 0.05 line quadratic
#reset_timestep		0

#minimize		0.000001 0.000001 10000 10000
#velocity        all create 300.0 497 dist gaussian

reset_timestep 0

velocity all scale 300.0

timestep           0.01

fix 2 all nve

compute vacf all vacf
fix 3 all ave/time 10 10 100 c_vacf[*] file vacf.dat

thermo	           10
thermo_style       custom step cpu temp pe ke epair ecoul  etotal density lx ly lz
dump		   1 all custom 5 dump.xyz id type x y z vx vy vz
dump_modify	   1 sort id

run 		   20000


Could anybody please assist in what could get wrong for getting this result?

I am very sorry if this question is not quite relevant to this forum. But I will really appreciate any thoughts on this.

Thank you very much for your assistance, in advance.

Note that

  1. The velocity autocorrelation function is usually dominated by an exponential decay signal (in the time domain). This makes sense – particles’ velocities do not instantaneously change randomly; seeing a velocity at time t gives you information about its velocity a short time later, with that information’s validity decaying over time.
  2. The Fourier transform of a signal with non-zero mean will have a non-zero peak at frequency zero.

I do not remember enough physics offhand to tell you how to interpret this peak, but it is mathematically expected and correct.

Dear @srtee thank you very much for your response.

Does this mean that I must equilibrate the system for longer times before running NVE for VACF analysis? Or this trend is normal? Because I have checked many papers in the literature and almost all plots have a zero value at zero frequency.

Thank you very much for your valuable comments,

Bests,

Mahsheed

I’m worried that your post shows a foundational misunderstanding of the velocity autocorrelation function.

Let’s say that in your simulation a particle is moving at time 0 with the velocity 0.1 angstrom/femtosecond in the x-direction.

One timestep later, the velocity will be almost exactly the same. Let’s say it’s (0.0999, 0.002, -0.004) instead of (0.1000, 0.0, 0.0).

A hundred timesteps later, the velocity should still be quite similar but have changed somewhat – let’s say (0.070, 0.010, -0.040).

By ten thousand timesteps the velocity would be quite different.

This is what I mean by “an exponential velocity autocorrelation function” – it will be 1 at time 0, 0.9 at a hundred timesteps, … The area under that curve will generally not be zero. Therefore, the Fourier transform of the VACF will generally not be zero at zero velocity.

I hope you can see that this has nothing to do with equilibration, NVE vs NVT integrators … this is a fundamental property of microscopic thermal motion.

This means that the papers you cite perform some additional transformation after the Fourier transform. From my offhand memory the FT is multiplied by frequency squared to get the VDOS – that would give zero at zero frequency. For best results, look at these papers you are looking at and search for a formula that shows you exactly what calculations to do.

1 Like

Dear @srtee,

Thank you very much for your informative response. I took your advice and checked some procedures to remove the peak at the zero frequency. However, my other concern is about the frequency value itself. From my analysis, the frequency at which peaks happen, is in the high frequency ranges (in the order of 100 to 200 THz). However, I should be able to see some vibrations in the smaller frequencies, corresponding to the acoustic movements. I was wondering if this could be related to the time steps of my simulation or not? I have actually tried different time steps, but they seem not to be capturing the small frequency vibrations. Is there any considerations that I need to take in mind?

I sincerely appreciate your guidance on this matter.
Bests,