I recently conducted a comparative analysis of energy calculations in LAMMPS simulations, specifically focusing on the utilization of the mix arithmetic command. In one case, I provided all mixed pair coefficients directly in the input file, while in the other, I employed mix arithmetic and only included diagonal pair coefficients (i.e., 1 1, 2 2, and so forth) in the input file. Both simulations utilized identical initial configurations and seed values for velocity generation.
During my investigation, I observed the following discrepancies:
Energy Deviations: When comparing the energies calculated in both cases, I noted a marginal difference on the order of 10^(-7) kcal/mol. I believe this variance is negligibly small and may be ignored.
Disparity in (U1-U0) Calculation: To compute the difference in potential energy (U1-U0) for a delta of 0.002, I employed the compute fep command. Surprisingly, I encountered a noticeable difference in the order of 0.1 kcal/mol between the two cases. Such a significant variance raises concerns regarding the reliability of our energy calculations, particularly when employing mix arithmetic.
To facilitate a comprehensive analysis, I have attached all relevant files. Any insights or feedback on this matter would be greatly appreciated.
NB: I’m using LAMMPS version 17 Apr 2024. I didn’t able to attach the energy.dat file for both case due to the number of file restrictions.
on Point 1:
Some base 10 floating-point numbers cannot be exactly represented in the base 2 floating point format that CPUs use internally (e.g. 0.1), so when entering mixed numbers for potential parameters, there is a good chance that there will be tiny differences, since the mixing is done with full double precision. The magnitude of energy differences you see is within that range.
on Point 2:
Please keep in mind that MD is doing numerical integration of coupled partial linear differential equations. These are essentially chaotic systems and thus the tiniest difference will grow exponentially and thus trajectories will diverge. This can happen due to the difference between mixing and explicitly providing the mixed terms, different random seeds, different number of processors. Different platforms/compilers/optimization levels. When you do “statistical” calculations like FEP, the error of your results is based on the convergence of your statistics. To gauge how well converged you are, it is common to multiple simulations with different initial conditions. The difference you see is thus not so much caused by the choice of mixing, but general convergence of the FEP methods. It should get smaller if you increase the statistical sample. You would likely also get difference with different starting settings. Bottom line, it is not a function of how the potential parameters are set, but of how well your results are converged.
As context, “chemical accuracy” is usually defined as 1-2 kcal/mol (John Pople’s Nobel Lecture, section VI, puts it at 0.02 kcal/mol per valence electron for DFT), and is still generally considered difficult (but within reach!) for even quantum chemistry methods.
It’s absolutely worthwhile to quantify the numerical accuracy of your simulations, but a difference of 0.1 kcal/mol between replicas is not likely to dominate either the random or systematic errors in your calculations.