  The behavior of polyatomic molecules around their equilibrium positions can
be regarded as quantum coupled anharmonic oscillators. Solving the
corresponding Schr\"odinger equations can interpret or predict experimental
spectra of molecules. In this study, we develop a novel approach to solve
excited states of anharmonic vibrational systems. The normal coordinates of
molecules are transformed into new coordinates through a normalizing flow
parameterized by a neural network, facilitating the construction of a set of
orthogonal many-body variational wavefunctions. Our methodology has been
validated on an exactly solvable $64$-dimensional coupled harmonic oscillator,
yielding numerical results with a relative error on the order of $10^{-6}$.
Furthermore, the neural canonical transformations are also applied to calculate
the energy levels of two specific molecules, acetonitrile
($\text{C}\text{H}_3\text{CN}$) and ethylene oxide
($\text{C}_2\text{H}_4\text{O}$) involving $12$ and $15$ vibrational modes
respectively, which are in agreement with experimental findings. One of the key
advantages of our approach is its flexibility concerning the potential energy
surface, requiring no specific form for its application. Moreover, this method
can be easily implemented on large-scale distributed computing platforms,
making it particularly suitable for investigating more complex vibrational
structures.
