  Significant differences exist among literature for thermal conductivity of
various systems computed using molecular dynamics simulation. In some cases,
unphysical results, for example, negative thermal conductivity, have been
found. Using GaN as an example case and the direct non-equilibrium method,
extensive molecular dynamics simulations and Monte Carlo analysis of the
results have been carried out to quantify the uncertainty level of the
molecular dynamics methods and to identify the conditions that can yield
sufficiently accurate calculations of thermal conductivity. We found that the
errors of the calculations are mainly due to the statistical thermal
fluctuations. Extrapolating results to the limit of an infinite-size system
tend to magnify the errors and occasionally lead to unphysical results. The
error in bulk estimates can be reduced by performing longer time averages using
properly selected systems over a range of sample lengths. If the errors in the
conductivity estimates associated with each of the sample lengths are kept
below a certain threshold, the likelihood of obtaining unphysical bulk values
becomes insignificant. Using a Monte-Carlo approach developed here, we have
determined the probability distributions for the bulk thermal conductivities
obtained using the direct method. We also have observed a nonlinear effect that
can become a source of significant errors. For the extremely accurate results
presented here, we predict a [0001] GaN thermal conductivity of 185 $\rm{W/K
\cdot m}$ at 300 K, 102 $\rm{W/K \cdot m}$ at 500 K, and 74 $\rm{W/K \cdot m}$
at 800 K. Using the insights obtained in the work, we have achieved a
corresponding error level (standard deviation) for the bulk (infinite sample
length) GaN thermal conductivity of less than 10 $\rm{W/K \cdot m}$, 5 $\rm{W/K
\cdot m}$, and 15 $\rm{W/K \cdot m}$ at 300 K, 500 K, and 800 K respectively.
