Those three contributions are fitted independently of each other using well-defined grids generated in terms of displacement with respect to a reference geometry, Q re f . A least-square procedure is used with the ∆V ab K (Q) values on the grid points. Those values are related to the ab initio energy values, V ab (Q). The full 3D-grid is split in four parts, G 0 , G 1 , G 2 and G 3 depending on the coupling between modes. The first grid, G 0 , contains a single point, the reference geometry and the corresponding ab initio energy, E 0 . The other points are generated with constraints, R G K , on the indices [u 1 , u 2 , u 3 ]: u 1 + u 2 + u 3 ≤ L G K , where L G K control the size of the grid. The constraints for K = 1, 2, 3 are the following: (i) One mode expansion (K=1), R G 1 : The possible values of the indices are [u 1 , 0, 0] or [0, u 2 , 0] or [0, 0, u 3 ] (with u 1 ≥ 1, u 2 ≥ 1 and u 3 ≥ 1). With L G 1 = 5, the number of grid points is 25 (10 for [u 1 , 0, 0] or [0, u 2 , 0] and only 5 for [0, 0, u 3 ]. Then, the ∆V ab 1 (Q) values are obtained as follows: ∆V ab 1 (Q) = V ab (Q) -E 0 . (ii) Two mode expansion (K=2), R G 2 : The possible values of the indices are [u 1 , u 2 , 0] or [u 1 , 0, u 3 ] or [0, u 2 , u 3 ] (with u 1 ≥ 1, u 2 ≥ 1 and u 3 ≥ 2). With L G 2 = 5, the number of grid points is 80 (40 for [u 1 , y 2 , 0] and 20 [u 1 , 0, u 3 ] or [0, u 2 , u 3 ]. The ∆V ab 2 (Q) is given by: