Can coccus work with RIGID package?

[iquasere]

I have been trying to create RIGID structures with the coccus atom style. The structures are created just fine, but at the first timestep they disappear. This is the simplest script I could put together to simulate that behaviour:

• minimal.lmp

# 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      minimal.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
group group_gac_sphere1 region gac_sphere1

# 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

# fix ID group-ID style bodystyle args keyword values
fix sphere1 group_gac_sphere1 rigid/nve single reinit yes force 1 on on on force 1 on on on

# 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

• minimal.in

NUFEB Simulation 

       0 atoms
       3 atom types

       0   2e-05  xlo xhi
       0   2e-05  ylo yhi
       0   2e-05  zlo zhi

The "gac" molecules (type 3 atoms) immediately disappear. I think this might be something LAMMPS and not purely NUFEB related, but since it's such an old version of LAMMPS, I'm sure Axel would just tell me to piss off.

I've tried for three days but haven't managed to put it to work. Trying with just the "gac" atoms alone has similar results, but ends immediately with a "Segmentation fault", which is likely because the simulation lost all atoms. The most likely explanation here is the extreme forces caused by the atoms so close together, with the fix rigid/nve command they cannot expand, and instead just pick a direction to move the combined force away...

I might try to do this with a "hybrid" type of atoms, but wanted to know if this has an easy fix before revamping my scripts...

[shelllbw]

We've chatted this issue, just want to reply here to make a note. The short answer is yes. The reason of atoms disappear from system is because you use a very large pairdt (physical timestep), so that atoms move outside the domain in a single iteration. The solution is to reduce the pairdt.

In your case, however, since you are coupling physical processes with biological processes and the system is non-steady state due to atom motion, make sure the physical and biological timesteps are consistent.
For example, if pairdt = 0.1s and biodt/timestep = 1000s, then pairmax should be 10000

[iquasere]

Changing the line run_style nufeb diffdt 1e-4 difftol 1e-12 pairdt 2 pairtol 0 pairmax 500 diffmax 50000
to run_style nufeb diffdt 1e-4 difftol 1e-12 pairdt 1e-3 pairtol 0 pairmax 1000e3 diffmax 50000
doesn't seem to impact the simulation. It does however help to see the problem.

For the first time, the little rigid ball is kept in the first timestep, however, the simulation becomes unstable afterward. Also, it's strange that the pressure, which was put back into a normal value by your suggestion on the espilon value of the Lennard-Jones (I put it at 1e-24), is back to -nan.