[ESPResSo] LJ potential shift in the sample scripts

[Peter Kosovan]

Dear all

I have noticed a disagreement between the calculation of the LJ shift in 
the sample scripts (both in versions 2.0.2n and 1.9.7h) and the 
documentation/implementation. In the sample pe_solution.tcl, the shift is 
calculated as

set lj2_shift [expr -4.0*$lj2_epsilon*(pow($lj2_cut,-12)-pow($lj2_cut,-6))];

but the documentation says

U_{LJ}(r)=4\epsilon ((sigma/r)^{12}-(sigma/r)^{6}+shift)

the implementation on line 237 of lj.h is

return 4.0*ia_params->LJ_eps*(SQR(frac6)-frac6+ia_params->LJ_shift);

which is the documented version but disagrees with how the shift is 
calculated in the sample scripts.

When the shift from the sample scripts is used with a cutoff 2.5*sigma, it 
produces a step in the energy of 0.05*epsilon (3/4*epsilon*lj2_shift). If 
one uses the Langevin thermostat, the discontinuity in the energy is 
presumably lost in the noise produced by the theromostat, so the error is 
difficult to notice.

With regards

Peter Kosovan