diff --git a/unittest/formats/test_atom_styles.cpp b/unittest/formats/test_atom_styles.cpp
index 058d92567f..2dcefab6dd 100644
--- a/unittest/formats/test_atom_styles.cpp
+++ b/unittest/formats/test_atom_styles.cpp
@@ -51,10 +51,49 @@
 
 using LAMMPS_NS::utils::split_words;
 
+static void create_molecule_files()
+{
+    // create molecule files
+    const char h2o_file[] = "# Water molecule. SPC/E model.\n\n3 atoms\n2 bonds\n1 angles\n\n"
+                            "Coords\n\n1 1.12456 0.09298 1.27452\n"
+                            "2 1.53683 0.75606 1.89928\n3 0.49482 0.56390 0.65678\n\n"
+                            "Types\n\n1 1\n2 2\n3 2\n\n"
+                            "Charges\n\n1 -0.8472\n2 0.4236\n3 0.4236\n\n"
+                            "Bonds\n\n1 1 1 2\n2 1 1 3\n\n"
+                            "Angles\n\n1 1 2 1 3\n\n"
+                            "Shake Flags\n\n1 1\n2 1\n3 1\n\n"
+                            "Shake Atoms\n\n1 1 2 3\n2 1 2 3\n3 1 2 3\n\n"
+                            "Shake Bond Types\n\n1 1 1 1\n2 1 1 1\n3 1 1 1\n\n"
+                            "Special Bond Counts\n\n1 2 0 0\n2 1 1 0\n3 1 1 0\n\n"
+                            "Special Bonds\n\n1 2 3\n2 1 3\n3 1 2\n\n";
+    const char co2_file[] = "# CO2 molecule file. TraPPE model.\n\n"
+                            "3 atoms\n2 bonds\n1 angles\n\n"
+                            "Coords\n\n1 0.0 0.0 0.0\n2 -1.16 0.0 0.0\n3 1.16 0.0 0.0\n\n"
+                            "Types\n\n1 1\n2 2\n3 2\n\n"
+                            "Charges\n\n1 0.7\n2 -0.35\n3 -0.35\n\n"
+                            "Bonds\n\n1 1 1 2\n2 1 1 3\n\n"
+                            "Angles\n\n1 1 2 1 3\n\n"
+                            "Special Bond Counts\n\n1 2 0 0\n2 1 1 0\n3 1 1 0\n\n"
+                            "Special Bonds\n\n1 2 3\n2 1 3\n3 1 2\n\n";
+
+    FILE *fp = fopen("tmp.h2o.mol", "w");
+    if (fp) {
+        fputs(h2o_file, fp);
+        fclose(fp);
+    }
+    rename("tmp.h2o.mol", "h2o.mol");
+    fp = fopen("tmp.co2.mol", "w");
+    if (fp) {
+        fputs(co2_file, fp);
+        fclose(fp);
+    }
+    rename("tmp.co2.mol", "co2.mol");
+}
+
 // whether to print verbose output (i.e. not capturing LAMMPS screen output).
 bool verbose = false;
 
-const double EPSILON=5.0e-14;
+const double EPSILON = 5.0e-14;
 
 namespace LAMMPS_NS {
 using ::testing::Eq;
@@ -3432,6 +3471,1024 @@ TEST_F(AtomStyleTest, body_nparticle)
     ASSERT_EQ(bonus[3].ilocal, 6);
 }
 
