diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp
index 6ab1508c..3b73cd34 100644
--- a/src/CabanaPD_Comm.hpp
+++ b/src/CabanaPD_Comm.hpp
@@ -230,12 +230,13 @@ struct HaloIds
     }
 };
 
-template <class ParticleType, class ModelType, class ThermalType>
+template <class ParticleType, class ModelType, class SpeciesType,
+          class ThermalType>
 class Comm;
 
 // FIXME: extract model from ParticleType instead.
 template <class ParticleType>
-class Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
     int mpi_size = -1;
@@ -320,11 +321,12 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
 };
 
 template <class ParticleType>
-class Comm<ParticleType, LPS, TemperatureIndependent>
-    : public Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, LPS, SingleSpecies, TemperatureIndependent>
+    : public Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
-    using base_type = Comm<ParticleType, PMB, TemperatureIndependent>;
+    using base_type =
+        Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>;
     using memory_space = typename base_type::memory_space;
     using halo_type = typename base_type::halo_type;
     using base_type::gather_u;
@@ -368,12 +370,52 @@ class Comm<ParticleType, LPS, TemperatureIndependent>
     }
 };
 
+template <class ParticleType, class ModelType>
+class Comm<ParticleType, ModelType, MultiSpecies, TemperatureIndependent>
+    : public Comm<ParticleType, ModelType, SingleSpecies,
+                  TemperatureIndependent>
+{
+  public:
+    using base_type =
+        Comm<ParticleType, ModelType, SingleSpecies, TemperatureIndependent>;
+    using memory_space = typename base_type::memory_space;
+    using halo_type = typename base_type::halo_type;
+    using base_type::halo;
+
+    using base_type::_init_timer;
+    using base_type::_timer;
+
+    using gather_species_type =
+        Cabana::Gather<halo_type, typename ParticleType::aosoa_species_type>;
+    std::shared_ptr<gather_species_type> gather_species;
+
+    Comm( ParticleType& particles, int max_export_guess = 100 )
+        : base_type( particles, max_export_guess )
+    {
+        _init_timer.start();
+
+        gather_species = std::make_shared<gather_species_type>(
+            *halo, particles._aosoa_species );
+
+        particles.resize( halo->numLocal(), halo->numGhost() );
+        _init_timer.stop();
+    }
+
+    void gatherSpecies()
+    {
+        _timer.start();
+        gather_species->apply();
+        _timer.stop();
+    }
+};
+
 template <class ParticleType>
-class Comm<ParticleType, PMB, TemperatureDependent>
-    : public Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, PMB, SingleSpecies, TemperatureDependent>
+    : public Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
-    using base_type = Comm<ParticleType, PMB, TemperatureIndependent>;
+    using base_type =
+        Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>;
     using memory_space = typename base_type::memory_space;
     using halo_type = typename base_type::halo_type;
     using base_type::halo;
diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp
index 83ed009c..42d456e6 100644
--- a/src/CabanaPD_ForceModels.hpp
+++ b/src/CabanaPD_ForceModels.hpp
@@ -23,6 +23,7 @@ struct BaseForceModel;
 template <>
 struct BaseForceModel<>
 {
+    using species_type = SingleSpecies;
     double delta;
 
     BaseForceModel(){};
@@ -40,6 +41,8 @@ struct BaseForceModel<>
 template <typename... MemorySpace>
 struct BaseForceModel
 {
+    using species_type = MultiSpecies;
+
     // Only allow one memory space.
     using memory_space =
         typename std::tuple_element<0, std::tuple<MemorySpace...>>::type;
@@ -86,6 +89,7 @@ struct BaseForceModel
 template <typename TemperatureType>
 struct BaseTemperatureModel
 {
+    using species_type = SingleSpecies;
     using memory_space = typename TemperatureType::memory_space;
     double alpha;
     double temp0;
diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
index c1151770..ec4f2cf7 100644
--- a/src/CabanaPD_Particles.hpp
+++ b/src/CabanaPD_Particles.hpp
@@ -101,12 +101,11 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     using vector_type = Cabana::MemberTypes<double[dim]>;
     // volume, dilatation, weighted_volume.
     using scalar_type = Cabana::MemberTypes<double>;
-    // no-fail.
+    // no-fail, type.
     using int_type = Cabana::MemberTypes<int>;
-    // v, W, rho, damage,  type.
+    // v, W, rho, damage.
     using other_types =
-        Cabana::MemberTypes<double[dim], double, double, double, int>;
-    // Potentially needed later: body force (b), ID.
+        Cabana::MemberTypes<double[dim], double, double, double>;
 
     // FIXME: add vector length.
     // FIXME: enable variable aosoa.
@@ -114,6 +113,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     using aosoa_y_type = Cabana::AoSoA<vector_type, memory_space, 1>;
     using aosoa_vol_type = Cabana::AoSoA<scalar_type, memory_space, 1>;
     using aosoa_nofail_type = Cabana::AoSoA<int_type, memory_space, 1>;
+    using aosoa_species_type = Cabana::AoSoA<int_type, memory_space, 1>;
     using aosoa_other_type = Cabana::AoSoA<other_types, memory_space>;
     // Using grid here for the particle init.
     using plist_x_type =
@@ -356,8 +356,11 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     {
         return Cabana::slice<0>( _aosoa_vol, "volume" );
     }
-    auto sliceType() { return Cabana::slice<4>( _aosoa_other, "type" ); }
-    auto sliceType() const { return Cabana::slice<4>( _aosoa_other, "type" ); }
+    auto sliceType() { return Cabana::slice<0>( _aosoa_species, "type" ); }
+    auto sliceType() const
+    {
+        return Cabana::slice<0>( _aosoa_species, "type" );
+    }
     auto sliceStrainEnergy()
     {
         return Cabana::slice<1>( _aosoa_other, "strain_energy" );
@@ -427,6 +430,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
         _plist_f.aosoa().resize( new_local );
         _aosoa_other.resize( new_local );
         _aosoa_nofail.resize( new_local + new_ghost );
+        _aosoa_species.resize( new_local + new_ghost );
         size = _plist_x.size();
         _timer.stop();
     };
@@ -457,7 +461,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
                 "particles", local_grid->globalGrid(), output_step, output_time,
                 0, n_local, getPosition( use_reference ), sliceStrainEnergy(),
                 sliceForce(), sliceDisplacement(), sliceVelocity(),
-                sliceDamage() );
+                sliceDamage(), sliceType() );
 #else
         log( std::cout, "No particle output enabled." );
 #endif
@@ -470,14 +474,17 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     auto timeOutput() { return _output_timer.time(); };
     auto time() { return _timer.time(); };
 
-    friend class Comm<self_type, PMB, TemperatureIndependent>;
-    friend class Comm<self_type, PMB, TemperatureDependent>;
+    friend class Comm<self_type, PMB, SingleSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, SingleSpecies, TemperatureDependent>;
+    friend class Comm<self_type, PMB, MultiSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, MultiSpecies, TemperatureDependent>;
 
   protected:
     aosoa_u_type _aosoa_u;
     aosoa_y_type _aosoa_y;
     aosoa_vol_type _aosoa_vol;
     aosoa_nofail_type _aosoa_nofail;
+    aosoa_species_type _aosoa_species;
     aosoa_other_type _aosoa_other;
 
     plist_x_type _plist_x;
@@ -619,7 +626,7 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, Dimension>
                 base_type::sliceStrainEnergy(), base_type::sliceForce(),
                 base_type::sliceDisplacement(), base_type::sliceVelocity(),
                 base_type::sliceDamage(), sliceWeightedVolume(),
-                sliceDilatation() );
+                sliceDilatation(), base_type::sliceType() );
 #else
         log( std::cout, "No particle output enabled." );
 #endif
@@ -628,8 +635,10 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, Dimension>
         _output_timer.stop();
     }
 
-    friend class Comm<self_type, PMB, TemperatureIndependent>;
-    friend class Comm<self_type, LPS, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, SingleSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, LPS, SingleSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, MultiSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, LPS, MultiSpecies, TemperatureIndependent>;
 
   protected:
     void init_lps()
