diff --git a/src/KOKKOS/neigh_bond_kokkos.cpp b/src/KOKKOS/neigh_bond_kokkos.cpp
index 2fdb09f880..8f054e46b5 100644
--- a/src/KOKKOS/neigh_bond_kokkos.cpp
+++ b/src/KOKKOS/neigh_bond_kokkos.cpp
@@ -75,30 +75,27 @@ void NeighBondKokkos<DeviceType>::init_topology_kk() {
 
   // 1st time allocation of topology lists
 
-  if (atom->molecular && atom->nbonds && maxbond == 0) {
-    if (nprocs == 1) maxbond = atom->nbonds;
-    else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / nprocs);
-    memoryKK->create_kokkos(k_bondlist,neighbor->bondlist,maxbond,3,"neigh:neighbor->bondlist");
-  }
-
-  if (atom->molecular && atom->nangles && maxangle == 0) {
-    if (nprocs == 1) maxangle = atom->nangles;
-    else maxangle = static_cast<int> (LB_FACTOR * atom->nangles / nprocs);
-    memoryKK->create_kokkos(k_anglelist,neighbor->anglelist,maxangle,4,"neigh:neighbor->anglelist");
-  }
-
-  if (atom->molecular && atom->ndihedrals && maxdihedral == 0) {
-    if (nprocs == 1) maxdihedral = atom->ndihedrals;
-    else maxdihedral = static_cast<int>
-           (LB_FACTOR * atom->ndihedrals / nprocs);
-    memoryKK->create_kokkos(k_dihedrallist,neighbor->dihedrallist,maxdihedral,5,"neigh:neighbor->dihedrallist");
-  }
-
-  if (atom->molecular && atom->nimpropers && maximproper == 0) {
-    if (nprocs == 1) maximproper = atom->nimpropers;
-    else maximproper = static_cast<int>
-           (LB_FACTOR * atom->nimpropers / nprocs);
-    memoryKK->create_kokkos(k_improperlist,neighbor->improperlist,maximproper,5,"neigh:neighbor->improperlist");
+  if (atom->molecular != Atom::ATOMIC) {
+    if (atom->nbonds && maxbond == 0) {
+      if (nprocs == 1) maxbond = atom->nbonds;
+      else maxbond = static_cast<int> (LB_FACTOR * atom->nbonds / nprocs);
+      memoryKK->create_kokkos(k_bondlist,neighbor->bondlist,maxbond,3,"neigh:neighbor->bondlist");
+    }
+    if (atom->nangles && maxangle == 0) {
+      if (nprocs == 1) maxangle = atom->nangles;
+      else maxangle = static_cast<int> (LB_FACTOR * atom->nangles / nprocs);
+      memoryKK->create_kokkos(k_anglelist,neighbor->anglelist,maxangle,4,"neigh:neighbor->anglelist");
+    }
+    if (atom->ndihedrals && maxdihedral == 0) {
+      if (nprocs == 1) maxdihedral = atom->ndihedrals;
+      else maxdihedral = static_cast<int> (LB_FACTOR * atom->ndihedrals / nprocs);
+      memoryKK->create_kokkos(k_dihedrallist,neighbor->dihedrallist,maxdihedral,5,"neigh:neighbor->dihedrallist");
+    }
+    if (atom->nimpropers && maximproper == 0) {
+      if (nprocs == 1) maximproper = atom->nimpropers;
+      else maximproper = static_cast<int> (LB_FACTOR * atom->nimpropers / nprocs);
+      memoryKK->create_kokkos(k_improperlist,neighbor->improperlist,maximproper,5,"neigh:neighbor->improperlist");
+    }
   }
 
   // set flags that determine which topology neighboring routines to use
diff --git a/src/KOKKOS/neighbor_kokkos.cpp b/src/KOKKOS/neighbor_kokkos.cpp
index cadfccfa6b..1795d87611 100644
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@@ -316,7 +316,7 @@ void NeighborKokkos::build_kokkos(int topoflag)
 
   // build topology lists for bonds/angles/etc
 
-  if (atom->molecular && topoflag) build_topology();
+  if ((atom->molecular != Atom::ATOMIC) && topoflag) build_topology();
 }
 
 template<class DeviceType>
diff --git a/src/REPLICA/fix_hyper_global.cpp b/src/REPLICA/fix_hyper_global.cpp
index 9a396d2987..811346d959 100644
--- a/src/REPLICA/fix_hyper_global.cpp
+++ b/src/REPLICA/fix_hyper_global.cpp
@@ -125,7 +125,7 @@ void FixHyperGlobal::init()
   if (force->newton_pair == 0)
     error->all(FLERR,"Hyper global requires newton pair on");
 
-  if (atom->molecular && me == 0)
+  if ((atom->molecular != Atom::ATOMIC) && (me == 0))
     error->warning(FLERR,"Hyper global for molecular systems "
                    "requires care in defining hyperdynamic bonds");
 
