Molecules diffuse, even when diffusion is set to 0

[iquasere]

I'm trying to create local scarcity of nutrients, using this script

# NUFEB simulation with methanogens and GAC floor

units si                                    # using si units
atom_style      coccus                      # using nufeb atom style

# map array: find atoms using indices sort 1000 5.0e-6: sort every 1000 steps with 5.0e-6 binsize
atom_modify     map array sort 1000 1e-6

# periodic boundaries in x and y fixed boundary in z
boundary        ff ff ff

# forces between local and ghost atoms are computed in each processor without communication
newton          off

processors      * * *                      # processor grid

# communicate velocities for ghost atoms
comm_modify     vel yes

# guarantee that enough atoms are communicated to correctly compute
comm_modify     cutoff 2e-6

read_data      gac_pure_culture.in

# Shift the lattice grid 0.5 unit so that atom can be created in the center of the grid
lattice sc 1e-6 origin 0.5 0.5 0.5
region reg block 0 20 0 20 1 20

create_atoms 1 random 80 31324 reg

# create the GAC floor - this is influenced by the lattice command
region gac_floor block 0 20 0 20 0 1
create_atoms 3 region gac_floor

region outside_gac_floor block 0 20 0 20 0 1

set type 1 density 150
set type 1 diameter 1.3e-6
set type 1 mass 1.725e-16
set type 1 outer_diameter 1.3e-6
set type 1 outer_density 30

# define attributes for type 3, set gac diameter to 1e-6, make it consistent with lattice size
set type 3 density 150
set type 3 mass 1.725e-16
set type 3 diameter 1e-6
set type 3 outer_diameter 1e-6
set type 3 outer_density 30

# set the groups
group met type 1
group eps type 2
group gac type 3

# setting neighbour skin distance and style
neighbor        7e-7 bin

# rebuild neighbour list if any atom had moved more than half the skin distance
neigh_modify    check yes

# select grid style
grid_style      nufeb/chemostat  4 h2 gco2 co2 ch4 2e-6

# set substrates initial concentration
grid_modify     set h2      nn nn nn 10 10
grid_modify     set gco2    nn nn nn 10 10
grid_modify     set co2     nn nn nn 0 0
grid_modify     set ch4     nn nn nn 0 0

# define pair styles
pair_style hybrid gran/hooke/history 1e-4 NULL 1e-5 NULL 0.0 1 lj/cut 2.5e-6
pair_coeff * * gran/hooke/history
pair_coeff *2 3 lj/cut 1.0e-20 1.0e-6 2.5e-6

# NVE integration with maximum distance limit -> only update position for non gac atoms
fix nve1 met nve/limit 1e-8
fix nve2 eps nve/limit 1e-8

# monod reaction fixes -> u max = 0.69 / h = 1.92e-4 / s
fix monod_met met nufeb/growth/methanogen h2 6e-6 co2 2e-4 ch4 gco2 growth 1.92e-4 yield 1.6 decay 0 epsyield 1 epsdens 30 gco2_flag 1
#fix monod_eps eps nufeb/growth/eps h2 decay 0

# diffusion reaction fixes
fix diff_h2  all nufeb/diffusion_reaction   h2      0
fix diff_ch4 all nufeb/diffusion_reaction   gco2    0
fix diff_co2 all nufeb/diffusion_reaction   co2     0
fix diff_ch4 all nufeb/diffusion_reaction   ch4     0

# biological model fixes -> apply division to met only
fix div met nufeb/division/limited 1.36e-6 30 1234 region_blocked gac_floor
fix eps_ext met nufeb/eps_extract 2 eps 1.3 30 5678

# mechanical model fixes
fix wall all wall/gran hooke/history 1e-3 NULL 1e-4 NULL 0 0 zplane 0.0 8e-5
fix eps_adh all nufeb/adhesion/eps eps 1e-6
fix vis all viscous 1e-5
# fix ID atomsID nufeb/shear shearRate viscosity +x (or -x or +y or -y)
fix s1 met nufeb/shear 5e-7 0.01 +x
fix s2 eps nufeb/shear 5e-7 0.01 +x

