# Bond order analysis in molpro

I'm trying to compute bond orders in a molecule, let's consider S2 as an example. I'm not entirely sure what's the right way to do that, but I assume NBO (Natural Bond Orbital) should give me what I need. For calculations I'm using MOLPRO. I've tried to launch it with UCCSD(T)-F12 as follows:

 symmetry, nosym
geometry = {
S1
S2, S1, bond
}
bond = 3.60529338 BOHR
basis = cc-pVTZ-F12
set, spin = 0

hf
uccsd(t)-f12a, scale_trip = 1
nbo, keep_wbi = 1


It seemed to work and gave me the following output:

  - - - - - - - - - - - - - - - - - - - - - - - -
NBO  OCCUPATION       LABEL     CENTERS
- - - - - - - - - - - - - - - - - - - - - - - -
1    2.000000         CR       1
2    1.999827         CR       1
3    1.999918         CR       1
4    1.999989         CR       1
5    1.999848         CR       1
6    2.000000         CR       2
7    1.999827         CR       2
8    1.999918         CR       2
9    1.999989         CR       2
10    1.999848         CR       2
11    1.999151         LP       1
12    1.989845         LP       1
13    1.999151         LP       2
14    1.989845         LP       2
15    2.000000         BD       1    2
16    2.000000         BD       1    2


CR seems to stand for core, LP - lone pair, BD - bonding. From this I can conclude the bond order is 2, since we have two bonding orbitals with 2 electrons on each. However, the spin state was set to 0 here, whereas it's supposed to be 2 (2 unpaired electrons). When I tried to set spin to 2, the same calculation failed:

No suitable structure found!
? Error
? The problem occurs in nbo_transform


Can anyone suggest what is wrong here? Are there any other ways to compute bond order?

The full outputs: spin=0; spin=2

• You could try the many methods by Multiwfn for bond orders. – pH13 - Yet another Philipp Jun 12 '18 at 17:01