The depth of focus of a microscopy imaging system can often fall below the sample thickness when high spatial resolution is demanded. It is known that multiple scatterings of waves in the sample volume is responsible for the constrained depth of focus, but the inversion of this effect is difficult because it relies on the knowledge of the content of the sample, which is exactly the unknown to be solved. Therefore, we have developed an iterative algorithm to address this issue. We initialize the sample function with a random guess, and propagate the incident beam through the heterogeneous sample using multislice propagation. By updating the sample function to minimize the discrepancy between the calculated intensity and the measured data, we may finally converge onto the reconstructed object. the optimization of the loss function is conducted using automatic differentiation, which saves us from the tedious derivations of the gradients, and makes the implementation easily adaptable for the different imaging techniques. Our codes have been deployed on ALCF Cooley and Theta, where the distribution of tasks on multiple nodes allows the computationally intensive reconstruction to complete with a reasonable amount of time.