Dear Mohammed Gharaibeh

It is possible to calculate BO empirically from rho at BCP, provided two participating nuclei are same. For this you should have few data of rho at BCP for pure single, double and triple bond. The you can fit the the data empirically to the expression

n=exp[A(rho-B)]; n= bond order and A and B are constant

Then with the your value of rho(BCP), you can guess the BO of the bond you want.

For more detail please read Chapter 10 in book: Atom in Molecule, An introduction by Prof. P. Popelier

Regards

Piyush

The known ways to estimate the bond order in the framework of  Bader theory are too semiempirical and not reliable. It is better to calculate this quantity using electronic density matrix in the basis of occupied molecular orbitals and overlap integrals (say, in terms of the Wiberg or Mayer indices).

A recent paper introduces a new method for computing bond orders that is more general and more robust than previous approaches: T. A. Manz, "Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders," RSC Advances, 7 (2017) 45552-45581 (http://dx.doi.org/10.1039/C7RA07400J ). The method introduced in this paper computes the bond order as a functional of the electron density (and the spin density for magnetic materials) on a set of grid points. Therefore, it should accomplish what you seek. You do not need to know the density matrix or orbitals to compute these bond orders.

The method is implemented in the free Chargemol program (http://ddec.sourceforge.net). It will compute and print the net atomic charges, bond orders, and other information. Chargemol can analyze the result of GAUSSIAN, VASP, and other quantum chemistry programs.