Dear Aidan,
Thank you for your email and suggestions.
How about for using the volume fluctuations (dv) in NPT ensemble to compute the bulk modulus and then the sound velocity as:
E=kB * T * V / <dv^2> and then c=(E/rho)^1/2.
The problem in LAMMPS is that strongly depends on choosing Pdamp factor (e.g. 0.1 or 0.01 in units metal) which results in different values of c.

Dear Aidan,
Thank you for your email and suggestions.
How about for using the volume fluctuations (dv) in NPT ensemble to compute
the bulk modulus and then the sound velocity as:
E=kB * T * V / <dv^2> and then c=(E/rho)^1/2.
The problem in LAMMPS is that <dv> strongly depends on choosing Pdamp factor
(e.g. 0.1 or 0.01 in units metal) which results in different values of c.

wouldn't it be technically easier to compute the bulk modulus from
finite differences?

i.e. equilibrate your system with npt and take a few (decorrelated!)
snapshot restarts from when the system is in equilibrium.

then run with fix nvt equilibrate a bit more and record pressure until
sufficiently converged
change box from restart to a slightly smaller volume, equilibrate a
bit and record pressure until sufficiently converged again.
change box from restart to a slightly larger volume, equilibrate a bit
and record pressure until sufficiently converged again.

compute: -V*(delta p)/(delta V) for all three points, fit to a
parabola and get -V*dp/dV from it.
repeat with additional snapshots until converged.

with the partition option in lammps this can be run in a very parallel
fashion and thus very effectively parallelized.

The idea of using volume fluctuations has been around since the early
days of MD. The only problem is, it does not work. Many papers have
been published confirming this fact, looking at lots of different
variations on the Nose-Hoover barostat, and in all cases, the
simulations do not reliably converge to a single value. It seems to
me that there is a problem with ergodicity.

Axel's method will work, if care is taken. The hard part is balancing
the need for a statistically accurate <deltaP>, with the requirement
that deltaP be small. At zero temperature, all measurements are
exact, so statistical efficiency is not an issue.

Getting elastic constants from pressure fluctuations in an NVT
simulation requires computing an additional term that is related to
the second derivative of energy i.e. d2E/dRi2. This is not easy to do,
because LAMMPS only computes F = -dE/dRi. But I believe it converges
well. See the paper by Yubao Zhen:

"A deformation-fluctuation hybrid method for fast evaluation of
elastic constants with many-body potentials"
Yubao Zhen *, Chengbiao Chu
Computer Physics Communications 183 (2012) 261-265

The second derivative terms will also add the extra complication that they would have to be uniquely defined. Potentials like reaxff I believe do not have continous second derivatives. If features like this were to be implemented this would probably need to show in a section of the manual where the limitations of the method will appear at the top of the page or in BIG RED letters at the bottom.
Carlos

Getting elastic constants from pressure fluctuations in an NVT
simulation requires computing an additional term that is related to
the second derivative of energy i.e. d2E/dRi2. This is not easy to do,
because LAMMPS only computes F = -dE/dRi. But I believe it converges
well. See the paper by Yubao Zhen:

“A deformation-fluctuation hybrid method for fast evaluation of
elastic constants with many-body potentials”
Yubao Zhen *, Chengbiao Chu
Computer Physics Communications 183 (2012) 261-265

I don't see any reason to be worried about discontinuities in
d2E/dRi2. If there are jumps, these will just be averaged over during
the sampling run. And why are you singling out ReaxFF for the
red-letter treatment? Even lj/cut has discontinuous *first
derivatives*.

Reaxff was just the example that came to mind. Certainly was not intending to single out reaxff at all. I guess maybe my answer should have been placed in a new thread as what I intended to point to was not really related to the current problem of the sound velocity which is totally out of my interest scope. My issue with discontinuous second derivatives applies more to features that are not yet implemented in Lammps. One of them is the computation of Hessian eigenvectors with small eigenvalues without explicitly computing the Hessian as it is done in hyperdynamics. Reaxff is a potential where I have noticed that using an iterative procedure to compute such Hessian eigenvectors fails to converge in many occasions. I am going to leave it here as I have no intention of making this a topic of discussion in the forum.

Returning Back to the calculation of the elastic constants I just wanna say that It would be really great (WISH LIST) if a FIX could be implemented. Something like “FIX econst” . i.e Calculations of the elastic constants in the background and the information be available into the LOG FIle or the Thermo Output . You could use the work by John Ray : http://prb.aps.org/abstract/PRB/v32/i2/p733_1 . Although its more easy to make your own post processing script …

This would indeed be nice. The problem is, Ray's method uses
"fluctuation formula involving the internal stress tensor," which
requires d2E/dR2, which would have to be coded up separately for each
potential. The paper by Yubao Chen provides a solution to this, but we
never got it working properly in LAMMPS.