While fast folding of small proteins has been relatively well characterized by experiments and theories, much less is known for slow folding of larger proteins, for which recent experiments suggested quite complex and rich folding behaviors. Here, we address how the energy landscape theory can be applied to these slow folding reactions. Combining the perfect-funnel approximation with a multiscale method, we first extended our previous atomic-interaction based coarse grained (AICG) model to take into account local flexibility of protein molecules. Using this model, we then investigated the energy landscapes and folding routes of two proteins with complex topologies: a multidomain protein adenylate kinase (AKE) and a knotted protein 2ouf-knot. In the AKE folding, consistent with experimental results, the kinetic free energy surface showed several substates between the fully unfolded and native states. We characterized the structural features of these substates and transitions among them, finding temperature-dependent multiroute folding. For protein 2ouf-knot, we found that the improved atomic-interaction based coarse-grained model can spontaneously tie a knot and fold the protein with a probability up to 96%. The computed folding rate of the knotted protein was much slower than that of its unknotted counterpart, in agreement with experimental findings. Similar to the AKE case, the 2ouf-knot folding exhibited several substates and transitions among them. Interestingly, we found a dead-end substate that lacks the knot, thus suggesting backtracking mechanisms.