How can I determine the vector parallel to the long molecular axis?

You may want to look at the principal axes of the molecule, which can be calculated from the inertia matrix in the body frame. The theory is discussed in rigid body dynamics and you should be able to find examples of the calculation on the internet. The calculation procedure can be outlined as follows:

First, calculate the center of mass of the rigid molecule. Find the x,y,z coordinates of the atoms of your molecule in the body frame, a reference frame whose origin is at the center of mass of the molecule.

Second, use the body frame coordinates to calculate the inertia matrix in the body frame. The inertia matrix must look like this: $$begin{pmatrix}I_{xx} & -I_{xy} & -I_{xz} -I_{yx} & I_{yy} & -I_{yz} -I_{zx} & -I_{zy} & I_{zz} end{pmatrix}$$

I put the minus signs on the products of inertia because they can be easily forgotten. Sometimes the minus signs are included in the definition of $I_{jk}$, in which case they should not appear again in the inertia matrix.

Third, find the eigenvalues and eigenvectors of the inertia matrix. The three eigenvectors will be parallel to the principal axes of the rigid molecule. The axes themselves are lines that pass through the center of mass of the rigid molecule. In your example molecule, the eigenvector that corresponds to the largest eigenvalue should be nearly parallel to the line through the blue nitrogen atoms.

Fourth, perform a sanity check on your calculations. This is optional, but recommended. Normalize the three eigenvectors so they have unit length. Use the normalized eigenvectors as columns of a matrix, T. Multiply the x,y,z body frame coordinates of each atom by matrix T to get new coordinates. Use these new coordinates to calculate a new inertia matrix. If all has gone well, the new inertia matrix will be diagonal (or have very small off diagonal values) and the values on its diagonal will be the same as the eigenvalues calculated for the first inertia matrix.