@@ -761,10 +770,10 @@ class Particles<MemorySpace, PMB, TemperatureDependent, Dimension>
 #endif
     }
 
-    friend class Comm<self_type, PMB, TemperatureIndependent>;
-    friend class Comm<self_type, LPS, TemperatureIndependent>;
-    friend class Comm<self_type, PMB, TemperatureDependent>;
-    friend class Comm<self_type, LPS, TemperatureDependent>;
+    friend class Comm<self_type, PMB, SingleSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, SingleSpecies, TemperatureDependent>;
+    friend class Comm<self_type, PMB, MultiSpecies, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, MultiSpecies, TemperatureDependent>;
 
   protected:
     void init_temp()
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index d01b00e0..acfa379e 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -103,6 +103,7 @@ class SolverElastic
     using force_model_type = ForceModel;
     using force_type = Force<exec_space, force_model_type>;
     using comm_type = Comm<particle_type, typename force_model_type::base_model,
+                           typename force_model_type::species_type,
                            typename force_model_type::thermal_type>;
     using neighbor_type =
         Cabana::VerletList<memory_space, Cabana::FullNeighborTag,
@@ -128,7 +129,10 @@ class SolverElastic
         // Add ghosts from other MPI ranks.
         comm = std::make_shared<comm_type>( *particles );
 
-        // Update temperature ghost size if needed.
+        // Update optional property ghost sizes if needed.
+        if constexpr ( std::is_same<typename force_model_type::species_type,
+                                    MultiSpecies>::value )
+            force_model.update( particles->sliceType() );
         if constexpr ( std::is_same<typename force_model_type::thermal_type,
                                     TemperatureDependent>::value )
             force_model.update( particles->sliceTemperature() );
@@ -213,7 +217,10 @@ class SolverElastic
         if ( !boundary_condition.forceUpdate() )
             boundary_condition.apply( exec_space(), *particles, 0.0 );
 
-        // Communicate temperature.
+        // Communicate optional properties.
+        if constexpr ( std::is_same<typename force_model_type::species_type,
+                                    MultiSpecies>::value )
+            comm->gatherSpecies();
         if constexpr ( std::is_same<typename force_model_type::thermal_type,
                                     TemperatureDependent>::value )
             comm->gatherTemperature();
@@ -253,6 +260,11 @@ class SolverElastic
             // Update ghost particles.
             comm->gatherDisplacement();
 
+            // Communicate optional type.
+            if constexpr ( std::is_same<typename force_model_type::species_type,
+                                        MultiSpecies>::value )
+                comm->gatherSpecies();
+
             // Compute internal forces.
             updateForce();
 
@@ -285,6 +297,11 @@ class SolverElastic
             // Update ghost particles.
             comm->gatherDisplacement();
 
+            // Communicate optional properties.
+            if constexpr ( std::is_same<typename force_model_type::species_type,
+                                        MultiSpecies>::value )
+                comm->gatherSpecies();
+
             // Compute internal forces.
             updateForce();
 
@@ -506,7 +523,10 @@ class SolverFracture
         if ( !boundary_condition.forceUpdate() )
             boundary_condition.apply( exec_space(), *particles, 0.0 );
 
-        // Communicate temperature.
+        // Communicate optional properties.
+        if constexpr ( std::is_same<typename force_model_type::species_type,
+                                    MultiSpecies>::value )
+            comm->gatherSpecies();
         if constexpr ( std::is_same<typename force_model_type::thermal_type,
                                     TemperatureDependent>::value )
             comm->gatherTemperature();
@@ -546,6 +566,11 @@ class SolverFracture
             // Update ghost particles.
             comm->gatherDisplacement();
 
+            // Communicate optional type.
+            if constexpr ( std::is_same<typename force_model_type::species_type,
+                                        MultiSpecies>::value )
+                comm->gatherSpecies();
+
             // Compute internal forces.
             updateForce();
 
@@ -581,6 +606,10 @@ class SolverFracture
 
             // Update ghost particles.
             comm->gatherDisplacement();
