The intracellular polymerization of deoxy-sickle cell hemoglobin (HbS) has been identified as the main cause of sickle cell disease. Therefore, the material properties and biomechanical behavior of polymerized HbS fibers is a topic of intense research interest. A solvent-free coarse-grain molecular dynamics (CGMD) model is developed to represent a single hemoglobin fiber as four tightly bonded chains. A finitely extensible nonlinear elastic (FENE) potential, a bending potential, a torsional potential, a truncated Lennard-Jones potential and a Lennard-Jones potential are implemented along with the Langevin thermostat to simulate the behavior of a polymerized HbS fiber in the cytoplasm. The parameters of the potentials are identified via comparison of the simulation results to the experimentally measured values of bending and torsional rigidity of single HbS fibers. After it is shown that the proposed model is able to very efficiently simulate the mechanical behavior of single HbS fibers, it is employed in the study of the interaction between HbS fibers. It is illustrated that frustrated fibers and fibers under compression require a much larger interaction force to zipper than free fibers resulting in partial unzippering of these fibers. Continuous polymerization of the unzippered fibers via heterogeneous nucleation and additional unzippering under compression can explain the formation of HbS fiber networks and consequently the wide variety of shapes of deoxygenated sickle cells.