Hello everyone,
I've made a crude patch to the function add_nonbonded_pair_virials() to calculate Helfand moments (sometimes called Einstein-Kubo-Helfand moments).
This is basically a slight adition to the stress tensor which allows direct calculation of shear viscosities in the form:
\eta \propto <G^2>
by analogy with the usual method for calculating diffusion constants D \propto <x^2>. I've been working from Erpenbeck PRE 1995 for the formula, and validating my answers against Meier et al. JCP 2004.
Espresso already provides a means to find the shear viscosity, of course, using a Green-Kubo method starting from the stress-tensor autocorrelation function, but this is quite a messy and difficult calculation, given the tendency of the acf to quite often have a long and noisy tail.
EKH doesn't seem to be any cheaper than GK, but it is simpler for the user: there is less onus to e.g. fit something to the acf and integrate the fitted function. (Or for that matter to notice in the Espresso docs that online acfs are available).
The disadvantage is that it needs information from the verlet nonbonded force calculation, including the *un-imaged* particle coordinates, as well as the force and the imaged pairwise distance. I've been using the function
get_particle_data(p1->p.identity, &pp1);
to retrieve the un-imaged coordinates, however this of course introduces a performance and scaling penalty, meaning that the functionality will have to be hash-deffed out by default if it is integrated with the main build.
Is there any interest in having this? Shall I create a branch and add it in?
Josh