/****************************************************************************** * * Copyright (C) 2019 by * The Salk Institute for Biological Studies and * Pittsburgh Supercomputing Center, Carnegie Mellon University * * Use of this source code is governed by an MIT-style * license that can be found in the LICENSE file or at * https://opensource.org/licenses/MIT. * ******************************************************************************/ #ifndef SRC4_MOLECULE_H_ #define SRC4_MOLECULE_H_ #include "bng/bng.h" #include "defines.h" namespace BNG { class SpeciesContainer; } namespace MCell { class Partition; // WARNING: do not change these values, checkpointed models use them // TODO: probably print this out in some reasonable form into checkpoints, it is printed as a number now enum molecule_flag_t { // volume/surface information is only cached from BNG CplxInstance MOLECULE_FLAG_SURF = 1 << 0, // originally TYPE_SURF MOLECULE_FLAG_VOL = 1 << 1, // originally TYPE_VOL MOLECULE_FLAG_MATURE = 1 << 2, // originally MATURE_MOLECULE MOLECULE_FLAG_ACT_CLAMPED = 1 << 3, // originally ACT_CLAMPED MOLECULE_FLAG_SCHEDULE_UNIMOL_RXN = 1 << 4, MOLECULE_FLAG_RESCHEDULE_UNIMOL_RXN_ON_NEXT_RXN_RATE_UPDATE = 1 << 5, // flags needed for concentration clamp handling, // only one of them may be set MOLECULE_FLAG_CLAMP_ORIENTATION_UP = 1 << 6, MOLECULE_FLAG_CLAMP_ORIENTATION_DOWN = 1 << 7, MOLECULE_FLAG_NO_NEED_TO_SCHEDULE = 1 << 14, MOLECULE_FLAG_DEFUNCT = 1 << 15, }; /** * Base class for all molecules. */ class Molecule { public: // Warning: ctors do not reset surf or vol data Molecule() : id(MOLECULE_ID_INVALID), species_id(SPECIES_ID_INVALID), flags(0), diffusion_time(TIME_INVALID), unimol_rxn_time(TIME_FOREVER), birthday(TIME_INVALID) { } Molecule(const Molecule& m) { *this = m; } // volume molecule Molecule( const molecule_id_t id_, const species_id_t species_id_, const Vec3& pos_, const double birthday_ ) : id(id_), species_id(species_id_), flags(MOLECULE_FLAG_VOL), diffusion_time(TIME_INVALID), unimol_rxn_time(TIME_INVALID), birthday(birthday_) { v.pos = pos_; v.subpart_index = SUBPART_INDEX_INVALID; v.reactant_subpart_index = SUBPART_INDEX_INVALID; v.counted_volume_index = COUNTED_VOLUME_INDEX_INVALID; v.previous_wall_index = WALL_INDEX_INVALID; } // surface molecule Molecule( const molecule_id_t id_, const species_id_t species_id_, const Vec2& pos2d, const double birthday_ ) : id(id_), species_id(species_id_), flags(MOLECULE_FLAG_SURF), diffusion_time(TIME_INVALID), unimol_rxn_time(TIME_INVALID), birthday(birthday_) { s.pos = pos2d; s.orientation = ORIENTATION_NONE; s.wall_index = WALL_INDEX_INVALID; s.grid_tile_index = TILE_INDEX_INVALID; } // assuming that this function has no virtual methods and has only POD types // gcc9 reports a warning but, the memcpy here is safe #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wclass-memaccess" void operator = (const Molecule& m) { memcpy(this, &m, sizeof(Molecule)); } #pragma GCC diagnostic pop void reset_vol_data() { v.pos = Vec3(POS_INVALID); v.subpart_index = SUBPART_INDEX_INVALID; v.reactant_subpart_index = SUBPART_INDEX_INVALID; v.counted_volume_index = COUNTED_VOLUME_INDEX_INVALID; v.previous_wall_index = WALL_INDEX_INVALID; } void reset_surf_data() { s.pos = Vec2(POS_INVALID); s.orientation = ORIENTATION_NONE; s.wall_index = WALL_INDEX_INVALID; s.grid_tile_index = TILE_INDEX_INVALID; } // may set flag for optimizations // called when a molecule is added to partition void set_no_need_to_schedule_flag(const BNG::SpeciesContainer& all_species); // data is ordered to avoid alignment holes (for 64-bit floats) molecule_id_t id; // unique molecule id (for now it is unique per partition but should be world-wide unique) species_id_t species_id; uint flags; // time for which it was scheduled, based on this value Partition creates 'ready list' // for DiffuseAndReactEvent double diffusion_time; // time assigned for unimol rxn, TIME_INVALID if time has not been set or molecule has no unimol rxn, // TIME_FOREVER if the probability of an existing unimol rxn is 0 double unimol_rxn_time; // - time when the molecule was released or created // - used when determining whether this molecule is mature // - release delay time is not added to the birthday time, a newly released molecule // with release delay will have its birthday at the beginning of the iteration double birthday; // update assignment operator when modifying this union { // volume molecule data struct { Vec3 pos; subpart_index_t subpart_index; // during diffusion the molecules' subpart index might change but the reactant_subpart_index // stays the same until its moved in the Partition's volume_molecule_reactants_per_subpart[] array subpart_index_t reactant_subpart_index; // do not assign directly, use set_counted_volume_and_compartment istead counted_volume_index_t counted_volume_index; // needed for clamp concentration handling wall_index_t previous_wall_index; } v; // surface molecule data struct { Vec2 pos; // we probably do not want subpart index, wall index serves this purpose //subpart_index_t subpart_index; orientation_t orientation; wall_index_t wall_index; tile_index_t grid_tile_index; } s; }; bool has_flag(uint flag) const { assert(__builtin_popcount(flag) == 1); return (flags & flag) != 0; } void set_flag(uint flag) { assert(__builtin_popcount(flag) == 1); flags = flags | flag; } void clear_flag(uint flag) { assert(__builtin_popcount(flag) == 1); flags = flags & ~flag; } void set_clamp_orientation(orientation_t value) { assert(value >= ORIENTATION_DOWN && value <= ORIENTATION_UP); if (value == ORIENTATION_NONE) { clear_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_UP); clear_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_DOWN); } else if (value == ORIENTATION_DOWN) { clear_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_UP); set_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_DOWN); } else { set_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_UP); clear_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_DOWN); } } int get_clamp_orientation() const { if (has_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_UP)) { return ORIENTATION_UP; } else if (has_flag(MOLECULE_FLAG_CLAMP_ORIENTATION_DOWN)) { return ORIENTATION_DOWN; } else { return ORIENTATION_NONE; } } bool is_vol() const { bool res = has_flag(MOLECULE_FLAG_VOL); if (res) { assert(!is_surf()); } return res; } bool is_surf() const { bool res = has_flag(MOLECULE_FLAG_SURF); if (res) { assert(!is_vol()); } return res; } bool is_defunct() const { // TODO LATER: molecules with release_delay > 0 were actually not released yet, // but is seems that mcell3 does not care, so let's keep it like this for now return has_flag(MOLECULE_FLAG_DEFUNCT) != 0; } void set_is_defunct() { assert(!is_defunct() && "We really should not be defuncting one molecule multiple times"); flags |= MOLECULE_FLAG_DEFUNCT; } orientation_t get_orientation() const { if (is_surf()) { return s.orientation; } else { return ORIENTATION_NONE; // or not set? } } void dump(const std::string ind="") const; void dump( const Partition& p, const std::string extra_comment, const std::string ind = " ", const uint64_t iteration = 0, const double time = 0, const bool print_position = true, const bool print_flags = false ) const; std::string to_string() const; static void dump_array(const std::vector& vec); }; } // namespace mcell #endif // SRC4_MOLECULE_H_