Skip to content

Commit

Permalink
Migrate mutation code from NUFEB:761bc36
Browse files Browse the repository at this point in the history 
* Package: MUTATION

* Src files: fix_growth_mutant fix_mutate

* Example: example/mutation

* Update install.sh and src/Makefile for the new package

To use the package:
$./install.sh --enable-mutation --enable-nufeb
cd to example/mutation
$mpirun -np 4 ../../nufeb_mpi -in inputscript.nufeb
  • Loading branch information
@shelllbw
shelllbw committed Apr 28, 2024
1 parent 6f83d2c commit 0d703b6
Show file tree
Hide file tree
Showing 9 changed files with 488 additions and 0 deletions.
4 changes: 4 additions & 0 deletions examples/mutation/Allclean.sh
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
rm -rf vtk
rm -rf hdf5
rm log.*
rm *.avi
52 changes: 52 additions & 0 deletions examples/mutation/atom.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
NUFEB Simulation

40 atoms
2 atom types

0.000000e-04 1e-04 xlo xhi
0.000000e-04 0.4e-04 ylo yhi
0.000000e-04 1e-04 zlo zhi

Atoms

1 1 1.0e-6 150 0.5e-5 0.5e-5 1e-6 1.0e-6
2 1 1.0e-6 150 1.5e-5 0.5e-5 1e-6 1.0e-6
3 1 1.0e-6 150 2.5e-5 0.5e-5 1e-6 1.0e-6
4 1 1.0e-6 150 3.5e-5 0.5e-5 1e-6 1.0e-6
5 1 1.0e-6 150 4.5e-5 0.5e-5 1e-6 1.0e-6
6 1 1.0e-6 150 5.5e-5 0.5e-5 1e-6 1.0e-6
7 1 1.0e-6 150 6.5e-5 0.5e-5 1e-6 1.0e-6
8 1 1.0e-6 150 7.5e-5 0.5e-5 1e-6 1.0e-6
9 1 1.0e-6 150 8.5e-5 0.5e-5 1e-6 1.0e-6
10 1 1.0e-6 150 9.5e-5 0.5e-5 1e-6 1.0e-6
11 1 1.0e-6 150 0.5e-5 1.5e-5 1e-6 1.0e-6
12 1 1.0e-6 150 1.5e-5 1.5e-5 1e-6 1.0e-6
13 1 1.0e-6 150 2.5e-5 1.5e-5 1e-6 1.0e-6
14 1 1.0e-6 150 3.5e-5 1.5e-5 1e-6 1.0e-6
15 1 1.0e-6 150 4.5e-5 1.5e-5 1e-6 1.0e-6
16 1 1.0e-6 150 5.5e-5 1.5e-5 1e-6 1.0e-6
17 1 1.0e-6 150 6.5e-5 1.5e-5 1e-6 1.0e-6
18 1 1.0e-6 150 7.5e-5 1.5e-5 1e-6 1.0e-6
19 1 1.0e-6 150 8.5e-5 1.5e-5 1e-6 1.0e-6
20 1 1.0e-6 150 9.5e-5 1.5e-5 1e-6 1.0e-6
21 1 1.0e-6 150 0.5e-5 2.5e-5 1e-6 1.0e-6
22 1 1.0e-6 150 1.5e-5 2.5e-5 1e-6 1.0e-6
23 1 1.0e-6 150 2.5e-5 2.5e-5 1e-6 1.0e-6
24 1 1.0e-6 150 3.5e-5 2.5e-5 1e-6 1.0e-6
25 1 1.0e-6 150 4.5e-5 2.5e-5 1e-6 1.0e-6
26 1 1.0e-6 150 5.5e-5 2.5e-5 1e-6 1.0e-6
27 1 1.0e-6 150 6.5e-5 2.5e-5 1e-6 1.0e-6
28 1 1.0e-6 150 7.5e-5 2.5e-5 1e-6 1.0e-6
29 1 1.0e-6 150 8.5e-5 2.5e-5 1e-6 1.0e-6
30 1 1.0e-6 150 9.5e-5 2.5e-5 1e-6 1.0e-6
31 1 1.0e-6 150 0.5e-5 3.5e-5 1e-6 1.0e-6
32 1 1.0e-6 150 1.5e-5 3.5e-5 1e-6 1.0e-6
33 1 1.0e-6 150 2.5e-5 3.5e-5 1e-6 1.0e-6
34 1 1.0e-6 150 3.5e-5 3.5e-5 1e-6 1.0e-6
35 1 1.0e-6 150 4.5e-5 3.5e-5 1e-6 1.0e-6
36 1 1.0e-6 150 5.5e-5 3.5e-5 1e-6 1.0e-6
37 1 1.0e-6 150 6.5e-5 3.5e-5 1e-6 1.0e-6
38 1 1.0e-6 150 7.5e-5 3.5e-5 1e-6 1.0e-6
39 1 1.0e-6 150 8.5e-5 3.5e-5 1e-6 1.0e-6
40 1 1.0e-6 150 9.5e-5 3.5e-5 1e-6 1.0e-6

