The current analysis aims to present the results of a QCD analysis of diffractive parton distribution functions (PDFs) at next-to-leading-order accuracy in perturbative QCD. In this new determination of diffractive PDFs, we use all available and up-to-date diffractive deep-inelastic scattering (DIS) datasets from H1 and ZEUS collaborations at HERA, including the most recent H1/ZEUS combined measurements. In this analysis, we consider the heavy quark contributions to the diffractive DIS in the so-called framework of the fonll general mass variable flavor number scheme. The uncertainties on the diffractive PDFs are calculated using the standard “Hessian error propagation,” which served to provide a more realistic estimate of the uncertainties. This analysis is enriched, for the first time, by including the nonperturbative higher twist (HT) effects in the calculation of diffractive DIS cross sections, which are particularly important at large −x and low Q2 regions. Then, the stability and reliability of the extracted diffractive PDFs are investigated upon the inclusion of HT effects. We discuss the novel aspects of the approach used in this QCD fit, namely, optimized and flexible parametrizations of diffractive PDFs, the inclusion of HT effects, and considering the recent H1/ZEUS combined dataset. Finally, we present the extracted diffractive PDFs with and without the presence of HT effects and discuss the fit quality and the stability upon variations of the kinematic cuts and the fitted datasets. We show that the inclusion of HT effects in diffractive DIS can improve the description of the data, which leads, in general, to a very good agreement between data and theory predictions.