trap() function

Dear Lammps users,

I have been performing Green-Kubo calculations with Lammps to calculate the thermal conductivity but have been getting different values than that of non-equilibrium calculations for certain systems.

I have been taking the auto-correlation function which is dumped by lammps during the ave/correlate and integrating it via trapezoidal integration and applying the relevant scaling factors via an extremely simple external Fortran90 program.

For a particular example of SPC/E water at ambient conditions the values which I am getting appear to be different from those that Lammps gets after integration by a factor of 1.15. I have read the relevant Lammps documentation to fully interpret exactly how Lammps executes the input lines thus evaluating the integral and am completely clueless as to why my calculation deviates from that of Lammps.

I have attached the Lammps input script "in.spce", a particular iteration of the heat flux auto-correlation function "jj300.dat" and the fortran90 program "intjj.f90" that I am using to integrate it in the hope that someone can deduce if and what I am doing wrong.

The program by nature is extremely short and I have fully annotated it to minimise the amount of time it should take to interpret.

Many thanks in advance, Jeff.

in.spce (1.91 KB)

jj300.dat (21.9 KB)

intjj.f90 (1.19 KB)

I have been performing various simple tests to determine the origin of the difference between my trapezoid integration calculation and that of Lammps for the heat flux autocorrelation function.

I have come to the conclusion that it is an effect of the difference of significant figures which Lammps and myself have for the function.

I noticed that the output for the thermal conductivity calculated internally by Lammps is given to 8 significant figures, whereas the output of the heat-flux autocorrelation function in the standard example of Green-Kubo thermal conductivity documented on the Lammps website gives the discrete values to 6 significant figures.

The differences of performing trapezoid integration with these differences in accuracy (especially for phase-points that contain negative correlation and thus summations whose result is orders of magnitude smaller than that of the average modulus of the function) appears to be very significant (I have seen deviations of ~18% for SPC/E water).

Obviously when dealing with Green-Kubo calculations there is no clear a priori way of determining a "good" upper integration limit so ideally one would want to perform the integration after inspection of the apparent region of convergence.

With this in mind, I was hoping to find out if Lammps has some sort of "flag" to increase the precision of the "ave/correlate" fix's output (similar to the case of the dump command), but could not see anything in the documentation for this particular fix.

Many thanks, Jeff.

Internally, fix ave/correlate is doing everything in full double precision.
There is no flag for the output precision to a file, though it wouldn't
be hard to add. However, you could easily look at individual
values by printing them with thermo output, where you can control
the precision (thermo_modify format). That would test your hypothesis.

Steve