+            // Communicate optional properties.
+            if constexpr ( std::is_same<typename force_model_type::species_type,
+                                        MultiSpecies>::value )
+                comm->gatherSpecies();
 
             // Compute internal forces.
             updateForce();
diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp
index 62076735..ac5d097a 100644
--- a/src/CabanaPD_Types.hpp
+++ b/src/CabanaPD_Types.hpp
@@ -20,6 +20,14 @@ struct Elastic
 struct Fracture
 {
 };
+
+struct SingleSpecies
+{
+};
+struct MultiSpecies
+{
+};
+
 struct TemperatureDependent
 {
 };
diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp
index 2a4eeeb9..b3231c40 100644
--- a/src/force/CabanaPD_ForceModels_LPS.hpp
+++ b/src/force/CabanaPD_ForceModels_LPS.hpp
@@ -23,6 +23,7 @@ template <>
 struct ForceModel<LPS, Elastic> : public BaseForceModel<Kokkos::HostSpace>
 {
     using base_type = BaseForceModel<Kokkos::HostSpace>;
+    using species_type = SingleSpecies;
     using base_model = LPS;
     using fracture_type = Elastic;
     using thermal_type = TemperatureIndependent;
@@ -60,6 +61,7 @@ template <>
 struct ForceModel<LPS, Fracture> : public ForceModel<LPS, Elastic>
 {
     using base_type = ForceModel<LPS, Elastic>;
+    using species_type = SingleSpecies;
     using base_model = typename base_type::base_model;
     using fracture_type = Fracture;
     using thermal_type = base_type::thermal_type;
@@ -81,11 +83,11 @@ struct ForceModel<LPS, Fracture> : public ForceModel<LPS, Elastic>
     {
         if ( influence_type == 1 )
         {
-            s0 = sqrt( 5.0 * G0 / 9.0 / K / delta(0) ); // 1/xi
+            s0 = sqrt( 5.0 * G0 / 9.0 / K / delta( 0 ) ); // 1/xi
         }
         else
         {
-            s0 = sqrt( 8.0 * G0 / 15.0 / K / delta(0) ); // 1
+            s0 = sqrt( 8.0 * G0 / 15.0 / K / delta( 0 ) ); // 1
         }
         bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 );
     }
@@ -95,6 +97,7 @@ template <>
 struct ForceModel<LinearLPS, Elastic> : public ForceModel<LPS, Elastic>
 {
     using base_type = ForceModel<LPS, Elastic>;
+    using species_type = SingleSpecies;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
     using thermal_type = base_type::thermal_type;
@@ -113,6 +116,7 @@ template <>
 struct ForceModel<LinearLPS, Fracture> : public ForceModel<LPS, Fracture>
 {
     using base_type = ForceModel<LPS, Fracture>;
+    using species_type = SingleSpecies;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;
     using thermal_type = base_type::thermal_type;
diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp
index b2d0709a..57c04d6f 100644
--- a/src/force/CabanaPD_ForceModels_PMB.hpp
+++ b/src/force/CabanaPD_ForceModels_PMB.hpp
@@ -23,6 +23,7 @@ struct ForceModel<PMB, Elastic, TemperatureIndependent>
     : public BaseForceModel<>
 {
     using base_type = BaseForceModel<>;
+    using species_type = typename base_type::species_type;
     using base_model = PMB;
     using fracture_type = Elastic;
     using thermal_type = TemperatureIndependent;
@@ -62,6 +63,7 @@ struct ForceModel<PMB, Elastic, TemperatureIndependent, ParticleType>
 {
     using memory_space = typename ParticleType::memory_space;
     using base_type = BaseForceModel<memory_space>;
+    using species_type = typename base_type::species_type;
     using base_model = PMB;
     using fracture_type = Elastic;
     using thermal_type = TemperatureIndependent;
@@ -121,6 +123,8 @@ struct ForceModel<PMB, Elastic, TemperatureIndependent, ParticleType>
     {
         return c( type( i ), type( j ) );
     }
+
+    void update( const ParticleType _type ) { type = _type; }
 };
 
 template <typename ParticleType>
@@ -181,6 +185,7 @@ struct ForceModel<PMB, Fracture, TemperatureIndependent, ParticleType>
         ForceModel<PMB, Elastic, TemperatureIndependent, ParticleType>;
     using memory_space = typename base_type::memory_space;
     using base_model = typename base_type::base_model;
+    using species_type = typename base_type::species_type;
     using fracture_type = Fracture;
     using thermal_type = typename base_type::thermal_type;
 
@@ -239,18 +244,18 @@ struct ForceModel<PMB, Fracture, TemperatureIndependent, ParticleType>
     }
 
     KOKKOS_INLINE_FUNCTION
-    auto criticalStretch( const int i )
+    auto criticalStretch( const int type_i )
     {
-        return sqrt( 5.0 * G0( i ) / 9.0 / K( i ) / delta( i ) );
+        return sqrt( 5.0 * G0( type_i ) / 9.0 / K( type_i ) / delta( type_i ) );
     }
 
     KOKKOS_INLINE_FUNCTION
-    auto criticalStretch( const int i, const int j ) const
+    auto criticalStretch( const int type_i, const int type_j ) const
     {
-        auto s0_i = s0( type( i ), type( i ) );
-        auto s0_j = s0( type( j ), type( j ) );
-        auto c_i = c( type( i ), type( i ) );
-        auto c_j = c( type( j ), type( j ) );
+        auto s0_i = s0( type_i, type_i );
+        auto s0_j = s0( type_j, type_j );
+        auto c_i = c( type_i, type_i );
+        auto c_j = c( type_j, type_j );
         return Kokkos::sqrt( ( s0_i * s0_i * c_i + s0_j * s0_j * c_j ) /
                              ( c_i + c_j ) );
     }
@@ -259,6 +264,10 @@ struct ForceModel<PMB, Fracture, TemperatureIndependent, ParticleType>
     bool criticalStretch( const int i, const int j, const double r,
                           const double xi ) const
     {
+        // if ( type( j ) != 0 )
+        //     std::cout << bond_break_coeff.size() << " " << type( i ) << " "
+        //               << type( j ) << "\n";
+        //<< type.size() << " " << i << " " << j << "\n";
         return r * r >= bond_break_coeff( type( i ), type( j ) ) * xi * xi;
     }
 };
@@ -281,6 +290,7 @@ struct ForceModel<LinearPMB, Elastic, TemperatureIndependent, ModelParams...>
     using base_type =
         ForceModel<PMB, Elastic, TemperatureIndependent, ModelParams...>;
     using base_model = typename base_type::base_model;
+    using species_type = typename base_type::species_type;
     using fracture_type = typename base_type::fracture_type;
     using thermal_type = typename base_type::thermal_type;
 
@@ -298,6 +308,7 @@ struct ForceModel<LinearPMB, Fracture, TemperatureIndependent, ModelParams...>
     using base_type =
         ForceModel<PMB, Fracture, TemperatureIndependent, ModelParams...>;
     using base_model = typename base_type::base_model;
+    using species_type = typename base_type::species_type;
     using fracture_type = typename base_type::fracture_type;
     using thermal_type = typename base_type::thermal_type;
 
@@ -319,6 +330,7 @@ struct ForceModel<PMB, Elastic, TemperatureDependent, TemperatureType>
 {
     using memory_space = typename TemperatureType::memory_space;
     using base_type = ForceModel<PMB, Elastic, TemperatureIndependent>;
+    using species_type = typename base_type::species_type;
     using base_temperature_type = BaseTemperatureModel<TemperatureType>;
     using base_model = PMB;
     using fracture_type = Elastic;
@@ -363,6 +375,7 @@ struct ForceModel<PMB, Fracture, TemperatureDependent, TemperatureType>
 {
     using memory_space = typename TemperatureType::memory_space;
     using base_type = ForceModel<PMB, Fracture, TemperatureIndependent>;
+    using species_type = typename base_type::species_type;
     using base_temperature_type = BaseTemperatureModel<TemperatureType>;
     using base_model = typename base_type::base_model;
     using fracture_type = typename base_type::fracture_type;