We construct theoretically "exact" and numerically "accurate" Beyond Born-Oppenheimer (BBO) based diabatic potential energy surfaces (PESs) of pyrazine (C 4 N 2 H 4 ) molecule involving lowest four excited adiabatic PESs (S 1 to S 4 ) and nonadiabatic coupling terms (NACTs) among those surfaces as functions of nonadiabatically active normal modes (Q 1 , Q 6a , Q 9a and Q 10a ) to compute its photoabsorption (PA) spectra. Those adiabatic PESs are calculated using CASSCF as well as MRCI based methodologies, where NACTs are obtained from CP-MCSCF approach. Employing ab initio quantities (adiabatic PESs and NACTs), it is possible to depict the conical intersections (CIs) and develop matrices of diabatic PESs over six normal mode planes. Once single-valued, smooth, symmetric and continuous 2 × 2 and 4 × 4 diabatic surface matrices are in hand for the first time, such matrices are used to perform multi-state multi-mode nuclear dynamics with the aid of Time-Dependent Discrete Variable Representation (TDDVR) methodology initializing the product type wavefunction on 1 1 B 1u (S 1 ) and 1 1 B 2u (S 2 ) states to obtain the corresponding PA spectra. TDDVR calculated spectra for those states (S 1 and S 2 ) obtained from BBO based 2 × 2 and 4 × 4 diabatic surface matrices show good and better agreement with the experimental results, respectively. Both of these calculated results depict better peak progression over the existing profiles of Multi-Configuration Time-Dependent Hartree (MCTDH) dynamics over 2 × 2 Vibronic Coupling Model (VCM) Hamiltonian.