[BUG] pair lubricate/poly+brownian/poly leads to strange attractions and a temperature increase

[JitingTian]

Description of the bug

When combining pair lubricate/poly and pair brownian/poly to simulate the dynamics of a binary hard-sphere (HS) mixture, I find that there are strange attractions between close particle pairs and that the system's temperature cannot be kept at the target one. I already realize where is wrong and write new codes to fix the bug, see details below.

LAMMPS Version and Platform

LAMMPS Version: 29Otc2020, the most recent stable version. Platform: Ubuntu 20.10, 8 processors.

Expected and Actual Behavior

(All scripts and results are provided in the attached package) Using Fast Lubricate Dynamics (FLD) to simulate the equilibrium dynamics of a binary hard-sphere mixture, it's expected that the system will go normal suspension Brownian dynamics at a target temperature. For a 50:50 binary HS mixture with identical diameters equilibrated at T=1, the mono (pair lubricate + pair brownian) and poly (pair lubricate/poly and pair brownian/poly) version of FLD should give the same results of a NVE simulation from this configuration.

However, as shown below, poly-FLD will result in a temperature increase in the system.

Observing the atomic trajectories, I find that in the poly-FLD simulation close particles tend to bound together and increase their energies mysteriously, as shown below. I notice this bug is already reported several years ago (https://matsci.org/t/mysterious-attraction-between-hard-spheres-with-pair-lubricate/23325/20), but it still exists in the most recent stable version of LAMMPS.

Fixing this bug

I already fix this bug by writing new codes.

I analyze the source codes of poly-FLD, i.e. pair_lubricate_poly.cpp and pair_brownian_poly.cpp, finding that they do not use newton on to ensure a pairwise lubricate random force. This will result in two different lubricate random forces for a single pair, violating the fluctuation-dissipation theorem, and this is why the temperature cannot be kept. On the other hand, in the source codes of mono-FLD, i.e. pair_lubricate.cpp and pair_brownian.cpp, newton on is used. This difference is shown below. This also means that if we use newton off in the mono-FLD, it will also create unphysical results, as shown below where close and very kinetic particles show up.

I use a two-particle configuration to show this difference. I create 2 atoms with different types but identical diameters, then fix particle #1 and move #2 with a constant velocity. Both particles are under FLD dynamics with flag parameters 0010, so they will experience lubricate force+random force. In principle, the forces of #1 and #2 should be opposite to each other, i.e. F_1=-F_2, satisfying Newton's third law. This is indeed what happens in the mono-FLD simulation, as shown below.

(thermo output of in.pairwise_force_mono) Step Time v_fx1 v_fx2 0 0 622.07516 -622.07516 10 0.001 -230.27202 230.27202 20 0.002 30.310019 -30.310019 30 0.003 -284.452 284.452 40 0.004 111.57893 -111.57893 50 0.005 -0.15459842 0.15459842 60 0.006 104.95937 -104.95937 70 0.007 318.7733 -318.7733 80 0.008 0 0 90 0.009 0 0 100 0.01 0 0 Loop time of 0.000102515 on 1 procs for 100 steps with 2 atoms

However, this does not happen in the poly-FLD simulation, as shown below.

(thermo output of in.pairwise_force_poly_bug) Step Time v_fx1 v_fx2 0 0 622.07516 823.21752 10 0.001 87.610697 -105.55415 20 0.002 190.70195 -475.3001 30 0.003 182.11018 -653.32028 40 0.004 574.92324 245.33342 50 0.005 3.6849458 -147.50467 60 0.006 -90.390638 84.347268 70 0.007 304.09305 -349.23358 80 0.008 0 0 90 0.009 0 0 100 0.01 0 0 Loop time of 9.6339e-05 on 1 procs for 100 steps with 2 atoms

My new codes of poly-FLD use newton on to ensure that F_1=-F_2 holds, creating exactly identical results of mono-FLD, as shown in the attached files. I'm trying to add this modification to LAMMPS, but I'm not familiar with github, so I'm still working on it.

Further Information, Files, and Links

An explanation for the attached package: /codes: new codes for pair lubricate/poly and pair brownian/poly. /verify/many_atoms/: scripts and results of the NVE relaxation of a binary HS mixture, including mono-FLD, old poly-FLD with bug, and my new poly-FLD without bug. /verify/two_atoms/: scripts and results of the 2-particle pairwise force calculation, including mono-FLD, old poly-FLD with bug, and my new poly-FLD without bug. debug_polyFLD.zip