This issue was already happening with previous versions of NUFEB, could it be a problem of using a Lennard-Jones with atoms of type coccus?
atom.in
NUFEB Simulation
0 atoms
3 atom types
0 2e-05 xlo xhi
0 2e-05 ylo yhi
0 2e-05 zlo zhi
- A small 8 atoms sphere in the middle of the simulation puts the pressure at
2.9760175e+17.
# 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 atom.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
region gac_sphere1 sphere 10 10 10 1
create_atoms 3 region gac_sphere1
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 HET 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
# defining grid sytle, substrate names, and grid size
grid_style nufeb/chemostat 4 sub o2 no2 no3 4e-6
# set diffusion boundary conditions and initial concentrations (liquid:kg/m3)
grid_modify set sub pp pp nd 1e-4
grid_modify set o2 pp pp nd 1e-4
grid_modify set no2 pp pp nd 1e-4
grid_modify set no3 pp pp nd 1e-4
# 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 *2 *2 gran/hooke/history
pair_coeff 3 3 lj/cut 1.0 1.0e-6 2.5e-6
pair_coeff *2 3 lj/cut 1.0e-6 1.0e-6 2.5e-6
# NVE integration with maximum distance limit -> only update position for non gac atoms
fix nve1 HET nve/limit 1e-8
fix nve2 EPS nve/limit 1e-8
# heterotrophs growth
fix growth_het HET nufeb/growth/het sub 3.5e-5 o2 0 no2 0 no3 0 &
growth 0.00028 yield 0.61 decay 0.0 epsyield 0.18 anoxic 0.0 epsdens 30
# heterotrophs division, division diameter: 1.36e-6
fix div HET nufeb/division/coccus 1.36e-6 1234
# EPS secretion from heterotrophs
fix eps_ext HET nufeb/eps_extract 2 EPS 1.3 30 2345
# diffusion reaction fixes
fix diff_sub all nufeb/diffusion_reaction sub 1.6e-9
fix div HET nufeb/division/coccus 1.36e-6 1234
fix eps_ext HET 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
# 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(HET)"
variable neps equal "count(EPS)"
compute mycon all nufeb/ave_conc
# file output
shell mkdir vtk_gac_rigids
dump 1 all vtk 10 vtk_gac_rigids/dump*.vtu id type diameter
dump 2 all grid/vtk 10 vtk_gac_rigids/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 300
- The weirdest one. Simply creating random atoms of type 3 puts the pressure at
-5.2855696e+15.
# 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 atom.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_atoms 3 random 80 31325 reg
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 HET 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
# defining grid sytle, substrate names, and grid size
grid_style nufeb/chemostat 4 sub o2 no2 no3 4e-6
# set diffusion boundary conditions and initial concentrations (liquid:kg/m3)
grid_modify set sub pp pp nd 1e-4
grid_modify set o2 pp pp nd 1e-4
grid_modify set no2 pp pp nd 1e-4
grid_modify set no3 pp pp nd 1e-4
# 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 *2 *2 gran/hooke/history
pair_coeff 3 3 lj/cut 1.0 1.0e-6 2.5e-6
pair_coeff *2 3 lj/cut 1.0e-6 1.0e-6 2.5e-6
# NVE integration with maximum distance limit -> only update position for non gac atoms
fix nve1 HET nve/limit 1e-8
fix nve2 EPS nve/limit 1e-8
# heterotrophs growth
fix growth_het HET nufeb/growth/het sub 3.5e-5 o2 0 no2 0 no3 0 &
growth 0.00028 yield 0.61 decay 0.0 epsyield 0.18 anoxic 0.0 epsdens 30
# heterotrophs division, division diameter: 1.36e-6
fix div HET nufeb/division/coccus 1.36e-6 1234
# EPS secretion from heterotrophs
fix eps_ext HET nufeb/eps_extract 2 EPS 1.3 30 2345
# diffusion reaction fixes
fix diff_sub all nufeb/diffusion_reaction sub 1.6e-9
fix div HET nufeb/division/coccus 1.36e-6 1234
fix eps_ext HET 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
# 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(HET)"
variable neps equal "count(EPS)"
compute mycon all nufeb/ave_conc
# file output
shell mkdir vtk_gac_rigids
dump 1 all vtk 10 vtk_gac_rigids/dump*.vtu id type diameter
dump 2 all grid/vtk 10 vtk_gac_rigids/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 300
- Mountains of immovable atoms of type 3 put the pressure at
-nan.
# 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 atom.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
region gac_floor block 0 20 0 20 0 1
create_atoms 3 region gac_floor
region gac_limits block 0 20 0 20 0 20
variable xx internal 0.0
variable yy internal 0.0
variable zz internal 0.0
variable v equal "v_zz < sin(sin(5e5*v_xx))*sin(cos(5e5*v_yy))/10e4"
create_atoms 3 region gac_limits var v set x xx set y yy set z zz
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 HET 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
# defining grid sytle, substrate names, and grid size
grid_style nufeb/chemostat 4 sub o2 no2 no3 4e-6
# set diffusion boundary conditions and initial concentrations (liquid:kg/m3)
grid_modify set sub pp pp nd 1e-4
grid_modify set o2 pp pp nd 1e-4
grid_modify set no2 pp pp nd 1e-4
grid_modify set no3 pp pp nd 1e-4
# 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 *2 *2 gran/hooke/history
pair_coeff 3 3 lj/cut 1.0 1.0e-6 2.5e-6
pair_coeff *2 3 lj/cut 1.0e-6 1.0e-6 2.5e-6
# NVE integration with maximum distance limit -> only update position for non gac atoms
fix nve1 HET nve/limit 1e-8
fix nve2 EPS nve/limit 1e-8
# heterotrophs growth
fix growth_het HET nufeb/growth/het sub 3.5e-5 o2 0 no2 0 no3 0 &
growth 0.00028 yield 0.61 decay 0.0 epsyield 0.18 anoxic 0.0 epsdens 30
# heterotrophs division, division diameter: 1.36e-6
fix div HET nufeb/division/coccus 1.36e-6 1234
# EPS secretion from heterotrophs
fix eps_ext HET nufeb/eps_extract 2 EPS 1.3 30 2345
# diffusion reaction fixes
fix diff_sub all nufeb/diffusion_reaction sub 1.6e-9
fix div HET nufeb/division/coccus 1.36e-6 1234
fix eps_ext HET 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
# 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(HET)"
variable neps equal "count(EPS)"
compute mycon all nufeb/ave_conc
# file output
shell mkdir vtk_gac_rigids
dump 1 all vtk 10 vtk_gac_rigids/dump*.vtu id type diameter
dump 2 all grid/vtk 10 vtk_gac_rigids/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 300
This issue was already happening with previous versions of NUFEB, could it be a problem of using a Lennard-Jones with atoms of type coccus?
atom.in2.9760175e+17.-5.2855696e+15.-nan.