General questions about espresso and its implementation of LBM

[rdapena]

Hello!

I've recently started studying the LBM and I was hoping to use ESPRESSO to run some simulations. However, after reading through the documentation, I have some questions:

• Regarding the units to be used in a simulation, it's my understanding that it is up to the user to choose a proper scaling for the basic magnitudes (length, time, etc.). So, if I were to set $C_L =10^{-3} \ m$ and then define an arbitrary variable $x=2$ that would imply that x's physical value is $x=2\ C_L$, right?
• Also related to the topic of units, in one of your samples, you define a lb fluid as follows:

  system.time_step = 0.01
  lbf = espressomd.lb.LBFluidWalberla(agrid=1.0, density=1.0, kinematic_viscosity=1.0, tau=0.01)
  

  However, it is stated in LBM's literature that $\tau > 0.5$, since
  $$\nu =c_s^{ * 2}\left(\tau^*-\frac{1}{2}\right)\frac{\Delta x^2}{\Delta t}$$

with the star * representing a non-dimensional/lattice quantity. So the above-mentioned definition would violate this rule, as well as the consistency check between the kinematic viscosity and space and time discretisations (I'm also assuming that system.time_step is $\Delta t$ and agrid is $\Delta x$). What am I missing here? Did I misinterpret these parameters?

• As for the interaction between particles and particles-walls, I could not find a potential that can model an elastic collision. The Weeks–Chandler–Andersen (WCA) potential, a special case of Lennard-Jones interaction, looks like a promising alternative, since it resembles a hard-sphere potential . Is there a better option?
• In an old user guide of espresso from 2016, back when part of the code was written in Tcl, constraints could be defined with a reflecting argument. Was this feature translated to the current version of espresso?
• When adding new particles to the system, is it possible to specify their size/radius? It doesn't appear listed in the documentation as an argument of the part.add function.

Thank you in advance for your time.

[jngrad]

Hi!

• The reduced unit system needs to be defined by the user, preferably with a tool like pint. Then dimensioned quantities like the particle position or the electrostatics prefactor (equation 1 in the user guide chapter on electrostatics) can be passed as arguments to ESPResSo functions, in the form of a dimensionless quantity, e.g. via position.to('reduced_length').magnitude. The pyMBE project provides a pint-based API that forwards dimensioned quantities to an ESPResSo system as dimensionless quantities (see e.g. pyMBE.py).
• Krüger et al. 2017 (doi:10.1007/978-3-319-44649-3) defines the kinematic shear viscosity as $\nu = c_s^2(\tau -\frac{\Delta t}{2})$ (eqn. 3.5) with $\tau$ the collision operator relaxation rate and $c_s^2 = \frac{1}{3}\left(\frac{\Delta x}{\Delta t}\right)^2$ the squared speed of sound. Here $\Delta x$ and $\Delta t$ represent agrid resp. tau, which are unity in LB units (thus $\tau$ is indeed larger than 0.5 in LB units), and are 1.0 resp 0.01 in MD units in the code example, while system.time_step is the MD time step, and is here 0.01 in MD units.
• The simulation setup you are describing sounds like hard spheres, sometimes also called hard disks. They are defined with a discontinuous potential $\Phi(\mathbf{r}_1,\mathbf{r}_2)$ outlined below, which is not differentiable and thus cannot be integrated via the velocity-Verlet scheme. This type of simulation requires custom integrators, such as those described in Ben Tsou, Kevin Wayne, Molecular Dynamics Simulation of Hard Spheres, 2004. In ESPResSo, we use the Lennard-Jones potential to simulate particle dynamics in the NVE ensemble. Collisions will redistribute momentum progressively over multiple time steps. Time discretization leads to precision loss and can cause tiny fluctuations in the total energy of the system. While there are corrective methods (see e.g. section 4.2.2 Energy and momentum conserving discretization in Schiebl and Romero 2021, doi:10.1007/s00466-020-01971-6), I'm not familiar with these techniques and cannot comment on them. In practice, many ESPResSo simulations use a thermostat to simulate the NVT ensemble.

$$\Phi(\mathbf{r}_1,\mathbf{r}_2)= \begin{cases} \,\,\,0 & \text{if} |\mathbf{r}_1-\mathbf{r}_2| \geq \sigma\\\ \infty & \text{if}|\mathbf{r}_1-\mathbf{r}_2| < \sigma \end{cases}$$

• I'm not aware of any reflecting boundary in the Python versions of ESPResSo.
• In ESPResSo, particles do not have a size or radius. Features that rely on a concept of radius have to define their own metric. For example, 3D visualizers use the Lennard-Jones parameter sigma as a proxy for the particle diameter. Monte Carlo schemes expect an exclusion_radius_per_type dictionary that maps particle types to user-defined hard sphere radii. Large nanoparticles are modeled using raspberries, i.e. a constellation of virtual particles whose forces are summed onto a central particle; these are used in lattice-Boltzmann simulations, because LB particle coupling relies on a point-coupling scheme. The raspberry will be coupled to the fluid on its entire surface (or volume if its interior is also filled with virtual particles), which is mathematically equivalent to treating particles as spherical objects rather than point-like objects.

[rdapena]

Thank you for your response. I'll go over the available potentials and integrators carefully. So far, I've been revising on the topic of units in LB following Krüger et al. 2017. But I'm still not sure if I undestand how espresso treats units. Going from what you said, all the parameters indicated in the definition of the fluid (LBFluidWalberla()) are given in MD units, which I assumed before were physical units (i.e dimensional), but they also seem to be non-dimensional units (correct me if I'm wrong). Are MD and LB units dimensionalised by the same conversion factors? If not, how are they related?

