Realism of models has two components: sophistication of implemented physics and realism of the overall setup. Current 3D radiation MHD simulations of the coupled solar atmosphere are presenting a compromise on a middle ground. Domains are at best large enough to capture the extent of small active regions, but the commonly adapted use of periodic boundary conditions leads in many setups to an unrealistic magnetic field connectivity. Relevant physics in terms of radiation transport, including non-equilibrium treatment and ion-neutral effects are implement in a few commonly used simulation codes, but in many of larger domain models the resolution is just barely enough to start resolving processes of interest. More detailed physics better describing reconnection and energy dissipation can be only studied in dedicated simulations that focus on either small domains or reduce the dimensionality of the domain. In this talk I will touch on these challenges and present a few future directions: (1) We need simulations in local domains that are better informed by the global magnetic field structure of the corona. (2) We need higher resolution and reduction of numerical diffusivity and adoption of sub-grid models that better capture underlying physics. (3) We need a better quantification of model errors. Any combination of these improvements will come with a high cost and will require utilization of the latest computing platforms, including the use of GPUs. Data analysis is developing into a major bottle neck and will require a change in how we conduct and share large numerical simulations.