||The thermal conductivity of graphite nano-fibers (GNFs) with different styles is predicted computationally. GNFs are formed as basal planes of graphene stacked based on the catalytic configuration. The large GNF thermal conductivity relative to a base phase change material (PCM) may lead to improved PCM performance when embedded with GNFs. Three different types of GNFs are modeled: platelet, ribbon, and herringbone. Molecular dynamics (MD) simulations are used in this study as a means to predict the thermal conductivity tensor based on atomic behavior. The in-house MD code, Molecular Dynamics in Arbitrary Geometries (MDAG), was updated with the features required to create the predictions. To model both interlayer van-der Waals and intralayer covalent bonding of carbon atoms in GNFs, a combination of the optimized Tersoff potential function for atoms within the layers and a pairwise Lennard-Jones (LJ) potential function to model the interactions between the layers was used. Tests of energy conservation in the NVE ensemble have been performed to validate the employed potential model. Nose-Hoover, Andersen, and Berendsen thermostats were also incorporated into MDAG to enable MD simulations in NVT ensembles, where the volume, number of atoms, and temperature of the system are conserved. Equilibrium MD with Green-Kubo (GK) relations was then employed to extract the thermal conductivity tensor for symmetric GNFs (platelet and ribbon). The thermal conductivity of solid argon at different temperatures was calculated and compared to other studies to validate the GK implementation. Different heat current formulations, as a result of using the three-body Tersoff potential, were considered and the discrepancy in the calculated thermal conductivity values of graphene using each formula was resolved by employing a novel comparative technique that identifies the most accurate formulation. The effect of stacking configuration on the thermal conductivity of platelet and ribbon GNFs was also investigated using equilibrium molecular dynamics (EMD) with GK relations. Simple Hexagonal (AAA), Bernal (ABA), and Rhombohedral (ABC) stacking forms were considered. The intralayer and interlayer thermal conductivity values were predicted in both zigzag and armchair directions to be in the range of 450-800 W/m.K and 17-55 W/m.K, respectively. Furthermore, non-equilibrium molecular dynamics (NEMD) simulations were used to investigate the thermal conductivity of herringbone graphite nanofibers (GNFs) at room temperature by breaking down the axial and transverse conductivity values into intralayer and interlayer components. The edge effect on a layer's thermal conductivity was investigated by computing the thermal conductivity values in both zigzag and armchair directions of the heat flow. The limiting case of a 90 degree crease angle was used to compare the results with those of single-layer graphene and few-layer graphene. The thermal conductivity values in the axial, transverse in the crease direction, and transverse normal to the crease directions for the case of a five-layer herringbone GNF with a 45-degree crease angle were calculated to be 27 W/m.K, 263 W/m.K, and 1500 W/m.K, respectively.