diff --git a/src/REPLICA/fix_hyper_local.cpp b/src/REPLICA/fix_hyper_local.cpp
index c26a9a6735..3dec6bdc10 100644
--- a/src/REPLICA/fix_hyper_local.cpp
+++ b/src/REPLICA/fix_hyper_local.cpp
@@ -274,7 +274,7 @@ void FixHyperLocal::init()
   if (force->newton_pair == 0)
     error->all(FLERR,"Hyper local requires newton pair on");
 
-  if (atom->molecular && me == 0)
+  if ((atom->molecular != Atom::ATOMIC) && (me == 0))
     error->warning(FLERR,"Hyper local for molecular systems "
                    "requires care in defining hyperdynamic bonds");
 
diff --git a/src/USER-BOCS/compute_pressure_bocs.cpp b/src/USER-BOCS/compute_pressure_bocs.cpp
index 46aabf7903..b91a790d15 100644
--- a/src/USER-BOCS/compute_pressure_bocs.cpp
+++ b/src/USER-BOCS/compute_pressure_bocs.cpp
@@ -152,10 +152,12 @@ void ComputePressureBocs::init()
   vptr = nullptr;
 
   if (pairflag && force->pair) nvirial++;
-  if (bondflag && atom->molecular && force->bond) nvirial++;
-  if (angleflag && atom->molecular && force->angle) nvirial++;
-  if (dihedralflag && atom->molecular && force->dihedral) nvirial++;
-  if (improperflag && atom->molecular && force->improper) nvirial++;
+  if (atom->molecular != Atom::ATOMIC) {
+    if (bondflag && force->bond) nvirial++;
+    if (angleflag && force->angle) nvirial++;
+    if (dihedralflag && force->dihedral) nvirial++;
+    if (improperflag && force->improper) nvirial++;
+  }
   if (fixflag)
     for (int i = 0; i < modify->nfix; i++)
       if (modify->fix[i]->virial_flag) nvirial++;
diff --git a/src/USER-INTEL/fix_intel.cpp b/src/USER-INTEL/fix_intel.cpp
index 21cc106e31..92b47179dd 100644
--- a/src/USER-INTEL/fix_intel.cpp
+++ b/src/USER-INTEL/fix_intel.cpp
@@ -441,7 +441,7 @@ void FixIntel::pair_init_check(const bool cdmessage)
   }
 
   #ifndef LMP_INTEL_NBOR_COMPAT