115 changes: 115 additions & 0 deletions examples/mutation/inputscript.nufeb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
#-----------------------------------------------------------------------------#
# NUFEB Simulation: Heterotrophic Biofilm Growth #
#-----------------------------------------------------------------------------#


#-----------------------------System Settings---------------------------------#

units si # using si units (m, s, kg)
atom_style coccus # using nufeb coccus atom style
atom_modify map array sort 10 0 # map array - find atoms using indices
# sort every 10

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

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

processors * * 1 # processor grid

comm_modify vel yes # communicate velocities for ghost atoms

read_data atom.in # read atom.in file which defines domain size
# and initial atoms


#--------------------Microbes and Functional Groups-------------------------#

group MUTA type 1 # assign type 1 atoms to MUTA group
group MUTB type 2 # assign type 2 atoms to MUTB group

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

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


#--------------------------Mesh Grid and Substrates--------------------------#

# defining grid sytle, substrate names, and grid size
grid_style nufeb/chemostat 2 sub inh 4e-6

# set diffusion boundary conditions and initial concentrations (liquid:kg/m3)
grid_modify set sub pp pp nd 1e-4
grid_modify set inh pp pp nd 1e-4


#--------------------------Biological Processes-------------------------------#

# mutant_a growth
fix f1 MUTA mutation/growth/mutant sub 1e-5 inh &
growth 0.0001 gamma 0.1 ic50 0.1 nl 0.1

# mutant_b growth
fix f2 MUTB mutation/growth/mutant sub 2e-5 inh &
growth 0.0002 gamma 0.2 ic50 0.2 nl 0.2

# mutate from MUTA to MUTB
fix f3 MUTA mutation/mutate MUTB 0.5 2024

# mutant division, division diameter: 1.36e-6
fix div all nufeb/division/coccus 1.36e-6 2025

#---------------------------Physical Processes--------------------------------#

pair_style gran/hooke/history 1e-4 NULL 1e-5 NULL 0.0 0 # pairwise interaction
pair_coeff * * # between atoms

# pairwise interaction between z-wall and atoms
fix wall all wall/gran hooke/history 1e-3 NULL 1e-4 NULL 0 0 zplane 0.0 1e-04

fix vis all viscous 1e-5 # viscous damping force

fix nve all nve/limit 1e-7 # NVE integration with
# maximum distance limit


#---------------------------Chemical Processes---------------------------------#

fix diff_sub all nufeb/diffusion_reaction sub 1.6e-9 # diffusion reaction for updating
# distribusion of substrate concentration


#--------------------------Computations and Outputs----------------------------#

variable mass equal "mass(all)" # total mass
variable nmuta equal "count(MUTA)" # total # of HET
variable nmutb equal "count(MUTB)" # total # of EPS

#shell mkdir vtk # dump vtk files to /vtk folder
#dump du1 all vtk 10 vtk/dump*.vtu id type diameter # require build NUFEB with vtk option
#dump du2 all grid/vtk 10 vtk/dump_%_*.vti con rea den gro


dump du3 all movie 10 movie.avi type diameter view 80 60 # dump video
dump_modify du3 acolor 1 blue acolor 2 lightgrey framerate 10


# dump hdf5 files to /hdf5 folder
#shell mkdir hdf5 # require build NUFEB with hdf5 option
#dump du3 all nufeb/hdf5 10 dump.h5 id type x y z vx vy vz fx fy fz radius conc reac

thermo_style custom step cpu atoms v_nmuta v_nmutb v_mass # screen and log outputs
thermo 1

#-----------------------------------Run----------------------------------------#

# issue run command, define timesteps for physical (pairdt) and chemical (diffdt) processes
run_style nufeb diffdt 1e-4 difftol 1e-6 pairdt 1e-2 pairtol 1 pairmax 1000 diffmax 5000

timestep 1000 # define biological timesteps (1000s)

run 900 # growing biofilm for 900x1000s

