Traditional molecular dynamics (MD) shock simulations to investigate spallation characteristics of materials are computationally expensive as they require a compressive pulse to be generated and propagated in a domain, followed by its reflection to a tensile pulse. We propose an efficient MD-based method to cause spallation by propagating two tensile shock fronts of equal amplitude from opposite ends of a domain and meet at the center where the stress doubles and causes deformation and fracture. We investigated the spallation characteristics of monocrystalline boron carbide along the [111] orientation (along the three-atom chain) and at 90° to [111]. In both cases, boron carbide exhibited a linear-elastic axial stress-strain response up to the yield stress followed by non-linear behavior that differed significantly in the two crystal orientations. Along [111], post-yield softening was noted followed by an almost perfectly plastic response. In contrast, a post-yield hardening was observed at 90° to [111] up to a peak stress followed by an abrupt loss of strength. Further, spallation along [111] was accompanied by micro-crack initiation normal to the loading direction while at 90° to [111], spallation was preceded by crack formation at ~45° to the loading direction. To explain these observations, a systematic analysis of bond deformation was conducted, wherein the average strain experienced by each of the 41 bonds in a unit cell of boron carbide was estimated. It was noted that the onset of plastic deformation along both crystal orientations is triggered by breakage of B-B inter-icosahedral bonds followed by B-C icosahedral chain bonds. The sequence in which the bonds break and the fraction of total number of broken bonds result in different macroscopic (domain-level) responses. This enabled identification of bonds which experience the highest strains (and damage), thus explaining for the observed differences.