-  if (force->pair->manybody_flag && atom->molecular) {
+  if (force->pair->manybody_flag && (atom->molecular != Atom::ATOMIC)) {
     int flag = 0;
     if (atom->nbonds > 0 && force->special_lj[1] == 0.0 &&
         force->special_coul[1] == 0.0) flag = 1;
@@ -456,7 +456,7 @@ void FixIntel::pair_init_check(const bool cdmessage)
   #endif
 
   int need_tag = 0;
-  if (atom->molecular) need_tag = 1;
+  if (atom->molecular != Atom::ATOMIC) need_tag = 1;
 
   // Clear buffers used for pair style
   char kmode[80];
@@ -507,8 +507,8 @@ void FixIntel::pair_init_check(const bool cdmessage)
 
 void FixIntel::bond_init_check()
 {
-  if (_offload_balance != 0.0 && atom->molecular &&
-      force->newton_pair != force->newton_bond)
+  if ((_offload_balance != 0.0) && (atom->molecular != Atom::ATOMIC)
+      && (force->newton_pair != force->newton_bond))
     error->all(FLERR,
       "USER-INTEL package requires same setting for newton bond and non-bond.");
 
diff --git a/src/USER-INTEL/npair_full_bin_ghost_intel.cpp b/src/USER-INTEL/npair_full_bin_ghost_intel.cpp
index ae7ac78bc8..50a3586d31 100644
--- a/src/USER-INTEL/npair_full_bin_ghost_intel.cpp
+++ b/src/USER-INTEL/npair_full_bin_ghost_intel.cpp
@@ -88,7 +88,7 @@ void NPairFullBinGhostIntel::fbi(NeighList * list,
                      _fix->nbor_pack_width());
 
   int need_ic = 0;
-  if (atom->molecular)
+  if (atom->molecular != Atom::ATOMIC)
     dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax,
                          neighbor->cutneighmax);
 
@@ -420,7 +420,7 @@ void NPairFullBinGhostIntel::fbi(const int offload, NeighList * list,
           }
         } // for u
 
-        if (molecular && i < nlocal) {
+        if ((molecular != Atom::ATOMIC) && (i < nlocal)) {
           int alln = n;
           n = 0;
           #if defined(LMP_SIMD_COMPILER)
diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp
index ae1bb4747d..f626718596 100644
--- a/src/compute_pressure.cpp
+++ b/src/compute_pressure.cpp
@@ -189,10 +189,12 @@ void ComputePressure::init()
 
   if (pairhybridflag && force->pair) nvirial++;
   if (pairflag && force->pair) nvirial++;
-  if (bondflag && atom->molecular && force->bond) nvirial++;
-  if (angleflag && atom->molecular && force->angle) nvirial++;
-  if (dihedralflag && atom->molecular && force->dihedral) nvirial++;
-  if (improperflag && atom->molecular && force->improper) nvirial++;
+  if (atom->molecular != Atom::ATOMIC) {
+    if (bondflag && force->bond) nvirial++;
+    if (angleflag && force->angle) nvirial++;
+    if (dihedralflag && force->dihedral) nvirial++;
+    if (improperflag && force->improper) nvirial++;
+  }
   if (fixflag)
     for (int i = 0; i < modify->nfix; i++)
       if (modify->fix[i]->thermo_virial) nvirial++;
diff --git a/src/domain.cpp b/src/domain.cpp
index b725495828..b5d6675a3f 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -727,7 +727,7 @@ void Domain::image_check()
   // if running verlet/split, don't check on KSpace partition since
   //    it has no ghost atoms and thus bond partners won't exist
 
-  if (!atom->molecular) return;
+  if (atom->molecular == Atom::ATOMIC) return;
   if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
   if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
       universe->iworld != 0) return;
@@ -837,7 +837,7 @@ void Domain::box_too_small_check()
   // if running verlet/split, don't check on KSpace partition since
   //    it has no ghost atoms and thus bond partners won't exist
 
-  if (!atom->molecular) return;
+  if (atom->molecular == Atom::ATOMIC) return;
   if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
   if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
       universe->iworld != 0) return;
diff --git a/src/finish.cpp b/src/finish.cpp
index 3f057a00b8..7525f5e632 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -332,7 +332,7 @@ void Finish::end(int flag)
     mpi_timings("Pair",timer,Timer::PAIR, world,nprocs,
                 nthreads,me,time_loop,screen,logfile);
 
-    if (atom->molecular)
+    if (atom->molecular != Atom::ATOMIC)
       mpi_timings("Bond",timer,Timer::BOND,world,nprocs,
                   nthreads,me,time_loop,screen,logfile);
 
@@ -391,7 +391,7 @@ void Finish::end(int flag)
       utils::logmesg(lmp,fmt::format(thr_fmt,me,thr_total,thr_total/time_loop*100.0));
 
       omp_times(fixomp,"Pair",Timer::PAIR,nthreads,screen,logfile);
-      if (atom->molecular)
+      if (atom->molecular != Atom::ATOMIC)
         omp_times(fixomp,"Bond",Timer::BOND,nthreads,screen,logfile);
       if (force->kspace)
         omp_times(fixomp,"Kspace",Timer::KSPACE,nthreads,screen,logfile);
@@ -585,7 +585,7 @@ void Finish::end(int flag)
       mesg += fmt::format("Total # of neighbors = {:.8g}\n",nall);
       if (atom->natoms > 0)
         mesg += fmt::format("Ave neighs/atom = {:.8}\n",nall/atom->natoms);
-      if (atom->molecular && atom->natoms > 0)
+      if ((atom->molecular != Atom::ATOMIC) && (atom->natoms > 0))
         mesg += fmt::format("Ave special neighs/atom = {:.8}\n",
                             nspec_all/atom->natoms);
       mesg += fmt::format("Neighbor list builds = {}\n",neighbor->ncalls);
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 4eabdcc7b4..58737b815d 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1296,7 +1296,7 @@ void Neighbor::init_topology()
 {
   int i,m;
 
-  if (!atom->molecular) return;
+  if (atom->molecular == Atom::ATOMIC) return;
 
   // set flags that determine which topology neighbor classes to use
   // these settings could change from run to run, depending on fixes defined
@@ -1801,7 +1801,7 @@ int Neighbor::choose_pair(NeighRequest *rq)
 
     if (molecular != Atom::ATOMIC) {
       if (mask & NP_ATOMONLY) continue;
-    } else if (!molecular) {
+    } else if (molecular == Atom::ATOMIC) {
       if (mask & NP_MOLONLY) continue;
     }
 
@@ -2097,7 +2097,7 @@ void Neighbor::build(int topoflag)
 
   // build topology lists for bonds/angles/etc
 
-  if (atom->molecular && topoflag) build_topology();
+  if ((atom->molecular != Atom::ATOMIC) && topoflag) build_topology();
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/reset_mol_ids.cpp b/src/reset_mol_ids.cpp
index 10337fd18f..970ac6beb4 100644
--- a/src/reset_mol_ids.cpp
+++ b/src/reset_mol_ids.cpp
@@ -64,7 +64,7 @@ void ResetMolIDs::command(int narg, char **arg)
     error->all(FLERR,"Reset_mol_ids command before simulation box is defined");
   if (atom->tag_enable == 0)
     error->all(FLERR,"Cannot use reset_mol_ids unless atoms have IDs");
-  if (atom->molecular != 1)
+  if (atom->molecular != Atom::MOLECULAR)
     error->all(FLERR,"Can only use reset_mol_ids on molecular systems");
 
   // process args