+TEST_F(AtomStyleTest, template)
+{
+    if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();
+    create_molecule_files();
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style template twomols");
+    lmp->input->one("newton on");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("template"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+    EXPECT_EQ(lmp->atom->natoms, 0);
+    EXPECT_EQ(lmp->atom->nlocal, 0);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_EQ(lmp->atom->nmax, 1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->nellipsoids, 0);
+    EXPECT_EQ(lmp->atom->nlines, 0);
+    EXPECT_EQ(lmp->atom->ntris, 0);
+    EXPECT_EQ(lmp->atom->nbodies, 0);
+    EXPECT_EQ(lmp->atom->nbonds, 0);
+    EXPECT_EQ(lmp->atom->nangles, 0);
+    EXPECT_EQ(lmp->atom->ndihedrals, 0);
+    EXPECT_EQ(lmp->atom->nimpropers, 0);
+    EXPECT_EQ(lmp->atom->ntypes, 0);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->ndihedraltypes, 0);
+    EXPECT_EQ(lmp->atom->nimpropertypes, 0);
+    EXPECT_EQ(lmp->atom->bond_per_atom, 0);
+    EXPECT_EQ(lmp->atom->angle_per_atom, 0);
+    EXPECT_EQ(lmp->atom->dihedral_per_atom, 0);
+    EXPECT_EQ(lmp->atom->improper_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_bond_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_angle_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_dihedral_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_improper_per_atom, 0);
+
+    EXPECT_EQ(lmp->atom->sphere_flag, 0);
+    EXPECT_EQ(lmp->atom->ellipsoid_flag, 0);
+    EXPECT_EQ(lmp->atom->line_flag, 0);
+    EXPECT_EQ(lmp->atom->tri_flag, 0);
+    EXPECT_EQ(lmp->atom->body_flag, 0);
+    EXPECT_EQ(lmp->atom->peri_flag, 0);
+    EXPECT_EQ(lmp->atom->electron_flag, 0);
+    EXPECT_EQ(lmp->atom->wavepacket_flag, 0);
+    EXPECT_EQ(lmp->atom->sph_flag, 0);
+    EXPECT_EQ(lmp->atom->molecule_flag, 1);
+    EXPECT_EQ(lmp->atom->molindex_flag, 0);
+    EXPECT_EQ(lmp->atom->molatom_flag, 0);
+    EXPECT_EQ(lmp->atom->q_flag, 0);
+    EXPECT_EQ(lmp->atom->mu_flag, 0);
+    EXPECT_EQ(lmp->atom->rmass_flag, 0);
+    EXPECT_EQ(lmp->atom->radius_flag, 0);
+    EXPECT_EQ(lmp->atom->omega_flag, 0);
+    EXPECT_EQ(lmp->atom->torque_flag, 0);
+    EXPECT_EQ(lmp->atom->angmom_flag, 0);
+    EXPECT_EQ(lmp->atom->vfrac_flag, 0);
+    EXPECT_EQ(lmp->atom->spin_flag, 0);
+    EXPECT_EQ(lmp->atom->eradius_flag, 0);
+    EXPECT_EQ(lmp->atom->ervel_flag, 0);
+    EXPECT_EQ(lmp->atom->erforce_flag, 0);
+    EXPECT_EQ(lmp->atom->cs_flag, 0);
+    EXPECT_EQ(lmp->atom->csforce_flag, 0);
+    EXPECT_EQ(lmp->atom->vforce_flag, 0);
+    EXPECT_EQ(lmp->atom->ervelforce_flag, 0);
+    EXPECT_EQ(lmp->atom->etag_flag, 0);
+    EXPECT_EQ(lmp->atom->rho_flag, 0);
+    EXPECT_EQ(lmp->atom->esph_flag, 0);
+    EXPECT_EQ(lmp->atom->cv_flag, 0);
+    EXPECT_EQ(lmp->atom->vest_flag, 0);
+    EXPECT_EQ(lmp->atom->dpd_flag, 0);
+    EXPECT_EQ(lmp->atom->edpd_flag, 0);
+    EXPECT_EQ(lmp->atom->tdpd_flag, 0);
+    EXPECT_EQ(lmp->atom->mesont_flag, 0);
+    EXPECT_EQ(lmp->atom->sp_flag, 0);
+    EXPECT_EQ(lmp->atom->x0_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_flag, 0);
+    EXPECT_EQ(lmp->atom->damage_flag, 0);
+    EXPECT_EQ(lmp->atom->contact_radius_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_data_9_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_stress_flag, 0);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_flag, 0);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_rate_flag, 0);
+    EXPECT_EQ(lmp->atom->pdscale, 1.0);
+
+    EXPECT_NE(lmp->atom->tag, nullptr);
+    EXPECT_NE(lmp->atom->type, nullptr);
+    EXPECT_NE(lmp->atom->mask, nullptr);
+    EXPECT_NE(lmp->atom->image, nullptr);
+    EXPECT_NE(lmp->atom->x, nullptr);
+    EXPECT_NE(lmp->atom->v, nullptr);
+    EXPECT_NE(lmp->atom->f, nullptr);
+    EXPECT_EQ(lmp->atom->q, nullptr);
+    EXPECT_EQ(lmp->atom->mu, nullptr);
+    EXPECT_EQ(lmp->atom->omega, nullptr);
+    EXPECT_EQ(lmp->atom->angmom, nullptr);
+    EXPECT_EQ(lmp->atom->torque, nullptr);
+    EXPECT_EQ(lmp->atom->radius, nullptr);
+    EXPECT_EQ(lmp->atom->rmass, nullptr);
+    EXPECT_EQ(lmp->atom->ellipsoid, nullptr);
+    EXPECT_EQ(lmp->atom->line, nullptr);
+    EXPECT_EQ(lmp->atom->tri, nullptr);
+    EXPECT_EQ(lmp->atom->body, nullptr);
+    EXPECT_NE(lmp->atom->molecule, nullptr);
+    EXPECT_NE(lmp->atom->molindex, nullptr);
+    EXPECT_NE(lmp->atom->molatom, nullptr);
+    EXPECT_EQ(lmp->atom->num_bond, nullptr);
+    EXPECT_EQ(lmp->atom->bond_type, nullptr);
+    EXPECT_EQ(lmp->atom->bond_atom, nullptr);
+    EXPECT_EQ(lmp->atom->num_angle, nullptr);
+    EXPECT_EQ(lmp->atom->angle_type, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->num_dihedral, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_type, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom4, nullptr);
+    EXPECT_EQ(lmp->atom->num_improper, nullptr);
+    EXPECT_EQ(lmp->atom->improper_type, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom4, nullptr);
+    EXPECT_EQ(lmp->atom->maxspecial, 1);
+    EXPECT_EQ(lmp->atom->nspecial, nullptr);
+    EXPECT_EQ(lmp->atom->special, nullptr);
+    EXPECT_EQ(lmp->atom->vfrac, nullptr);
+    EXPECT_EQ(lmp->atom->s0, nullptr);
+    EXPECT_EQ(lmp->atom->x0, nullptr);
+    EXPECT_EQ(lmp->atom->sp, nullptr);
+    EXPECT_EQ(lmp->atom->fm, nullptr);
+    EXPECT_EQ(lmp->atom->fm_long, nullptr);
+    EXPECT_EQ(lmp->atom->spin, nullptr);
+    EXPECT_EQ(lmp->atom->eradius, nullptr);
+    EXPECT_EQ(lmp->atom->ervel, nullptr);
+    EXPECT_EQ(lmp->atom->erforce, nullptr);
+    EXPECT_EQ(lmp->atom->ervelforce, nullptr);
+    EXPECT_EQ(lmp->atom->cs, nullptr);
+    EXPECT_EQ(lmp->atom->csforce, nullptr);
+    EXPECT_EQ(lmp->atom->vforce, nullptr);
+    EXPECT_EQ(lmp->atom->etag, nullptr);
+    EXPECT_EQ(lmp->atom->uCond, nullptr);
+    EXPECT_EQ(lmp->atom->uMech, nullptr);
+    EXPECT_EQ(lmp->atom->uChem, nullptr);
+    EXPECT_EQ(lmp->atom->uCG, nullptr);
+    EXPECT_EQ(lmp->atom->uCGnew, nullptr);
+    EXPECT_EQ(lmp->atom->duChem, nullptr);
+    EXPECT_EQ(lmp->atom->dpdTheta, nullptr);
+    EXPECT_EQ(lmp->atom->cc, nullptr);
+    EXPECT_EQ(lmp->atom->cc_flux, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_temp, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_flux, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_cv, nullptr);
+    EXPECT_EQ(lmp->atom->length, nullptr);
+    EXPECT_EQ(lmp->atom->buckling, nullptr);
+    EXPECT_EQ(lmp->atom->bond_nt, nullptr);
+    EXPECT_EQ(lmp->atom->contact_radius, nullptr);
+    EXPECT_EQ(lmp->atom->smd_data_9, nullptr);
+    EXPECT_EQ(lmp->atom->smd_stress, nullptr);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain, nullptr);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_rate, nullptr);
+    EXPECT_EQ(lmp->atom->damage, nullptr);
+    EXPECT_EQ(lmp->atom->rho, nullptr);
+    EXPECT_EQ(lmp->atom->drho, nullptr);
+    EXPECT_EQ(lmp->atom->esph, nullptr);
+    EXPECT_EQ(lmp->atom->desph, nullptr);
+    EXPECT_EQ(lmp->atom->cv, nullptr);
+    EXPECT_EQ(lmp->atom->vest, nullptr);
+    EXPECT_EQ(lmp->atom->nmolecule, 2);
+    EXPECT_NE(lmp->atom->molecules, nullptr);
+    EXPECT_EQ(lmp->atom->nivector, 0);
+    EXPECT_EQ(lmp->atom->ndvector, 0);
+    EXPECT_EQ(lmp->atom->iname, nullptr);
+    EXPECT_EQ(lmp->atom->dname, nullptr);
+    EXPECT_EQ(lmp->atom->mass, nullptr);
+    EXPECT_EQ(lmp->atom->mass_setflag, nullptr);
+    EXPECT_EQ(lmp->atom->nextra_grow, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart, 0);
+    EXPECT_EQ(lmp->atom->nextra_border, 0);
+    EXPECT_EQ(lmp->atom->nextra_grow_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_border_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_store, 0);
+    EXPECT_EQ(lmp->atom->extra_grow, nullptr);
+    EXPECT_EQ(lmp->atom->extra_restart, nullptr);
+    EXPECT_EQ(lmp->atom->extra_border, nullptr);
+    EXPECT_EQ(lmp->atom->extra, nullptr);
+    EXPECT_EQ(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->map_style, 3);
+    EXPECT_EQ(lmp->atom->map_user, 0);
+    EXPECT_EQ(lmp->atom->map_tag_max, -1);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("create_box 4 box bond/types 2 angle/types 2 ");
+    lmp->input->one("create_atoms 0 single -2.0  2.0  0.1 mol twomols 65234");
+    lmp->input->one("create_atoms 0 single -2.0 -2.0 -0.1 mol twomols 62346");
+    lmp->input->one("create_atoms 0 single  2.0  2.0 -0.1 mol twomols 61354");
+    lmp->input->one("create_atoms 3 single  2.0 -2.0  0.1");
+    lmp->input->one("create_atoms 3 single  2.0  2.0 -2.1");
+    lmp->input->one("create_atoms 4 single  2.0 -2.0  2.1");
+    lmp->input->one("mass 1 16.0");
+    lmp->input->one("mass 2 1.0");
+    lmp->input->one("mass 3 12.0");
+    lmp->input->one("mass 4 16.0");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("bond_coeff 1 1.0");
+    lmp->input->one("bond_coeff 2 1.16");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("angle_coeff * 109.0");
+    lmp->input->one("pair_coeff * *");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("template"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nellipsoids, 0);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_EQ(lmp->atom->nextra_grow, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart, 0);
+    EXPECT_EQ(lmp->atom->nextra_border, 0);
+    EXPECT_EQ(lmp->atom->nextra_grow_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_border_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_store, 0);
+    EXPECT_EQ(lmp->atom->extra_grow, nullptr);
+    EXPECT_EQ(lmp->atom->extra_restart, nullptr);
+    EXPECT_EQ(lmp->atom->extra_border, nullptr);
+    EXPECT_EQ(lmp->atom->extra, nullptr);
+
+    EXPECT_NE(lmp->atom->mass, nullptr);
+    EXPECT_NE(lmp->atom->mass_setflag, nullptr);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("write_data test_atom_styles.data");
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("newton off");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style template twomols");
+    lmp->input->one("pair_style zero 4.0");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("atom_modify map array");
+    lmp->input->one("read_data test_atom_styles.data");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("template"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_NE(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    EXPECT_EQ(lmp->atom->map_user, 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    auto molecule = lmp->atom->molecule;
+    auto molindex = lmp->atom->molindex;
+    auto molatom  = lmp->atom->molatom;
+
+    EXPECT_EQ(molecule[GETIDX(1)], 1);
+    EXPECT_EQ(molecule[GETIDX(2)], 1);
+    EXPECT_EQ(molecule[GETIDX(3)], 1);
+    EXPECT_EQ(molecule[GETIDX(4)], 2);
+    EXPECT_EQ(molecule[GETIDX(5)], 2);
+    EXPECT_EQ(molecule[GETIDX(6)], 2);
+    EXPECT_EQ(molecule[GETIDX(7)], 3);
+    EXPECT_EQ(molecule[GETIDX(8)], 3);
+    EXPECT_EQ(molecule[GETIDX(9)], 3);
+    EXPECT_EQ(molecule[GETIDX(10)], 0);
+    EXPECT_EQ(molecule[GETIDX(11)], 0);
+    EXPECT_EQ(molecule[GETIDX(12)], 0);
+    EXPECT_EQ(molindex[GETIDX(1)], 0);
+    EXPECT_EQ(molindex[GETIDX(2)], 0);
+    EXPECT_EQ(molindex[GETIDX(3)], 0);
+    EXPECT_EQ(molindex[GETIDX(4)], 0);
+    EXPECT_EQ(molindex[GETIDX(5)], 0);
+    EXPECT_EQ(molindex[GETIDX(6)], 0);
+    EXPECT_EQ(molindex[GETIDX(7)], 0);
+    EXPECT_EQ(molindex[GETIDX(8)], 0);
+    EXPECT_EQ(molindex[GETIDX(9)], 0);
+    EXPECT_EQ(molindex[GETIDX(10)], -1);
+    EXPECT_EQ(molindex[GETIDX(11)], -1);
+    EXPECT_EQ(molindex[GETIDX(12)], -1);
+    EXPECT_EQ(molatom[GETIDX(1)], 0);
+    EXPECT_EQ(molatom[GETIDX(2)], 1);
+    EXPECT_EQ(molatom[GETIDX(3)], 2);
+    EXPECT_EQ(molatom[GETIDX(4)], 0);
+    EXPECT_EQ(molatom[GETIDX(5)], 1);
+    EXPECT_EQ(molatom[GETIDX(6)], 2);
+    EXPECT_EQ(molatom[GETIDX(7)], 0);
+    EXPECT_EQ(molatom[GETIDX(8)], 1);
+    EXPECT_EQ(molatom[GETIDX(9)], 2);
+    EXPECT_EQ(molatom[GETIDX(10)], -1);
+    EXPECT_EQ(molatom[GETIDX(11)], -1);
+    EXPECT_EQ(molatom[GETIDX(12)], -1);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style template twomols");
+    lmp->input->one("pair_style zero 4.0");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("atom_modify map array");
+    lmp->input->one("read_data test_atom_styles.data");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("template"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_NE(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    EXPECT_EQ(lmp->atom->map_user, 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+
+    EXPECT_EQ(molindex[GETIDX(1)], 0);
+    EXPECT_EQ(molindex[GETIDX(2)], 0);
+    EXPECT_EQ(molindex[GETIDX(3)], 0);
+    EXPECT_EQ(molindex[GETIDX(4)], 0);
+    EXPECT_EQ(molindex[GETIDX(5)], 0);
+    EXPECT_EQ(molindex[GETIDX(6)], 0);
+    EXPECT_EQ(molindex[GETIDX(7)], 0);
+    EXPECT_EQ(molindex[GETIDX(8)], 0);
+    EXPECT_EQ(molindex[GETIDX(9)], 0);
+    EXPECT_EQ(molindex[GETIDX(10)], -1);
+    EXPECT_EQ(molindex[GETIDX(11)], -1);
+    EXPECT_EQ(molindex[GETIDX(12)], -1);
+    EXPECT_EQ(molatom[GETIDX(1)], 0);
+    EXPECT_EQ(molatom[GETIDX(2)], 1);
+    EXPECT_EQ(molatom[GETIDX(3)], 2);
+    EXPECT_EQ(molatom[GETIDX(4)], 0);
+    EXPECT_EQ(molatom[GETIDX(5)], 1);
+    EXPECT_EQ(molatom[GETIDX(6)], 2);
+    EXPECT_EQ(molatom[GETIDX(7)], 0);
+    EXPECT_EQ(molatom[GETIDX(8)], 1);
+    EXPECT_EQ(molatom[GETIDX(9)], 2);
+    EXPECT_EQ(molatom[GETIDX(10)], -1);
+    EXPECT_EQ(molatom[GETIDX(11)], -1);
+    EXPECT_EQ(molatom[GETIDX(12)], -1);
+
+    auto x    = lmp->atom->x;
+    auto v    = lmp->atom->v;
+    auto type = lmp->atom->type;
+
+    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_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(7)], 1);
+    EXPECT_EQ(type[GETIDX(8)], 2);
+    EXPECT_EQ(type[GETIDX(9)], 2);
+    EXPECT_EQ(type[GETIDX(10)], 3);
+    EXPECT_EQ(type[GETIDX(11)], 3);
+    EXPECT_EQ(type[GETIDX(12)], 4);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("group two id 7:10");
+    lmp->input->one("delete_atoms group two compress no");
+    lmp->input->one("write_restart test_atom_styles.restart");
+    lmp->input->one("clear");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("read_restart test_atom_styles.restart");
+    // FIXME.
+    // lmp->input->one("replicate 1 1 2 bbox");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("template"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+    EXPECT_EQ(lmp->atom->natoms, 8);
+    EXPECT_EQ(lmp->atom->nlocal, 8);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->q_flag, 0);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 0);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    type     = lmp->atom->type;
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+    EXPECT_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(11)], 3);
+    EXPECT_EQ(type[GETIDX(12)], 4);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("reset_atom_ids");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 8);
+
+    type = lmp->atom->type;
+    EXPECT_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(7)], 4);
+    EXPECT_EQ(type[GETIDX(8)], 3);
+}
+
+TEST_F(AtomStyleTest, template_charge)
+{
+    if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();
+    create_molecule_files();
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style hybrid template twomols charge");
+    lmp->input->one("newton on");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    auto hybrid = (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(lmp->atom->avec, nullptr);
+    ASSERT_NE(hybrid->styles[0], nullptr);
+    ASSERT_NE(hybrid->styles[1], nullptr);
+
+    EXPECT_EQ(lmp->atom->natoms, 0);
+    EXPECT_EQ(lmp->atom->nlocal, 0);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_EQ(lmp->atom->nmax, 1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->nellipsoids, 0);
+    EXPECT_EQ(lmp->atom->nlines, 0);
+    EXPECT_EQ(lmp->atom->ntris, 0);
+    EXPECT_EQ(lmp->atom->nbodies, 0);
+    EXPECT_EQ(lmp->atom->nbonds, 0);
+    EXPECT_EQ(lmp->atom->nangles, 0);
+    EXPECT_EQ(lmp->atom->ndihedrals, 0);
+    EXPECT_EQ(lmp->atom->nimpropers, 0);
+    EXPECT_EQ(lmp->atom->ntypes, 0);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->ndihedraltypes, 0);
+    EXPECT_EQ(lmp->atom->nimpropertypes, 0);
+    EXPECT_EQ(lmp->atom->bond_per_atom, 0);
+    EXPECT_EQ(lmp->atom->angle_per_atom, 0);
+    EXPECT_EQ(lmp->atom->dihedral_per_atom, 0);
+    EXPECT_EQ(lmp->atom->improper_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_bond_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_angle_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_dihedral_per_atom, 0);
+    EXPECT_EQ(lmp->atom->extra_improper_per_atom, 0);
+
+    EXPECT_EQ(lmp->atom->sphere_flag, 0);
+    EXPECT_EQ(lmp->atom->ellipsoid_flag, 0);
+    EXPECT_EQ(lmp->atom->line_flag, 0);
+    EXPECT_EQ(lmp->atom->tri_flag, 0);
+    EXPECT_EQ(lmp->atom->body_flag, 0);
+    EXPECT_EQ(lmp->atom->peri_flag, 0);
+    EXPECT_EQ(lmp->atom->electron_flag, 0);
+    EXPECT_EQ(lmp->atom->wavepacket_flag, 0);
+    EXPECT_EQ(lmp->atom->sph_flag, 0);
+    EXPECT_EQ(lmp->atom->molecule_flag, 1);
+    EXPECT_EQ(lmp->atom->molindex_flag, 0);
+    EXPECT_EQ(lmp->atom->molatom_flag, 0);
+    EXPECT_EQ(lmp->atom->q_flag, 1);
+    EXPECT_EQ(lmp->atom->mu_flag, 0);
+    EXPECT_EQ(lmp->atom->rmass_flag, 0);
+    EXPECT_EQ(lmp->atom->radius_flag, 0);
+    EXPECT_EQ(lmp->atom->omega_flag, 0);
+    EXPECT_EQ(lmp->atom->torque_flag, 0);
+    EXPECT_EQ(lmp->atom->angmom_flag, 0);
+    EXPECT_EQ(lmp->atom->vfrac_flag, 0);
+    EXPECT_EQ(lmp->atom->spin_flag, 0);
+    EXPECT_EQ(lmp->atom->eradius_flag, 0);
+    EXPECT_EQ(lmp->atom->ervel_flag, 0);
+    EXPECT_EQ(lmp->atom->erforce_flag, 0);
+    EXPECT_EQ(lmp->atom->cs_flag, 0);
+    EXPECT_EQ(lmp->atom->csforce_flag, 0);
+    EXPECT_EQ(lmp->atom->vforce_flag, 0);
+    EXPECT_EQ(lmp->atom->ervelforce_flag, 0);
+    EXPECT_EQ(lmp->atom->etag_flag, 0);
+    EXPECT_EQ(lmp->atom->rho_flag, 0);
+    EXPECT_EQ(lmp->atom->esph_flag, 0);
+    EXPECT_EQ(lmp->atom->cv_flag, 0);
+    EXPECT_EQ(lmp->atom->vest_flag, 0);
+    EXPECT_EQ(lmp->atom->dpd_flag, 0);
+    EXPECT_EQ(lmp->atom->edpd_flag, 0);
+    EXPECT_EQ(lmp->atom->tdpd_flag, 0);
+    EXPECT_EQ(lmp->atom->mesont_flag, 0);
+    EXPECT_EQ(lmp->atom->sp_flag, 0);
+    EXPECT_EQ(lmp->atom->x0_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_flag, 0);
+    EXPECT_EQ(lmp->atom->damage_flag, 0);
+    EXPECT_EQ(lmp->atom->contact_radius_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_data_9_flag, 0);
+    EXPECT_EQ(lmp->atom->smd_stress_flag, 0);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_flag, 0);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_rate_flag, 0);
+    EXPECT_EQ(lmp->atom->pdscale, 1.0);
+
+    EXPECT_NE(lmp->atom->tag, nullptr);
+    EXPECT_NE(lmp->atom->type, nullptr);
+    EXPECT_NE(lmp->atom->mask, nullptr);
+    EXPECT_NE(lmp->atom->image, nullptr);
+    EXPECT_NE(lmp->atom->x, nullptr);
+    EXPECT_NE(lmp->atom->v, nullptr);
+    EXPECT_NE(lmp->atom->f, nullptr);
+    EXPECT_NE(lmp->atom->q, nullptr);
+    EXPECT_EQ(lmp->atom->mu, nullptr);
+    EXPECT_EQ(lmp->atom->omega, nullptr);
+    EXPECT_EQ(lmp->atom->angmom, nullptr);
+    EXPECT_EQ(lmp->atom->torque, nullptr);
+    EXPECT_EQ(lmp->atom->radius, nullptr);
+    EXPECT_EQ(lmp->atom->rmass, nullptr);
+    EXPECT_EQ(lmp->atom->ellipsoid, nullptr);
+    EXPECT_EQ(lmp->atom->line, nullptr);
+    EXPECT_EQ(lmp->atom->tri, nullptr);
+    EXPECT_EQ(lmp->atom->body, nullptr);
+    EXPECT_NE(lmp->atom->molecule, nullptr);
+    EXPECT_NE(lmp->atom->molindex, nullptr);
+    EXPECT_NE(lmp->atom->molatom, nullptr);
+    EXPECT_EQ(lmp->atom->num_bond, nullptr);
+    EXPECT_EQ(lmp->atom->bond_type, nullptr);
+    EXPECT_EQ(lmp->atom->bond_atom, nullptr);
+    EXPECT_EQ(lmp->atom->num_angle, nullptr);
+    EXPECT_EQ(lmp->atom->angle_type, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->angle_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->num_dihedral, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_type, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->dihedral_atom4, nullptr);
+    EXPECT_EQ(lmp->atom->num_improper, nullptr);
+    EXPECT_EQ(lmp->atom->improper_type, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom1, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom2, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom3, nullptr);
+    EXPECT_EQ(lmp->atom->improper_atom4, nullptr);
+    EXPECT_EQ(lmp->atom->maxspecial, 1);
+    EXPECT_EQ(lmp->atom->nspecial, nullptr);
+    EXPECT_EQ(lmp->atom->special, nullptr);
+    EXPECT_EQ(lmp->atom->vfrac, nullptr);
+    EXPECT_EQ(lmp->atom->s0, nullptr);
+    EXPECT_EQ(lmp->atom->x0, nullptr);
+    EXPECT_EQ(lmp->atom->sp, nullptr);
+    EXPECT_EQ(lmp->atom->fm, nullptr);
+    EXPECT_EQ(lmp->atom->fm_long, nullptr);
+    EXPECT_EQ(lmp->atom->spin, nullptr);
+    EXPECT_EQ(lmp->atom->eradius, nullptr);
+    EXPECT_EQ(lmp->atom->ervel, nullptr);
+    EXPECT_EQ(lmp->atom->erforce, nullptr);
+    EXPECT_EQ(lmp->atom->ervelforce, nullptr);
+    EXPECT_EQ(lmp->atom->cs, nullptr);
+    EXPECT_EQ(lmp->atom->csforce, nullptr);
+    EXPECT_EQ(lmp->atom->vforce, nullptr);
+    EXPECT_EQ(lmp->atom->etag, nullptr);
+    EXPECT_EQ(lmp->atom->uCond, nullptr);
+    EXPECT_EQ(lmp->atom->uMech, nullptr);
+    EXPECT_EQ(lmp->atom->uChem, nullptr);
+    EXPECT_EQ(lmp->atom->uCG, nullptr);
+    EXPECT_EQ(lmp->atom->uCGnew, nullptr);
+    EXPECT_EQ(lmp->atom->duChem, nullptr);
+    EXPECT_EQ(lmp->atom->dpdTheta, nullptr);
+    EXPECT_EQ(lmp->atom->cc, nullptr);
+    EXPECT_EQ(lmp->atom->cc_flux, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_temp, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_flux, nullptr);
+    EXPECT_EQ(lmp->atom->edpd_cv, nullptr);
+    EXPECT_EQ(lmp->atom->length, nullptr);
+    EXPECT_EQ(lmp->atom->buckling, nullptr);
+    EXPECT_EQ(lmp->atom->bond_nt, nullptr);
+    EXPECT_EQ(lmp->atom->contact_radius, nullptr);
+    EXPECT_EQ(lmp->atom->smd_data_9, nullptr);
+    EXPECT_EQ(lmp->atom->smd_stress, nullptr);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain, nullptr);
+    EXPECT_EQ(lmp->atom->eff_plastic_strain_rate, nullptr);
+    EXPECT_EQ(lmp->atom->damage, nullptr);
+    EXPECT_EQ(lmp->atom->rho, nullptr);
+    EXPECT_EQ(lmp->atom->drho, nullptr);
+    EXPECT_EQ(lmp->atom->esph, nullptr);
+    EXPECT_EQ(lmp->atom->desph, nullptr);
+    EXPECT_EQ(lmp->atom->cv, nullptr);
+    EXPECT_EQ(lmp->atom->vest, nullptr);
+    EXPECT_EQ(lmp->atom->nmolecule, 2);
+    EXPECT_NE(lmp->atom->molecules, nullptr);
+    EXPECT_EQ(lmp->atom->nivector, 0);
+    EXPECT_EQ(lmp->atom->ndvector, 0);
+    EXPECT_EQ(lmp->atom->iname, nullptr);
+    EXPECT_EQ(lmp->atom->dname, nullptr);
+    EXPECT_EQ(lmp->atom->mass, nullptr);
+    EXPECT_EQ(lmp->atom->mass_setflag, nullptr);
+    EXPECT_EQ(lmp->atom->nextra_grow, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart, 0);
+    EXPECT_EQ(lmp->atom->nextra_border, 0);
+    EXPECT_EQ(lmp->atom->nextra_grow_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_border_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_store, 0);
+    EXPECT_EQ(lmp->atom->extra_grow, nullptr);
+    EXPECT_EQ(lmp->atom->extra_restart, nullptr);
+    EXPECT_EQ(lmp->atom->extra_border, nullptr);
+    EXPECT_EQ(lmp->atom->extra, nullptr);
+    EXPECT_EQ(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->map_style, 3);
+    EXPECT_EQ(lmp->atom->map_user, 0);
+    EXPECT_EQ(lmp->atom->map_tag_max, -1);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("create_box 4 box bond/types 2 angle/types 2 ");
+    lmp->input->one("create_atoms 0 single -2.0  2.0  0.1 mol twomols 65234");
+    lmp->input->one("create_atoms 0 single -2.0 -2.0 -0.1 mol twomols 62346");
+    lmp->input->one("create_atoms 0 single  2.0  2.0 -0.1 mol twomols 61354");
+    lmp->input->one("create_atoms 3 single  2.0 -2.0  0.1");
+    lmp->input->one("create_atoms 3 single  2.0  2.0 -2.1");
+    lmp->input->one("create_atoms 4 single  2.0 -2.0  2.1");
+    lmp->input->one("mass 1 16.0");
+    lmp->input->one("mass 2 1.0");
+    lmp->input->one("mass 3 12.0");
+    lmp->input->one("mass 4 16.0");
+    lmp->input->one("set atom 10 charge 0.7");
+    lmp->input->one("set atom 11 charge -0.35");
+    lmp->input->one("set atom 12 charge -0.35");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("bond_coeff 1 1.0");
+    lmp->input->one("bond_coeff 2 1.16");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("angle_coeff * 109.0");
+    lmp->input->one("pair_coeff * *");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    ASSERT_NE(lmp->atom->avec, nullptr);
+    hybrid = (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);
+
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nellipsoids, 0);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->q_flag, 1);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_EQ(lmp->atom->nextra_grow, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart, 0);
+    EXPECT_EQ(lmp->atom->nextra_border, 0);
+    EXPECT_EQ(lmp->atom->nextra_grow_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_restart_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_border_max, 0);
+    EXPECT_EQ(lmp->atom->nextra_store, 0);
+    EXPECT_EQ(lmp->atom->extra_grow, nullptr);
+    EXPECT_EQ(lmp->atom->extra_restart, nullptr);
+    EXPECT_EQ(lmp->atom->extra_border, nullptr);
+    EXPECT_EQ(lmp->atom->extra, nullptr);
+
+    EXPECT_NE(lmp->atom->mass, nullptr);
+    EXPECT_NE(lmp->atom->mass_setflag, nullptr);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("write_data test_atom_styles.data");
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("newton off");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style hybrid template twomols charge");
+    lmp->input->one("pair_style zero 4.0");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("atom_modify map array");
+    lmp->input->one("read_data test_atom_styles.data");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_NE(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    EXPECT_EQ(lmp->atom->map_user, 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    auto molecule = lmp->atom->molecule;
+    auto molindex = lmp->atom->molindex;
+    auto molatom  = lmp->atom->molatom;
+
+    EXPECT_EQ(molecule[GETIDX(1)], 1);
+    EXPECT_EQ(molecule[GETIDX(2)], 1);
+    EXPECT_EQ(molecule[GETIDX(3)], 1);
+    EXPECT_EQ(molecule[GETIDX(4)], 2);
+    EXPECT_EQ(molecule[GETIDX(5)], 2);
+    EXPECT_EQ(molecule[GETIDX(6)], 2);
+    EXPECT_EQ(molecule[GETIDX(7)], 3);
+    EXPECT_EQ(molecule[GETIDX(8)], 3);
+    EXPECT_EQ(molecule[GETIDX(9)], 3);
+    EXPECT_EQ(molecule[GETIDX(10)], 0);
+    EXPECT_EQ(molecule[GETIDX(11)], 0);
+    EXPECT_EQ(molecule[GETIDX(12)], 0);
+    EXPECT_EQ(molindex[GETIDX(1)], 0);
+    EXPECT_EQ(molindex[GETIDX(2)], 0);
+    EXPECT_EQ(molindex[GETIDX(3)], 0);
+    EXPECT_EQ(molindex[GETIDX(4)], 0);
+    EXPECT_EQ(molindex[GETIDX(5)], 0);
+    EXPECT_EQ(molindex[GETIDX(6)], 0);
+    EXPECT_EQ(molindex[GETIDX(7)], 0);
+    EXPECT_EQ(molindex[GETIDX(8)], 0);
+    EXPECT_EQ(molindex[GETIDX(9)], 0);
+    EXPECT_EQ(molindex[GETIDX(10)], -1);
+    EXPECT_EQ(molindex[GETIDX(11)], -1);
+    EXPECT_EQ(molindex[GETIDX(12)], -1);
+    EXPECT_EQ(molatom[GETIDX(1)], 0);
+    EXPECT_EQ(molatom[GETIDX(2)], 1);
+    EXPECT_EQ(molatom[GETIDX(3)], 2);
+    EXPECT_EQ(molatom[GETIDX(4)], 0);
+    EXPECT_EQ(molatom[GETIDX(5)], 1);
+    EXPECT_EQ(molatom[GETIDX(6)], 2);
+    EXPECT_EQ(molatom[GETIDX(7)], 0);
+    EXPECT_EQ(molatom[GETIDX(8)], 1);
+    EXPECT_EQ(molatom[GETIDX(9)], 2);
+    EXPECT_EQ(molatom[GETIDX(10)], -1);
+    EXPECT_EQ(molatom[GETIDX(11)], -1);
+    EXPECT_EQ(molatom[GETIDX(12)], -1);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("atom_style hybrid template twomols charge");
+    lmp->input->one("pair_style zero 4.0");
+    lmp->input->one("bond_style zero");
+    lmp->input->one("angle_style zero");
+    lmp->input->one("atom_modify map array");
+    lmp->input->one("read_data test_atom_styles.data");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+
+    EXPECT_EQ(lmp->atom->natoms, 12);
+    EXPECT_EQ(lmp->atom->nlocal, 12);
+    EXPECT_EQ(lmp->atom->nbonds, 6);
+    EXPECT_EQ(lmp->atom->nangles, 3);
+    EXPECT_EQ(lmp->atom->nbondtypes, 2);
+    EXPECT_EQ(lmp->atom->nangletypes, 2);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_NE(lmp->atom->sametag, nullptr);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_style, Atom::MAP_ARRAY);
+    EXPECT_EQ(lmp->atom->map_user, 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+
+    EXPECT_EQ(molindex[GETIDX(1)], 0);
+    EXPECT_EQ(molindex[GETIDX(2)], 0);
+    EXPECT_EQ(molindex[GETIDX(3)], 0);
+    EXPECT_EQ(molindex[GETIDX(4)], 0);
+    EXPECT_EQ(molindex[GETIDX(5)], 0);
+    EXPECT_EQ(molindex[GETIDX(6)], 0);
+    EXPECT_EQ(molindex[GETIDX(7)], 0);
+    EXPECT_EQ(molindex[GETIDX(8)], 0);
+    EXPECT_EQ(molindex[GETIDX(9)], 0);
+    EXPECT_EQ(molindex[GETIDX(10)], -1);
+    EXPECT_EQ(molindex[GETIDX(11)], -1);
+    EXPECT_EQ(molindex[GETIDX(12)], -1);
+    EXPECT_EQ(molatom[GETIDX(1)], 0);
+    EXPECT_EQ(molatom[GETIDX(2)], 1);
+    EXPECT_EQ(molatom[GETIDX(3)], 2);
+    EXPECT_EQ(molatom[GETIDX(4)], 0);
+    EXPECT_EQ(molatom[GETIDX(5)], 1);
+    EXPECT_EQ(molatom[GETIDX(6)], 2);
+    EXPECT_EQ(molatom[GETIDX(7)], 0);
+    EXPECT_EQ(molatom[GETIDX(8)], 1);
+    EXPECT_EQ(molatom[GETIDX(9)], 2);
+    EXPECT_EQ(molatom[GETIDX(10)], -1);
+    EXPECT_EQ(molatom[GETIDX(11)], -1);
+    EXPECT_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);
+    EXPECT_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(7)], 1);
+    EXPECT_EQ(type[GETIDX(8)], 2);
+    EXPECT_EQ(type[GETIDX(9)], 2);
+    EXPECT_EQ(type[GETIDX(10)], 3);
+    EXPECT_EQ(type[GETIDX(11)], 3);
+    EXPECT_EQ(type[GETIDX(12)], 4);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("group two id 7:10");
+    lmp->input->one("delete_atoms group two compress no");
+    lmp->input->one("write_restart test_atom_styles.restart");
+    lmp->input->one("clear");
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    lmp->input->one("read_restart test_atom_styles.restart");
+    // FIXME.
+    // lmp->input->one("replicate 1 1 2 bbox");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_THAT(std::string(lmp->atom->atom_style), Eq("hybrid"));
+    EXPECT_NE(lmp->atom->avec, nullptr);
+    EXPECT_EQ(lmp->atom->natoms, 8);
+    EXPECT_EQ(lmp->atom->nlocal, 8);
+    EXPECT_EQ(lmp->atom->nghost, 0);
+    EXPECT_NE(lmp->atom->nmax, -1);
+    EXPECT_EQ(lmp->atom->tag_enable, 1);
+    EXPECT_EQ(lmp->atom->molecular, Atom::TEMPLATE);
+    EXPECT_EQ(lmp->atom->ntypes, 4);
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 0);
+    EXPECT_EQ(lmp->atom->map_tag_max, 12);
+
+    type     = lmp->atom->type;
+    molecule = lmp->atom->molecule;
+    molindex = lmp->atom->molindex;
+    molatom  = lmp->atom->molatom;
+    EXPECT_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(11)], 3);
+    EXPECT_EQ(type[GETIDX(12)], 4);
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("reset_atom_ids");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+    EXPECT_EQ(lmp->atom->tag_consecutive(), 1);
+    EXPECT_EQ(lmp->atom->map_tag_max, 8);
+
+    type = lmp->atom->type;
+    EXPECT_EQ(type[GETIDX(1)], 1);
+    EXPECT_EQ(type[GETIDX(2)], 2);
+    EXPECT_EQ(type[GETIDX(3)], 2);
+    EXPECT_EQ(type[GETIDX(4)], 1);
+    EXPECT_EQ(type[GETIDX(5)], 2);
+    EXPECT_EQ(type[GETIDX(6)], 2);
+    EXPECT_EQ(type[GETIDX(7)], 4);
+    EXPECT_EQ(type[GETIDX(8)], 3);
+}
+
 TEST_F(AtomStyleTest, bond)
 {
     if (!LAMMPS::is_installed_pkg("MOLECULE")) GTEST_SKIP();