Hello everyone,
I’m encountering some difficulties in computing the velocity autocorrelation function (VACF) using LAMMPS.
As mentioned in this discussion: Which one is better, VACF or ave/correlate of velocity?,
"fix ave/correlate" is a newer, more general method in LAMMPS that allows an arbitrary number of prior reference values to be compared to the current value, whereas "compute vacf" only uses the initial value as the reference.
However, I am unsure how fix ave/correlate actually computes the VACF. To correctly average over different starting times, it should also compute the product of the velocity at the current time and the velocity at the starting time, then average over all atoms.
Another discussion, fix ave/correlate for VACF, suggests using compute vacf. However, as noted in the previous discussion, compute vacf does not perform averaging over different starting times, leading to significant noise in the results. This can cause issues when computing the phonon density of states (PDOS). Furthermore, using fix 1 all ave/time 1 10 100 c_vacf[*] file vacf.dat does not achieve averaging over different starting times either.
One common workaround is to output the velocity of all atoms at all time steps, as done in this open-source project: VACF_PDOS. However, this approach becomes impractical for large systems, as it results in massive output files, making storage, downloading, and processing very time-consuming.
Is there a more efficient way to compute VACF within LAMMPS that I might have overlooked? Or am I misunderstanding how fix ave/correlate works in this context? Below is an example script I used for computing VACF. The file cnt.data is attached
file.lmp (1.2 KB)
cnt.data (1.4 MB)
# LAMMPS (29 Aug 2024 - Update 1)
log log.lmp
units metal
dimension 3
atom_style atomic
boundary m m m
read_data cnt.data
pair_style airebo 3 0 0
pair_coeff * * /usr/share/lammps/potentials/CH.airebo C
region rEdgeBottom block INF INF INF INF INF 10
region rEdgeTop block INF INF INF INF 190 INF
region rPDOS block INF INF INF INF 80 120
region rBottom block INF INF INF INF 10 30
region rTop block INF INF INF INF 170 190
group gEdgeBottom region rEdgeBottom
group gEdgeTop region rEdgeTop
group gEdge union gEdgeBottom gEdgeTop
group gPDOS region rPDOS
group gBottom region rBottom
group gTop region rTop
group gMove subtract all gEdge
group gCenter subtract all gEdge gBottom gTop
velocity all create 300 1 rot yes mom yes dist gaussian
fix 1 all nvt temp 300 300 $(100.0*dt)
thermo 1000
run 10000
unfix 1
fix 1 gMove nve
fix 2 gEdge setforce 0 0 0
velocity gEdge set 0 0 0 sum no
thermo 1000
run 10000
# way 1
compute 2 gPDOS vacf
fix fv gPDOS ave/time 1 10 10 c_2[*] file vacf.dat
# way 2
# fix 1 gPDOS ave/correlate 1 10 10000 ? ? ?
# way 3
# https://github.com/lingtikong/VACF_PDOS/
dump 1 gPDOS custom 5 dump.vc vx vy vz
dump_modify 1 sort id
thermo 1000
run 10000
Thank you in advance for your help!
Hang Yu.