diff --git a/examples/mutation/Allclean.sh b/examples/mutation/Allclean.sh new file mode 100755 index 0000000000..b6c52426b9 --- /dev/null +++ b/examples/mutation/Allclean.sh @@ -0,0 +1,4 @@ +rm -rf vtk +rm -rf hdf5 +rm log.* +rm *.avi diff --git a/examples/mutation/atom.in b/examples/mutation/atom.in new file mode 100644 index 0000000000..db87140d6e --- /dev/null +++ b/examples/mutation/atom.in @@ -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 + diff --git a/examples/mutation/inputscript.nufeb b/examples/mutation/inputscript.nufeb new file mode 100644 index 0000000000..6c7ac83ce8 --- /dev/null +++ b/examples/mutation/inputscript.nufeb @@ -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 + diff --git a/install.sh b/install.sh index 36170204e9..744a19e69d 100755 --- a/install.sh +++ b/install.sh @@ -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 @@ -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 diff --git a/src/MUTATION/fix_growth_mutant.cpp b/src/MUTATION/fix_growth_mutant.cpp new file mode 100644 index 0000000000..b63b3bba4b --- /dev/null +++ b/src/MUTATION/fix_growth_mutant.cpp @@ -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 +#include +#include + +#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(); +} diff --git a/src/MUTATION/fix_growth_mutant.h b/src/MUTATION/fix_growth_mutant.h new file mode 100644 index 0000000000..dc5522a557 --- /dev/null +++ b/src/MUTATION/fix_growth_mutant.h @@ -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: +*/ diff --git a/src/MUTATION/fix_mutate.cpp b/src/MUTATION/fix_mutate.cpp new file mode 100644 index 0000000000..52d71fa749 --- /dev/null +++ b/src/MUTATION/fix_mutate.cpp @@ -0,0 +1,100 @@ +/* ---------------------------------------------------------------------- + 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_mutate.h" + +#include +#include "atom.h" +#include "error.h" +#include "modify.h" +#include "group.h" +#include "update.h" +#include "atom_masks.h" +#include "random_park.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixMutate::FixMutate(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 6) + error->all(FLERR, "Illegal fix mutation/mutate command"); + + imutant = group->find(arg[3]); + if (imutant < 0) + error->all(FLERR, "Can't find group in fix mutation/mutate"); + imutant = 1 | group->bitmask[imutant]; + + prob = utils::numeric(FLERR,arg[4],true,lmp); + if (prob > 1 || prob < 0) + error->all(FLERR, "Illegal fix mutation/mutate command"); + + seed = utils::inumeric(FLERR,arg[5],true,lmp); + if (seed < 0) + error->all(FLERR, "Illegal fix mutation/mutate command"); + + // Random number generator, same for all procs + random = new RanPark(lmp, seed); +} + +/* ---------------------------------------------------------------------- */ + +int FixMutate::setmask() +{ + int mask = 0; + mask |= BIOLOGY_NUFEB; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +int FixMutate::modify_param(int narg, char **arg) +{ + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "nevery") == 0) { + nevery = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + if (nevery <= 0) error->all(FLERR,"Illegal fix_modify command"); + iarg += 2; + } else { + error->all(FLERR, "Illegal fix_modify command"); + } + } + return iarg; +} + +/* ---------------------------------------------------------------------- */ + +void FixMutate::biology_nufeb() +{ + if (update->ntimestep % nevery) return; + compute(); +} + +/* ---------------------------------------------------------------------- */ + +void FixMutate::compute() +{ + int *mask = atom->mask; + + for (int i = 0; i < atom->nlocal; i++) { + if (mask[i] & groupbit) { + if (random->uniform() < prob) { + mask[i] = imutant; + }MUT + } + } +} diff --git a/src/MUTATION/fix_mutate.h b/src/MUTATION/fix_mutate.h new file mode 100644 index 0000000000..eb5fb4481e --- /dev/null +++ b/src/MUTATION/fix_mutate.h @@ -0,0 +1,49 @@ +/* -*- 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/mutate,FixMutate) + +#else + +#ifndef LMP_FIX_MUTATE_H +#define LMP_FIX_MUTATE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixMutate : public Fix { + public: + FixMutate(class LAMMPS *, int, char **); + ~FixMutate() {} + int modify_param(int, char **); + + void init() {} + int setmask(); + void biology_nufeb(); + void compute(); + + private: + class RanPark *random; + + int imutant; // group index of target species + double prob; // mutation probability + int seed; // random seed +}; + +} + +#endif +#endif diff --git a/src/Makefile b/src/Makefile index 2d0c844f1b..6abd220285 100644 --- a/src/Makefile +++ b/src/Makefile @@ -141,6 +141,7 @@ PACKAGE = \ phonon \ nufeb \ hdf5 \ + mutation \ plasmid # NUFEB-specific # NOTE: the last four packages must remain at the end since @@ -204,6 +205,7 @@ PACKMOST = \ phonon \ nufeb \ hdf5 \ + mutation \ plasmid # NOTE ^^^^^: the last three packages must remain at the end since