We improve the perturbative QCD (PQCD) formalism for multi-body charmless hadronic B meson decays, such as (Brightarrow VP_3 rightarrow P_1P_2P_3), by resolving the possible discrepancy in parametrizing the contribution of the P-wave resonance V to the two-meson distribution amplitudes (DAs) associated with the pairs (P_1P_2=pi pi , Kpi , KK). The determination of the Gegenbauer moments in the two-meson DAs is then updated in the global fit of the improved PQCD factorization formulas at leading order in the strong coupling (alpha _s) to available data for branching ratios and polarization fractions of three- and four-body B decays. The convergence of the Gegenbauer expansion of the resultant two-meson DAs is manifest. The satisfactory quality of the fit implies the consistency of the PQCD framework for multi-body B decays and the universality of the nonperturbative two-meson DAs. In particular, the predicted longitudinal polarization fraction (f_0(B_s^0rightarrow K^{*0} {bar{K}}^{*0})=28.2^{+8.8}_{-9.5} %) with the updated Gegenbauer moments matches well the measurement. The observable (L_{K^{*0}{bar{K}}^{*0}}=7.7^{+4.9}_{-3.8}), defined as the ratio of the longitudinal amplitudes of the two U-spin related channels (B_s^0rightarrow K^{*0} {bar{K}}^{*0}) and (B^0rightarrow K^{*0} {bar{K}}^{*0}), accommodates the current data within errors. It is found that the direct CP asymmetries ({{mathcal {A}}}^{0,||,bot }_{textrm{CP}}) in the polarization states of some four-body decays (Brightarrow V_1V_2rightarrow (P_1P_2)(P_3P_4)) might be large, but the destruction among them result in small net CP violation. Our predictions can be confronted with LHCb and Belle-II data in the future.