Short description

Calculates the power spectra from the velocity autocorrelation function. In addition, it uses forceconstants to fourier-interpolate the q-resolved spectra non the non-exact q-points. What that means is that it can generate continuous spectralfunctions as a function of wavevector not limited to the commensurate modes.

Command line options:

Optional switches:

• --noforceconstant
  default value .false.
  Do not use forceconstants for Fourier interpolation.

• --dos
  default value .false.
  Outputs the DOS. This requires a specification of a q-point mesh.

• --unit value, value in: thz,mev,icm
  default value thz
  Choose the output unit. The options are terahertz (in frequency, not angular frequency), inverse cm or meV.

• --loto
  default value .false.
  Add the corrections due to long-range electrostatic interactions, requires infile.lotosplitting Format is specified here.

• --nq_on_path value, -nq value
  default value 100
  Number of q-points between each high symmetry point

• --readpath, -rp
  default value .false.
  Read the q-point path from infile.qpoints_dispersion. Use crystal structure into to generate an example.

• --qpoint_grid value#1 value#2 value#3, -qg value#1 value#2 value#3
  default value 26 26 26
  Density of q-point mesh for DOS interpolation.

• --maxfrequency value
  default value 1.3
  Changes the maximum frequency. Default is 1.3 times the maximum frequency determined from the supplied forceconstants.

• --help, -h
  Print this help message

• --version, -v
  Print version

Examples

autocorrelation

autocorrelation

Long description

The (discrete) velocity autocorrelation function is defined for atom $i$ and Cartesian component $\alpha$ as

$$ C^{i\alpha}_v(\tau) = \sum_{t=1}^{N_t} v^{\alpha}_{i}(t)v^{\alpha}_{i}(t+\tau) $$

Where $\tau$ is the lag time and $v_i(t)$ the velocity of atom $i$ at time $t $. We want this resolved in reciprocal space in the frequency domain:

$$ C(\mathbf{q},\omega) = \sum_\alpha \sum_{i} \mathcal{F}\left\{\ C^{i\alpha}_v \right\}^2 e^{i\mathbf{q}\cdot \mathbf{R}_i} $$

where $\mathcal{F}$ is the discrete Fourier transform over lag time. Here we exploited the fact that the Fourier transform of a convolution is a product in the frequency domain.

Without loss of information, the velocities can be transformed to normal mode coordinates:

$$ v_{\mathbf{q}s}^{\alpha} = \sum_i \epsilon_{\mathbf{q}s}^{i\alpha} v_{i}^{\alpha} $$

where the modes $\mathbf{q}s$ are confined to the exact q-points and modes of the simulation cell. This gives us the mode-resolved velocity autocorrelation function

$$ C^{\mathbf{q}s}_v(\tau) = \sum_{\alpha} \sum_{t=1}^{N_t} v_{\mathbf{q}s}^{\alpha}(t) v_{\mathbf{q}s}^{\alpha}(t+\tau) $$

and the mode resolved power spectra

$$ C(\mathbf{q},s,\omega) = \mathcal{F}\left\{\ C^{\mathbf{q}s}_v(\tau) \right\}^2 $$

which is easily interpolated to an arbitrary q-vector. The code will output the raw power spectra as well as the interpolated.

Practical considerations

Please note that autocorrelation function spectra require careful calculations. It is far better suited for classical molecular dynamics rather than Born-Oppenheimer MD. The frequency resolution scales as the simulation time, requiring very long simulation times. The (typical) small sizes of first principles MD introduce odd artifacts. Some simple convergence tests shows that the phonon lineshapes from perturbation theory and correlation functions agree at about 10000 atoms and 1ns simulation times.

Do not use this code to try to determine phonon lineshapes/lifetimes unless you really know what you are doing.

Input files

These files are necesarry

• infile.ucposcar
• infile.ssposcar
• infile.meta
• infile.stat
• infile.positions
• infile.forces
• infile.forceconstant

Optional file to specify a path for the phonon dispersions

• infile.qpoints_dispersion

For including the electrostatic corrections

• infile.lotosplitting

Output files

outfile.acf_sqe.hdf5

The file is more or less self-explainatory, the following is a matlab snippet that produces a neat plot:

% read everything from file
fn=('outfile.acf_sqe.hdf5');
x=h5read(fn,'/q_values');
y=h5read(fn,'/energy_values');
gz=h5read(fn,'/intensity');
xtck=h5read(fn,'/q_ticks');
xtcklabel=strsplit(h5readatt(fn,'/','q_tick_labels'));
energyunit=h5readatt(fn,'/','energy_unit');

% plot the results

figure(1); clf; hold on; box on;

[gy,gx]=meshgrid(y,x);
s=pcolor(gx,gy,log10(gz+1E-2));
set(s,'edgecolor','none','facecolor','interp')
set(gca,'xtick',xtck,'xticklabel',xtcklabel)
ylabel(['Energy (' energyunit ')'])
xlim([0 max(x)])
ylim([0 max(y)])

outfile.acf_highsymmetrypoint_??

In case the exact q-points of the simulation cell coincide with high-symmetry points, the spectralfunction at those points will be printed:

That is, the first column is the energy, the following are intensities per mode.

outfile.acf_raw_dos

With option --dos this will be written, containing the raw power spectrum without any interpolation procedure. It usually looks awful.

The first column is energy, the second the total DOS. In case there is more than one atom per unit cell, the site projected DOS is printed as well.

outfile.acf_interpolated_dos

The format is identical to the raw dos above.