diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp
index 90c8b07b02..921d469e31 100644
--- a/unittest/formats/test_atom_styles.cpp
+++ b/unittest/formats/test_atom_styles.cpp
@@ -2743,7 +2743,6 @@ TEST_F(AtomStyleTest, template)
     expected.has_x         = true;
     expected.has_v         = true;
     expected.has_f         = true;
-    expected.q_flag        = 1;
     expected.molecule_flag = 1;
     expected.molindex_flag = 1;
     expected.molatom_flag  = 1;
@@ -3014,7 +3013,7 @@ TEST_F(AtomStyleTest, template)
     ASSERT_EQ(lmp->atom->nghost, 0);
     ASSERT_NE(lmp->atom->nmax, -1);
     ASSERT_EQ(lmp->atom->tag_enable, 1);
-    ASSERT_EQ(lmp->atom->q_flag, 1);
+    ASSERT_EQ(lmp->atom->q_flag, 0);
     ASSERT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
     ASSERT_EQ(lmp->atom->ntypes, 4);
     ASSERT_EQ(lmp->atom->tag_consecutive(), 0);
@@ -3118,6 +3117,434 @@ TEST_F(AtomStyleTest, template)
     ASSERT_EQ(molatom[GETIDX(16)], -1);
 }
 
+TEST_F(AtomStyleTest, template_charge)
+{
+    if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();
+    BEGIN_HIDE_OUTPUT();
+    command("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    command("atom_style hybrid template twomols charge");
+    command("newton on");
+    END_HIDE_OUTPUT();
+
+    AtomState expected;
+    expected.atom_style    = "hybrid";
+    expected.molecular     = Atom::TEMPLATE;
+    expected.nbondtypes    = 2;
+    expected.nangletypes   = 2;
+    expected.tag_enable    = 1;
+    expected.has_type      = true;
+    expected.has_mask      = true;
+    expected.has_image     = true;
+    expected.has_x         = true;
+    expected.has_v         = true;
+    expected.has_f         = true;
+    expected.molecule_flag = 1;
+    expected.molindex_flag = 1;
+    expected.molatom_flag  = 1;
+    expected.q_flag        = 1;
+    expected.nmolecule     = 2;
+    expected.map_style     = 3;
+
+    ASSERT_ATOM_STATE_EQ(lmp->atom, expected);
+
+    auto *hybrid = dynamic_cast<AtomVecHybrid *>(lmp->atom->avec);
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_EQ(hybrid->nstyles, 2);
+    ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("template"));
+    ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("charge"));
+    ASSERT_NE(hybrid->styles[0], nullptr);
+    ASSERT_NE(hybrid->styles[1], nullptr);
+
+    BEGIN_HIDE_OUTPUT();
+    command("create_box 4 box bond/types 2 angle/types 2 ");
+    command("create_atoms 0 single -2.0  2.0  0.1 mol twomols 65234");
+    command("create_atoms 0 single -2.0 -2.0 -0.1 mol twomols 62346");
+    command("create_atoms 0 single  2.0  2.0 -0.1 mol twomols 61354");
+    command("create_atoms 3 single  2.0 -2.0  0.1");
+    command("create_atoms 3 single  2.0  2.0 -2.1");
+    command("create_atoms 4 single  2.0 -2.0  2.1");
+    command("mass 1 16.0");
+    command("mass 2 1.0");
+    command("mass 3 12.0");
+    command("mass 4 16.0");
+    command("set atom 10 charge 0.7");
+    command("set atom 11 charge -0.35");
+    command("set atom 12 charge -0.35");
+    command("bond_style zero");
+    command("bond_coeff 1 1.0");
+    command("bond_coeff 2 1.16");
+    command("angle_style zero");
+    command("angle_coeff * 109.0");
+    command("pair_coeff * *");
+    END_HIDE_OUTPUT();
+    ASSERT_NE(lmp->atom->avec, nullptr);
+    hybrid = dynamic_cast<AtomVecHybrid *>(lmp->atom->avec);
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_EQ(hybrid->nstyles, 2);
+    ASSERT_THAT(std::string(hybrid->keywords[0]), Eq("template"));
+    ASSERT_THAT(std::string(hybrid->keywords[1]), Eq("charge"));
+    ASSERT_NE(hybrid->styles[0], nullptr);
+    ASSERT_NE(hybrid->styles[1], nullptr);
+
+    ASSERT_EQ(lmp->atom->natoms, 12);
+    ASSERT_EQ(lmp->atom->nbonds, 6);
+    ASSERT_EQ(lmp->atom->nbondtypes, 2);
+    ASSERT_EQ(lmp->atom->nangles, 3);
+    ASSERT_EQ(lmp->atom->nangletypes, 2);
+    ASSERT_EQ(lmp->atom->nellipsoids, 0);
+    ASSERT_EQ(lmp->atom->nlocal, 12);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    ASSERT_EQ(lmp->atom->q_flag, 1);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_EQ(lmp->atom->nextra_grow, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart, 0);
+    ASSERT_EQ(lmp->atom->nextra_border, 0);
+    ASSERT_EQ(lmp->atom->nextra_grow_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_restart_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_border_max, 0);
+    ASSERT_EQ(lmp->atom->nextra_store, 0);
+    ASSERT_EQ(lmp->atom->extra_grow, nullptr);
+    ASSERT_EQ(lmp->atom->extra_restart, nullptr);
+    ASSERT_EQ(lmp->atom->extra_border, nullptr);
+    ASSERT_EQ(lmp->atom->extra, nullptr);
+
+    ASSERT_NE(lmp->atom->mass, nullptr);
+    ASSERT_NE(lmp->atom->mass_setflag, nullptr);
+
+    BEGIN_HIDE_OUTPUT();
+    command("write_data test_atom_styles.data");
+    command("clear");
+    command("units real");
+    command("newton off on");
+    command("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    command("atom_style hybrid template twomols charge");
+    command("pair_style zero 4.0");
+    command("bond_style zero");
+    command("angle_style zero");
+    command("atom_modify map array");
+    command("read_data test_atom_styles.data");
+    END_HIDE_OUTPUT();
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_NE(lmp->atom->avec, nullptr);
+
+    ASSERT_EQ(lmp->atom->natoms, 12);
+    ASSERT_EQ(lmp->atom->nlocal, 12);
+    ASSERT_EQ(lmp->atom->nbonds, 6);
+    ASSERT_EQ(lmp->atom->nangles, 3);
+    ASSERT_EQ(lmp->atom->nbondtypes, 2);
+    ASSERT_EQ(lmp->atom->nangletypes, 2);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_NE(lmp->atom->sametag, nullptr);
+    ASSERT_EQ(lmp->atom->tag_consecutive(), 1);
+    ASSERT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    ASSERT_EQ(lmp->atom->map_user, 1);
+    ASSERT_EQ(lmp->atom->map_tag_max, 12);
+
+    auto *molecule = lmp->atom->molecule;
+    auto *molindex = lmp->atom->molindex;
+    auto *molatom  = lmp->atom->molatom;
+
+    ASSERT_EQ(molecule[GETIDX(1)], 1);
+    ASSERT_EQ(molecule[GETIDX(2)], 1);
+    ASSERT_EQ(molecule[GETIDX(3)], 1);
+    ASSERT_EQ(molecule[GETIDX(4)], 2);
+    ASSERT_EQ(molecule[GETIDX(5)], 2);
+    ASSERT_EQ(molecule[GETIDX(6)], 2);
+    ASSERT_EQ(molecule[GETIDX(7)], 3);
+    ASSERT_EQ(molecule[GETIDX(8)], 3);
+    ASSERT_EQ(molecule[GETIDX(9)], 3);
+    ASSERT_EQ(molecule[GETIDX(10)], 0);
+    ASSERT_EQ(molecule[GETIDX(11)], 0);
+    ASSERT_EQ(molecule[GETIDX(12)], 0);
+    ASSERT_EQ(molindex[GETIDX(1)], 0);
+    ASSERT_EQ(molindex[GETIDX(2)], 0);
+    ASSERT_EQ(molindex[GETIDX(3)], 0);
+    ASSERT_EQ(molindex[GETIDX(4)], 0);
+    ASSERT_EQ(molindex[GETIDX(5)], 0);
+    ASSERT_EQ(molindex[GETIDX(6)], 0);
+    ASSERT_EQ(molindex[GETIDX(7)], 0);
+    ASSERT_EQ(molindex[GETIDX(8)], 0);
+    ASSERT_EQ(molindex[GETIDX(9)], 0);
+    ASSERT_EQ(molindex[GETIDX(10)], -1);
+    ASSERT_EQ(molindex[GETIDX(11)], -1);
+    ASSERT_EQ(molindex[GETIDX(12)], -1);
+    ASSERT_EQ(molatom[GETIDX(1)], 0);
+    ASSERT_EQ(molatom[GETIDX(2)], 1);
+    ASSERT_EQ(molatom[GETIDX(3)], 2);
+    ASSERT_EQ(molatom[GETIDX(4)], 0);
+    ASSERT_EQ(molatom[GETIDX(5)], 1);
+    ASSERT_EQ(molatom[GETIDX(6)], 2);
+    ASSERT_EQ(molatom[GETIDX(7)], 0);
+    ASSERT_EQ(molatom[GETIDX(8)], 1);
+    ASSERT_EQ(molatom[GETIDX(9)], 2);
+    ASSERT_EQ(molatom[GETIDX(10)], -1);
+    ASSERT_EQ(molatom[GETIDX(11)], -1);
+    ASSERT_EQ(molatom[GETIDX(12)], -1);
+
+    BEGIN_HIDE_OUTPUT();
+    command("clear");
+    command("units real");
+    command("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    command("atom_style hybrid template twomols charge");
+    command("pair_style zero 4.0");
+    command("bond_style zero");
+    command("angle_style zero");
+    command("atom_modify map array");
+    command("read_data test_atom_styles.data");
+    END_HIDE_OUTPUT();
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_NE(lmp->atom->avec, nullptr);
+
+    ASSERT_EQ(lmp->atom->natoms, 12);
+    ASSERT_EQ(lmp->atom->nlocal, 12);
+    ASSERT_EQ(lmp->atom->nbonds, 6);
+    ASSERT_EQ(lmp->atom->nangles, 3);
+    ASSERT_EQ(lmp->atom->nbondtypes, 2);
+    ASSERT_EQ(lmp->atom->nangletypes, 2);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_NE(lmp->atom->sametag, nullptr);
+    ASSERT_EQ(lmp->atom->tag_consecutive(), 1);
+    ASSERT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    ASSERT_EQ(lmp->atom->map_user, 1);
+    ASSERT_EQ(lmp->atom->map_tag_max, 12);
+
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+
+    ASSERT_EQ(molindex[GETIDX(1)], 0);
+    ASSERT_EQ(molindex[GETIDX(2)], 0);
+    ASSERT_EQ(molindex[GETIDX(3)], 0);
+    ASSERT_EQ(molindex[GETIDX(4)], 0);
+    ASSERT_EQ(molindex[GETIDX(5)], 0);
+    ASSERT_EQ(molindex[GETIDX(6)], 0);
+    ASSERT_EQ(molindex[GETIDX(7)], 0);
+    ASSERT_EQ(molindex[GETIDX(8)], 0);
+    ASSERT_EQ(molindex[GETIDX(9)], 0);
+    ASSERT_EQ(molindex[GETIDX(10)], -1);
+    ASSERT_EQ(molindex[GETIDX(11)], -1);
+    ASSERT_EQ(molindex[GETIDX(12)], -1);
+    ASSERT_EQ(molatom[GETIDX(1)], 0);
+    ASSERT_EQ(molatom[GETIDX(2)], 1);
+    ASSERT_EQ(molatom[GETIDX(3)], 2);
+    ASSERT_EQ(molatom[GETIDX(4)], 0);
+    ASSERT_EQ(molatom[GETIDX(5)], 1);
+    ASSERT_EQ(molatom[GETIDX(6)], 2);
+    ASSERT_EQ(molatom[GETIDX(7)], 0);
+    ASSERT_EQ(molatom[GETIDX(8)], 1);
+    ASSERT_EQ(molatom[GETIDX(9)], 2);
+    ASSERT_EQ(molatom[GETIDX(10)], -1);
+    ASSERT_EQ(molatom[GETIDX(11)], -1);
+    ASSERT_EQ(molatom[GETIDX(12)], -1);
+
+    auto *x    = lmp->atom->x;
+    auto *v    = lmp->atom->v;
+    auto *type = lmp->atom->type;
+    auto *q    = lmp->atom->q;
+
+    EXPECT_NEAR(x[GETIDX(10)][0], 2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(10)][1], -2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(10)][2], 0.1, EPSILON);
+    EXPECT_NEAR(x[GETIDX(11)][0], 2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(11)][1], 2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(11)][2], -2.1, EPSILON);
+    EXPECT_NEAR(x[GETIDX(12)][0], 2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(12)][1], -2.0, EPSILON);
+    EXPECT_NEAR(x[GETIDX(12)][2], 2.1, EPSILON);
+    EXPECT_NEAR(v[GETIDX(1)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(1)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(1)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(2)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(3)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(4)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(5)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(6)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(7)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(8)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(9)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(10)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(11)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(11)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(11)][2], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(12)][0], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(12)][1], 0.0, EPSILON);
+    EXPECT_NEAR(v[GETIDX(12)][2], 0.0, EPSILON);
+    EXPECT_NEAR(q[GETIDX(1)], -0.8472, EPSILON);
+    EXPECT_NEAR(q[GETIDX(2)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(3)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(4)], -0.8472, EPSILON);
+    EXPECT_NEAR(q[GETIDX(5)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(6)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(7)], -0.8472, EPSILON);
+    EXPECT_NEAR(q[GETIDX(8)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(9)], 0.4236, EPSILON);
+    EXPECT_NEAR(q[GETIDX(10)], 0.7, EPSILON);
+    EXPECT_NEAR(q[GETIDX(11)], -0.35, EPSILON);
+    EXPECT_NEAR(q[GETIDX(12)], -0.35, EPSILON);
+    ASSERT_EQ(type[GETIDX(1)], 1);
+    ASSERT_EQ(type[GETIDX(2)], 2);
+    ASSERT_EQ(type[GETIDX(3)], 2);
+    ASSERT_EQ(type[GETIDX(4)], 1);
+    ASSERT_EQ(type[GETIDX(5)], 2);
+    ASSERT_EQ(type[GETIDX(6)], 2);
+    ASSERT_EQ(type[GETIDX(7)], 1);
+    ASSERT_EQ(type[GETIDX(8)], 2);
+    ASSERT_EQ(type[GETIDX(9)], 2);
+    ASSERT_EQ(type[GETIDX(10)], 3);
+    ASSERT_EQ(type[GETIDX(11)], 3);
+    ASSERT_EQ(type[GETIDX(12)], 4);
+
+    BEGIN_HIDE_OUTPUT();
+    command("group two id 7:10");
+    command("delete_atoms group two compress no");
+    command("write_restart test_atom_styles.restart");
+    command("clear");
+    command("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    command("read_restart test_atom_styles.restart");
+    command("replicate 1 1 2 bbox");
+    END_HIDE_OUTPUT();
+    ASSERT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    ASSERT_NE(lmp->atom->avec, nullptr);
+    ASSERT_EQ(lmp->atom->natoms, 16);
+
+    ASSERT_EQ(lmp->atom->nlocal, 16);
+    ASSERT_EQ(lmp->atom->nghost, 0);
+    ASSERT_NE(lmp->atom->nmax, -1);
+    ASSERT_EQ(lmp->atom->tag_enable, 1);
+    ASSERT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    ASSERT_EQ(lmp->atom->ntypes, 4);
+    ASSERT_EQ(lmp->atom->nbonds, 8);
+    ASSERT_EQ(lmp->atom->nangles, 4);
+    ASSERT_EQ(lmp->atom->tag_consecutive(), 0);
+    ASSERT_EQ(lmp->atom->map_tag_max, 24);
+
+    type     = lmp->atom->type;
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+    ASSERT_EQ(type[GETIDX(1)], 1);
+    ASSERT_EQ(type[GETIDX(2)], 2);
+    ASSERT_EQ(type[GETIDX(3)], 2);
+    ASSERT_EQ(type[GETIDX(4)], 1);
+    ASSERT_EQ(type[GETIDX(5)], 2);
+    ASSERT_EQ(type[GETIDX(6)], 2);
+    ASSERT_EQ(type[GETIDX(11)], 3);
+    ASSERT_EQ(type[GETIDX(12)], 4);
+    ASSERT_EQ(type[GETIDX(13)], 1);
+    ASSERT_EQ(type[GETIDX(14)], 2);
+    ASSERT_EQ(type[GETIDX(15)], 2);
+    ASSERT_EQ(type[GETIDX(16)], 1);
+    ASSERT_EQ(type[GETIDX(17)], 2);
+    ASSERT_EQ(type[GETIDX(18)], 2);
+    ASSERT_EQ(type[GETIDX(23)], 3);
+    ASSERT_EQ(type[GETIDX(24)], 4);
+
+    ASSERT_EQ(molindex[GETIDX(1)], 0);
+    ASSERT_EQ(molindex[GETIDX(2)], 0);
+    ASSERT_EQ(molindex[GETIDX(3)], 0);
+    ASSERT_EQ(molindex[GETIDX(4)], 0);
+    ASSERT_EQ(molindex[GETIDX(5)], 0);
+    ASSERT_EQ(molindex[GETIDX(6)], 0);
+    ASSERT_EQ(molindex[GETIDX(11)], -1);
+    ASSERT_EQ(molindex[GETIDX(12)], -1);
+    ASSERT_EQ(molindex[GETIDX(13)], 0);
+    ASSERT_EQ(molindex[GETIDX(14)], 0);
+    ASSERT_EQ(molindex[GETIDX(15)], 0);
+    ASSERT_EQ(molindex[GETIDX(16)], 0);
+    ASSERT_EQ(molindex[GETIDX(17)], 0);
+    ASSERT_EQ(molindex[GETIDX(18)], 0);
+    ASSERT_EQ(molindex[GETIDX(23)], -1);
+    ASSERT_EQ(molindex[GETIDX(24)], -1);
+    ASSERT_EQ(molatom[GETIDX(1)], 0);
+    ASSERT_EQ(molatom[GETIDX(2)], 1);
+    ASSERT_EQ(molatom[GETIDX(3)], 2);
+    ASSERT_EQ(molatom[GETIDX(4)], 0);
+    ASSERT_EQ(molatom[GETIDX(5)], 1);
+    ASSERT_EQ(molatom[GETIDX(6)], 2);
+    ASSERT_EQ(molatom[GETIDX(11)], -1);
+    ASSERT_EQ(molatom[GETIDX(12)], -1);
+    ASSERT_EQ(molatom[GETIDX(13)], 0);
+    ASSERT_EQ(molatom[GETIDX(14)], 1);
+    ASSERT_EQ(molatom[GETIDX(15)], 2);
+    ASSERT_EQ(molatom[GETIDX(16)], 0);
+    ASSERT_EQ(molatom[GETIDX(17)], 1);
+    ASSERT_EQ(molatom[GETIDX(18)], 2);
+    ASSERT_EQ(molatom[GETIDX(23)], -1);
+    ASSERT_EQ(molatom[GETIDX(24)], -1);
+
+    BEGIN_HIDE_OUTPUT();
+    command("reset_atoms id");
+    END_HIDE_OUTPUT();
+    ASSERT_EQ(lmp->atom->tag_consecutive(), 1);
+    ASSERT_EQ(lmp->atom->map_tag_max, 16);
+
+    type     = lmp->atom->type;
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+    ASSERT_EQ(type[GETIDX(1)], 1);
+    ASSERT_EQ(type[GETIDX(2)], 2);
+    ASSERT_EQ(type[GETIDX(3)], 2);
+    ASSERT_EQ(type[GETIDX(4)], 1);
+    ASSERT_EQ(type[GETIDX(5)], 2);
+    ASSERT_EQ(type[GETIDX(6)], 2);
+    ASSERT_EQ(type[GETIDX(7)], 4);
+    ASSERT_EQ(type[GETIDX(8)], 3);
+    ASSERT_EQ(type[GETIDX(9)], 1);
+    ASSERT_EQ(type[GETIDX(10)], 2);
+    ASSERT_EQ(type[GETIDX(11)], 2);
+    ASSERT_EQ(type[GETIDX(12)], 1);
+    ASSERT_EQ(type[GETIDX(13)], 2);
+    ASSERT_EQ(type[GETIDX(14)], 2);
+    ASSERT_EQ(type[GETIDX(15)], 4);
+    ASSERT_EQ(type[GETIDX(16)], 3);
+    ASSERT_EQ(molatom[GETIDX(1)], 0);
+    ASSERT_EQ(molatom[GETIDX(2)], 1);
+    ASSERT_EQ(molatom[GETIDX(3)], 2);
+    ASSERT_EQ(molatom[GETIDX(4)], 0);
+    ASSERT_EQ(molatom[GETIDX(5)], 1);
+    ASSERT_EQ(molatom[GETIDX(6)], 2);
+    ASSERT_EQ(molatom[GETIDX(7)], -1);
+    ASSERT_EQ(molatom[GETIDX(8)], -1);
+    ASSERT_EQ(molatom[GETIDX(9)], 0);
+    ASSERT_EQ(molatom[GETIDX(10)], 1);
+    ASSERT_EQ(molatom[GETIDX(11)], 2);
+    ASSERT_EQ(molatom[GETIDX(12)], 0);
+    ASSERT_EQ(molatom[GETIDX(13)], 1);
+    ASSERT_EQ(molatom[GETIDX(14)], 2);
+    ASSERT_EQ(molatom[GETIDX(15)], -1);
+    ASSERT_EQ(molatom[GETIDX(16)], -1);
+}
+
 TEST_F(AtomStyleTest, bond)
 {
     if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();