3 changes: 3 additions & 0 deletions install.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ do
vtk_hdf=$((vtk_hdf+1))
elif [ $var == "--enable-misc" ]; then continue
elif [ $var == "--enable-plasmid" ]; then continue
elif [ $var == "--enable-mutation" ]; then continue
elif [ $var == "--gpu" ]; then continue
elif [ $var == "--static" ]; then continue
elif [ $var == "--shared" ]; then continue
Expand Down Expand Up @@ -61,6 +62,8 @@ do
make yes-misc
elif [ $var == "--enable-plasmid" ]; then
make yes-plasmid
elif [ $var == "--enable-mutation" ]; then
make yes-mutation
elif [ $var == "--gpu" ]; then
make yes-kokkos
fi
Expand Down
111 changes: 111 additions & 0 deletions src/MUTATION/fix_growth_mutant.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "fix_growth_mutant.h"

#include <cstdio>
#include <cstring>
#include <cmath>

#include "atom.h"
#include "error.h"
#include "grid.h"
#include "group.h"
#include "grid_masks.h"
#include "math_const.h"

using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;

/* ---------------------------------------------------------------------- */

FixGrowthMutant::FixGrowthMutant(LAMMPS *lmp, int narg, char **arg) :
FixGrowth(lmp, narg, arg)
{
if (narg < 6)
error->all(FLERR, "Illegal fix mutation/growth/mutant command");

if (!grid->chemostat_flag)
error->all(FLERR, "fix mutation/growth/mutant requires grid_style nufeb/chemostat");

isub = -1;
iinh = -1;

sub_affinity = 0.0;

growth = 0.0;
ic50 = 1.0;
nl = 1.0;

isub = grid->find(arg[3]);
if (isub < 0)
error->all(FLERR, "Can't find substrate name");
sub_affinity = utils::numeric(FLERR,arg[4],true,lmp);

iinh = grid->find(arg[5]);
if (iinh < 0)
error->all(FLERR, "Can't find inhibitor name");

int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg], "growth") == 0) {
growth = utils::numeric(FLERR,arg[iarg+1],true,lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "gamma") == 0) {
gamma = utils::numeric(FLERR,arg[iarg+1],true,lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "ic50") == 0) {
ic50 = utils::numeric(FLERR,arg[iarg+1],true,lmp);
iarg += 2;
} else if (strcmp(arg[iarg], "nl") == 0) {
nl = utils::numeric(FLERR,arg[iarg+1],true,lmp);
iarg += 2;
} else {
error->all(FLERR, "Illegal fix mutation/growth/mutant command");
}
}
}

/* ---------------------------------------------------------------------- */

void FixGrowthMutant::update_cells()
{
double **conc = grid->conc;
double **reac = grid->reac;
double **dens = grid->dens;

for (int i = 0; i < grid->ncells; i++) {
if (grid->mask[i] & GRID_MASK) {
double tmp1 = growth * (conc[isub][i] / (sub_affinity + conc[isub][i])) / (1 + pow((conc[iinh][i] / ic50), nl));

reac[isub][i] += -gamma * tmp1 * dens[igroup][i];
reac[iinh][i] += 0;
}
}
}

/* ---------------------------------------------------------------------- */

void FixGrowthMutant::update_atoms()
{
double **conc = grid->conc;

for (int i = 0; i < grid->ncells; i++) {
double tmp1 = growth * (conc[isub][i] / (sub_affinity + conc[isub][i])) / (1 + pow ((conc[iinh][i] / ic50),nl));

grid->growth[igroup][i][0] = tmp1;
}

update_atoms_coccus();
}
52 changes: 52 additions & 0 deletions src/MUTATION/fix_growth_mutant.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#ifdef FIX_CLASS

FixStyle(mutation/growth/mutant,FixGrowthMutant)

#else

#ifndef LMP_FIX_GROWTH_MUTANT_H
#define LMP_FIX_GROWTH_MUTANT_H

#include "fix_growth.h"

namespace LAMMPS_NS {

class FixGrowthMutant: public FixGrowth {
public:
FixGrowthMutant(class LAMMPS *, int, char **);
virtual ~FixGrowthMutant() {}

virtual void update_atoms();
virtual void update_cells();

protected:
int isub; // substrate group index
int iinh; // inhibitor group index

double growth; // growth rate
double sub_affinity; // ks
double gamma;
double ic50;
double nl;
};

}

#endif
#endif

/* ERROR/WARNING messages:
*/
Loading

0 comments on commit 0d703b6

Please sign in to comment.