[jngrad]

We use MD units everywhere in the Python interface, including in the LB code (this is not the case in the C++ LB code). Users need to carry out the appropriate conversion in their Python scripts, since our API only accepts dimensionless quantities. To convert a MD velocity to LB units, you would like to divide it by the LB characteristic speed $\left(\Delta x / \Delta t \right)$. Here are the conversion factors for all physical properties:

espresso/src/script_interface/walberla/LBFluid.cpp

Lines 166 to 173 in 03033c4

	m_conv_dist = 1. / agrid;
	m_conv_visc = Utils::int_pow<1>(tau) / Utils::int_pow<2>(agrid);
	m_conv_energy = Utils::int_pow<2>(tau) / Utils::int_pow<2>(agrid);
	m_conv_dens = Utils::int_pow<3>(agrid);
	m_conv_speed = Utils::int_pow<1>(tau) / Utils::int_pow<1>(agrid);
	m_conv_press = Utils::int_pow<2>(tau) * Utils::int_pow<1>(agrid);
	m_conv_force = Utils::int_pow<2>(tau) / Utils::int_pow<1>(agrid);
	m_conv_force_dens = Utils::int_pow<2>(tau) * Utils::int_pow<2>(agrid);

We multiply MD units by these factors to yield LB units, and divide LB units by these factors to yield MD units.

After the summer school ends, I'll improve the section on units conversion in the LB theory tutorial.

[rdapena]

Thank you for pointing me to the file with the conversion factors. Going back to the sample I mentioned in my original question, the parameters of the fluid in lattice units would then be: agrid=1.0, density=1.0, kinematic_viscosity=0.01, tau=1 (using m_conv_dist=1, m_conv_visc= 0.01 and m_conv_dens=1).

The condition of $\tau$ being greater than 0.5 in LB units is then met, but according to Krüger et al. 2017 (doi:10.1007/978-3-319-44649-3) the kinematic viscosity in lattice units must also verify:

$$\nu^* =c_s^{ * 2}\left(\tau^{ *}-\frac{1}{2}\right)$$

(follows from eqn. 7.19 or eqn.7.22, or non-dimensionalizing eqn. 7.14), which does not hold for the aforementioned parameters in lattice units for a $c_s^{ * 2}=\frac{1}{3}$. What am I missing?

[jngrad]

That's a good point. The kinematic viscosity should be fully determined when the agrid, tau and time_step parameters are provided. According to the ESPResSo source code, the kinematic_viscosity parameter only affects the $\omega_{\nu}$ relaxation rate (eq. 10.34) and the expression of the pressure tensor. We use the TRT method with parameters $\omega_{e} = \omega^{+} = \omega_{\nu}$ the bulk viscosity relaxation rate, even relaxation rate and shear viscosity relaxation rate, respectively, and $\omega^{-}$ the odd relaxation rate (section 10.6.1), which in ESPResSo is derived from $\omega_{\nu}$.

@RudolfWeeber I think our documentation could better explain the role of the kinematic_viscosity parameter. We had discussions about it and the role of the waLBerla conversion factor between bulk and shear relaxation rates, but I never wrote our conclusions down in the user guide.

[rdapena]

I see, I thought ESPResSo used the BGK model with a single relaxation rate. I'll take a closer look at the TRT method then. Thank you for all the help!