We adopt the neural network (NN) flow method to study the Berezinskii-Kosterlitz-Thouless (BKT) phase transitions of the two-dimensional q-state clock model with q≥4. The NN flow consists of a sequence of the same units that proceed with the flow. This unit is a variational autoencoder trained by the data of Monte Carlo configurations in unsupervised learning. To gauge the difference among the ensembles of Monte Carlo configurations at different temperatures and the uniqueness of the ensemble of NN-flow states, we adopt the Jensen-Shannon divergence (JSD) as the information-distance measure "thermometer." This JSD thermometer compares the probability distribution functions of the mean spin value of two ensembles of states. Our results show that the NN flow will flow an arbitrary spin state to some state in a fixed-point ensemble of states. The corresponding JSD of the fixed-point ensemble takes a unique profile with peculiar features, which can help to identify the critical temperatures of BKT phase transitions of the underlying Monte Carlo configurations.