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_gacrigids/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.