Significant computing costs can be saved by first constructing a high-quality PES based on accurate but expensive ab initio points covering all dynamically relevant regions. However, it is not an easy task to obtain reliable full-dimensional PESs, whose topology is complicated with reactants and products, saddle points and intermediates. As discussed above, a set of molecular configurations need be sampled first and their electronic energies (and / or forces) are computed at an appropriate level of theory. For a molecular system with Natom atoms (Natom ≥ 3), the full-dimensional PES is dependent on 3Natom -6 molecular intramolecular coordinates. The sampling size increases remarkably with the size of the system. Thanks to the huge development in electronic structure theory and the rapid progress in computing power, it is now straightforward to calculate reliable energies of hundreds of thousands of configurations. There are many electronic structure methods of varying accuracies and computing speeds, such as the Hartree-Fock method, the Moller-Plesset (MPn) perturbation methods, and more sophisticated electron correlation methods, including the highly accurate coupled cluster method (e.g. CCSD(T)), the configuration interaction method, 38 the complete-active-space self-consistent field (CASSCF) or its variants, CASPT2, CASPT3, and the multi-reference configuration interaction method (MRCI). Generally, with a bigger basis set and more sophisticated post-Hartree-Fock electron correlation treatment, the result is generally more reliable. Recently, Dawes and his coworkers have reviewed recent developments in single-and multireference electronic structure methods that are used to obtain accurate energies for constructing global molecular PESs. Nonetheless, one has to perform comprehensive tests to choose an appropriate electronic structure method for good and consistent description of all sampling data with acceptable calculation cost. The next step is to sample sufficiently a large number of points in the configuration space. The sampling data points should have sufficient wide ranges with appropriate sampling densities in both the configuration and energy spaces for specified research targets. This is because that with any method to represent the PES, the performance of the PES can't be guaranteed for regions where no data is sampled. Generating configurations on a random or regular grid is an efficient and convenient way for small molecules (under four atoms) with only a few degrees of freedom. However, this sampling approach becomes unfeasible for larger systems due to the exponential growth of the number of data points with the system size. For a specified system, we need to consider only regions where the quantum/classical reaction dynamics simulations would access. Running classical dynamical simulations at some low-level methods, such as density functional theory (DFT), is efficient to sample points in the dynamically related regions. Note that the sampling is related to a specified research target. For instance, some regions are classically forbidden, but maybe accessed in QD simulations. Various techniques, such as active learning and concurrent learning, were developed to select the most and probably the smallest representative samples effectively. Then the sampling data points, which are represented as a set of coordinates and its associated reliable electronic energies, are used to obtain an analytical function of the nuclear degrees of freedom, namely, the PES, whose functional form can't be derived generally based on physical considerations. Many forms of analytical functions have been suggested to represent the PES and summarized in the review by Jiang et al. Among them, the permutation invariant polynomial (PIP) method is powerful to represent the full-dimensional PESs with a reasonable accuracy and exact invariant property with respect to the permutation of like atoms. Neural network (NN) approach, a nonlinear regression method, is a universal approximator as it can in principle approximate any unknown multidimensional function with arbitrary accuracy by learning from a set of known data, namely, the sampling data points. Many NN methods have demonstrated excellent fitting performance in representing PESs of various systems in the gas phase, and large condensed systems from bulk solids via surfaces and interfaces to solutions. Please refer to these reviews and literatures cites therein for further detailed descriptions and discussions.