The direct estimation of thermal conductivity in solids using the constant flux non-equilibrium molecular dynamics (NEMD) simulation method has proven to be difficult with significant differences still existing for various, even simple, systems. The application of NEMD to study the thermal transport in molecular solids was rarely attempted. We report the first application of NEMD based on the Ikeshoji and Hafskjold (NEMD-IH) scheme to study the thermal conductivity of crystalline hexahydro-1,3,5-trinitro-s-triazine (RDX), an important energetic material, at temperatures T = 230, 275, and 300 K for different sample sizes and different crystallographic directions. We found that the size effects arising from a finite distance between sources of heat in the NEMD-IH simulation show a well-defined linear dependence on the inverse of the sample length within the errors of simulation. The thermal conductivity exhibits anisotropy with the largest value found along the [0 0 1] crystallographic direction and the smallest value along the [0 1 0] direction. The extrapolated value of the directionally averaged thermal conductivity coefficient (kappa) over bar to the sample of infinite length (e.g. 0.355 W/Km at T = 300 K) is higher, but still in reasonable agreement with existing experimental data for powdered samples. The NEMD-IH heat transport properties were reproduced in a broad range of values of heat flux and sizes of heat source regions. Surprisingly, the thermal conductivity coefficient at T = 275 K is found to be higher by about 13% than at T = 300 K and 230 K, which may suggest a location of the temperature maximum above 230 K. A similar location of the low temperature peak was reported from experiments on powdered samples of RDX. Through the use of additional calculations based on the Green-Kubo approach, we discuss the possible origins of such a behavior in the NEMD-IH approach. Published by Elsevier Ltd.