Locating the low free energy paths (LFEPs) connecting different conformational states is among the major tasks for the simulations of complex biomolecules as the pathways encode the physical essence and, therefore, the underlying mechanism for their functional dynamics. Finding the LFEPs is yet challenging due to the numerous degrees of freedom of the molecules and expensive force calculations. To alleviate this issue, we have previously introduced a Traveling-Salesman-based Automated Path Searching (TAPS) approach that requires minimal input information to locate the LFEP closest to a given initial guess path. Despite its high efficiency for large biomolecules, it remains, as all path-searching methods, incapable of revealing multiple parallel LFEPs simultaneously, which are, however, near-ubiquitous. This work describes a comprehensive protocol that offers parallel LFEPs efficiently. Our protocol starts with a modified version of the parallel cascade approach, which extensively searches for a large pile of geometrically distinct paths of the target molecule in implicit solvents. These paths are clustered and then filtered by their cumulative barriers, yielding a smaller set of initial paths for subsequent optimization by TAPS in explicit solvents. Through this protocol, we successfully sampled eight LFEPs for the transition of Met-enkephalin from its 3