# pressure computation
compute vol all nufeb/volume
compute ke all ke
variable one equal 1.0
compute press all pressure NULL pair vol v_one
variable press equal "(c_ke + c_press) / (3.0 * c_vol)"
variable mass equal "mass(all)"

variable nmet equal "count(met)"
variable neps equal "count(eps)"
compute mycon all nufeb/ave_conc

# file output
shell mkdir vtk_gac_sheet
dump 1 all vtk 10 vtk_gac_sheet/dump*.vtu id type diameter
dump 2 all grid/vtk 10 vtk_gac_sheet/dump_%_*.vti con

# thermo output
thermo_style custom step atoms v_press v_mass v_nmet v_neps c_mycon[*]
thermo 1
thermo_modify lost ignore flush yes

# issue run command
run_style nufeb diffdt 1e-4 difftol 1e-12 pairdt 2 pairtol 0 pairmax 500 diffmax 50000
timestep 1000
run 500

and pointing at the lines

fix diff_h2  all nufeb/diffusion_reaction   h2      0
fix diff_ch4 all nufeb/diffusion_reaction   gco2    0
fix diff_co2 all nufeb/diffusion_reaction   co2     0
fix diff_ch4 all nufeb/diffusion_reaction   ch4     0

I thought I would have no diffusion of all metabolites. However, I still get quite a lot of it.

https://user-images.githubusercontent.com/16226371/220677536-36b5f29a-0593-4e7b-bca5-8d6c183f68e6.mp4

In this example, I would expect the methane produced to be kept at the bottom, however it fills everywhere at almost the same rate. Shouldn't it be kept at the bottom?

[shelllbw]

This is a bug. I will fix that. Instead of set to 0, you can remove fix diffusion command to avoid any diffusion for now. Regarding keeping methane at the bottom, I think if diffusion coeff > 0, the substrate can always diffuse to the entire box with sufficient long time (e.g, biological timestep)

[iquasere]

Thank you, Bowen! Didn't think of the timestep, since I'm only printing every 10th, it might be diffusing too fast for that.

[iquasere]

@shelllbw , seems to me this is still not fixed, is that correct? Does that also impact creating different local pHs, i.e., the pH is always a global homogeneous value?

I think if diffusion coeff > 0, the substrate can always diffuse to the entire box with sufficient long time

If I set the diffusion coeff to a small enough value (when it is fixed), it should be in line with biological processes, no?

[shelllbw]

seems to me this is still not fixed, is that correct?

The averaging is switched off in the previous commit. Is your code update-to-date?

Does that also impact creating different local pHs, i.e., the pH is always a global homogeneous value?

There is no pH in the system, the fix has not been migrated to the -dev version. fix nufeb/growth/energy and pH would be separate fixes.

If I set the diffusion coeff to a small enough value (when it is fixed), it should be in line with biological processes, no?

Sorry, I don't understand the question. Could you explain more?

[iquasere]

The averaging is switched off in the previous commit. Is your code update-to-date?

You are absolutely correct, I forgot to pull the code on the remote server, only did it in my local laptop! It works perfectly now!

There is no pH in the system, the fix has not been migrated to the -dev version. fix nufeb/growth/energy and pH would be separate fixes.

I'll keep it out for now, then. Makes sense to keep them separate, and I might hack something now that diffusion can be used to create local values of metabolites...

Sorry, I don't understand the question. Could you explain more?

Chemical species diffusion is much quicker than biological processes, as you said. But I want to simulate the sequestration of substrates by nanomaterials, and setting the diffusion to much lower values (justified by the effects of the nanoparticles) might allow to simulate this. Also, mechanisms of quorum sensing might also be simulated by release of qs signals with much lower diffusion coefficients.

You've given me a lot of tools to work here, is what I mean. Thank you very much!

[shelllbw]

Now I think I understand better. The consistency between biological and diffusion processes would be the timesteps rather than diffusion coeff. Say the timestep for bio process is 100s and the timestep for diff is 1e-4s. In theory, we need to run diffusion for 1e6 steps in order to make the two processes consistent. But in the code, we use the mechanism of pseudo-steady state and frozen state to speed up the simulation for both open and closed systems.