PyPI 官网下载 | nufeb_tools-0.0.52.tar.gz

/* ---------------------------------------------------------------------- NUFEB package - A LAMMPS user package for Individual-based Modelling of Microbial Communities Contributing authors: Bowen Li & Denis Taniguchi (Newcastle University, UK) Email: bowen.li2@newcastle.ac.uk & denis.taniguchi@newcastle.ac.uk 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. ------------------------------------------------------------------------- */ #include <math.h> #include <string.h> #include <algorithm> #include <cstdio> #include <iostream> #include <iterator> #include <vector> #include "atom.h" #include "domain.h" #include "error.h" #include "force.h" #include "input.h" #include "lammps.h" #include "math_const.h" #include "memory.h" #include "bio.h" #include "atom_vec_bio.h" #include "fix_bio_kinetics.h" #include "fix_bio_kinetics_monod.h" #include "modify.h" #include "pointers.h" #include "update.h" #include "variable.h" #include "group.h" using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; using namespace std; enum{HET, AOB, NOB, ANA, COM, EPS, DEAD, CYANO, ECW}; /* ---------------------------------------------------------------------- */ FixKineticsMonod::FixKineticsMonod(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { avec = (AtomVecBio *) atom->style_match("bio"); if (!avec) error->all(FLERR, "Fix kinetics requires atom style bio"); if (narg < 3) error->all(FLERR, "Not enough arguments in fix kinetics/growth/monod command"); species = NULL; growrate = NULL; eps_dens = 30; eta_het = 0; suc_exp = 0; gco2_flag = 0; kinetics = NULL; external_gflag = 1; int iarg = 3; while (iarg < narg){ if (strcmp(arg[iarg],"gflag") == 0) { external_gflag = force->inumeric(FLERR, arg[iarg+1]); if (external_gflag != 0 && external_gflag != 1) error->all(FLERR, "Illegal fix kinetics/growth/monod command: gflag"); iarg += 2; } else if (strcmp(arg[iarg],"epsdens") == 0){ eps_dens = force->numeric(FLERR, arg[iarg+1]); if (eps_dens <= 0) error->all(FLERR, "Illegal fix kinetics/growth/monod command: eps_dens cannot be less or equal than zero"); iarg += 2; } else if (strcmp(arg[iarg],"etahet") == 0){ eta_het = force->numeric(FLERR, arg[iarg+1]); if (eta_het < 0) error->all(FLERR, "Illegal fix kinetics/growth/monod command: eta_het cannot be less than zero"); iarg += 2; } else if (strcmp(arg[iarg], "sucexp") == 0){ suc_exp = force->numeric(FLERR, arg[iarg+1]);//suc_exp = input->variable->compute_equal(arg[iarg+1]); if (suc_exp < 0) error->all(FLERR, "Illegal fix kinetics/growth/monod command: suc_exp cannot be less than zero"); iarg += 2; } else if (strcmp(arg[iarg], "gco2flag") == 1){ gco2_flag = force->numeric(FLERR, arg[iarg+1]); if (gco2_flag != 0 && gco2_flag != 1) error->all(FLERR, "Illegal fix kinetics/growth/monod command: gco2_flag"); iarg += 2; } else error->all(FLERR, "Illegal fix kinetics/growth/monod command"); } } /* ---------------------------------------------------------------------- */ FixKineticsMonod::~FixKineticsMonod() { memory->destroy(species); memory->destroy(growrate); } /* ---------------------------------------------------------------------- */ int FixKineticsMonod::setmask() { int mask = 0; mask |= PRE_FORCE; return mask; } /* ---------------------------------------------------------------------- */ void FixKineticsMonod::init() { if (!atom->radius_flag) error->all(FLERR, "Fix requires atom attribute diameter"); // register fix kinetics with this class kinetics = NULL; int nfix = modify->nfix; for (int j = 0; j < nfix; j++) { if (strcmp(modify->fix[j]->style, "kinetics") == 0) { kinetics = static_cast<FixKinetics *>(lmp->modify->fix[j]); break; } } if (kinetics == NULL) lmp->error->all(FLERR, "fix kinetics command is required for running IbM simulation"); bio = kinetics->bio; if (bio->nnu == 0) error->all(FLERR, "fix_kinetics/monod requires Nutrients input"); else if (bio->maintain == NULL) error->all(FLERR, "fix_kinetics/monod requires Maintenance input"); else if (bio->decay == NULL) error->all(FLERR, "fix_kinetics/monod requires Decay input"); else if (bio->ks == NULL) error->all(FLERR, "fix_kinetics/monod requires Ks input"); else if (bio->yield == NULL) error->all(FLERR, "fix_kinetics/monod requires Yield input"); else if (bio->mu == NULL) error->all(FLERR, "fix_kinetics/monod requires Growth Rate input"); nx = kinetics->nx; ny = kinetics->ny; nz = kinetics->nz; if (species == NULL) { species = memory->create(species, atom->ntypes+1, "monod:species"); growrate = memory->create(growrate, atom->ntypes+1, 2, kinetics->ngrids, "monod:growrate"); } //Get computational domain size if (domain->triclinic == 0) { xlo = domain->boxlo[0]; xhi = domain->boxhi[0]; ylo = domain->boxlo[1]; yhi = domain->boxhi[1]; zlo = domain->boxlo[2]; zhi = domain->boxhi[2]; } else { xlo = domain->boxlo_bound[0]; xhi = domain->boxhi_bound[0]; ylo = domain->boxlo_bound[1]; yhi = domain->boxhi_bound[1]; zlo = domain->boxlo_bound[2]; zhi = domain->boxhi_bound[2]; } stepx = (xhi - xlo) / nx; stepy = (yhi - ylo) / ny; stepz = (zhi - zlo) / nz; vol = stepx * stepy * stepz; init_param(); } /* ---------------------------------------------------------------------- initialize parameters for monod growth ------------------------------------------------------------------------- */ void FixKineticsMonod::init_param() { mu = bio->mu; decay = bio->decay; maintain = bio->maintain; yield = bio->yield; ks = bio->ks; ntypes = atom->ntypes; isub = 0; io2 = 0; inh4 = 0; ino2 = 0; ino3 = 0; isuc = 0; ico2 = 0; igco2 = 0; // initialize nutrients for (int nu = 1; nu <= bio->nnu; nu++) { if (strcmp(bio->nuname[nu], "sub") == 0) isub = nu; else if (strcmp(bio->nuname[nu], "o2") == 0) io2 = nu; else if (strcmp(bio->nuname[nu], "nh4") == 0) inh4 = nu; else if (strcmp(bio->nuname[nu], "no2") == 0) ino2 = nu; else if (strcmp(bio->nuname[nu], "no3") == 0) ino3 = nu; else if (strcmp(bio->nuname[nu], "suc") == 0) isuc = nu; else if (strcmp(bio->nuname[nu], "co2") == 0) ico2 = nu; else if (strcmp(bio->nuname[nu], "gco2") == 0) igco2 = nu; } // initialize species for (int i = 1; i <= ntypes; i++) { // take the first three chars from type name; char *name = new char[4]; strncpy(name, bio->tname[i], 3); name[3] = 0; if (strcmp(name, "het") == 0 || strcmp(name, "HET") == 0) { species[i] = HET; if (isub == 0) error->all(FLERR, "het growth requires nutrient 'sub' (substrate) defined in Nutrients section"); if (io2 == 0) error->all(FLERR, "het growth requires nutrient 'o2' defined in Nutrients section"); if (eta_het > 0) { if (ino2 == 0) error->all(FLERR, "het anaerobic growth requires nutrient 'no2' defined in Nutrients section"); if (ino3 == 0) error->all(FLERR, "het anaerobic growth requires nutrient 'no3' defined in Nutrients section"); } } else if (strcmp(name, "aob") == 0 || strcmp(name, "AOB") == 0) { species[i] = AOB; if (inh4 == 0) error->all(FLERR, "aob growth requires nutrient 'nh4' defined in Nutrients section"); if (ino2 == 0) error->all(FLERR, "aob growth requires nutrient 'no2' defined in Nutrients section"); if (io2 == 0) error->all(FLERR, "aob growth requires nutrient 'o2' defined in Nutrients section"); } else if (strcmp(name, "nob") == 0 || strcmp(name, "NOB") == 0) { species[i] = NOB; if (ino2 == 0) error->all(FLERR, "nob growth