program bee1
	dimension atom(50),aa(40,40),ss(40,40)
	dimension title(20)
	common/jac01/n,a(40,40),s(40,40)
    2 format(i5,i5,f10.4,25x,25x,i2)
    3 format(1h0,4x,i5,i5,e20.8)
    5 format(1h0,2x,10heigenvalue,i3,3h = ,e15.8,5h beta)
    6 format(1h0,2x,30hcoefficients of wave function , i3)
    7 format(18a4)
    8 format(1h0)
    9 format(1h ,e12.5)
   66 format(1h0,20f6.2)
   77 format(1h1,18a4)
  199 format(1h ,'pi-bond order',2i2,' is ',e12.5,4x,'electron
     1density atom',i2,4h is ,e12.5)
  204 format(1h ,38x,21helectron density atom,i2,4h is ,e12.5)
  270 format(1h ,'number of atoms in this molecule   
     1is',2x,i3,2x,' and the number of electrons is',2x,i3)
  299 format(30i2)
  600 format(56heigenvalues listed in decreasing energy,in units 
     1of beta)
      open(1,file='bee1.dat',status='old')
  135 read(1,7)(title(i),i=1,18)
  170 continue
        read(1,*)n
	read(1,*)ne
	read(1,*)(atom(i),i=1,n)
	do 150 i=1,n
	do 150 j=1,n
	a(i,j)=0.0
	if(i-j)100,101,100 
  100 s(i,j)=0.0
	go to 150
  101 s(i,j)=1.0
  150 continue
  152 read(1,*)i,j,fa,ind
	a(i,j)=fa
	a(j,i)=fa
	if(ind)151,152,151
  151 continue
	open(9,file='bee1.out')
	write(9,77)(title(i),i=1,18)
	write(6,77)(title(i),i=1,18)
	write(9,270)n,ne
	write(6,270)n,ne
	do 500 i=1,n
  500 write(9,66)(a(i,j),j=1,n)
      write(6,66)(a(i,j),j=1,n)
	call jacb
	i=1
	amin=a(1,1)
	in=2
  300 do 301 id=in,n
	if(a(id,id)-amin)302,301,301
  302 amax=amin
	amin=a(id,id)
	a(i,i)=amin
	a(id,id)=amax
	do 301 ip=1,n
	as=s(ip,i)
	bs=s(ip,id)
	s(ip,id)=as
	s(ip,i)=bs
  301 continue
	i=i+1
	if(i-n)305,304,304
  305 amin=a(in,in)
	in=in+1
	go to 300
  304 continue
	do 130 i=1,n
	do 130 j=1,n
	if(i-j)130,131,130
  131	write(9,600)
   	write(6,600)
        write(9,5)i,a(i,j)
        write(6,5)i,a(i,j)
	write(9,6)j
	do 208 ip=1,n
  208 write(9,3)ip,n-j+1,s(ip,n-j+1)
  130 continue
	ii=1
  201 i=atom(ii)
	if(n-ii)203,203,202
  202 k=atom(ii+1)
  203 d=0
	pbo=0
	do 99 j=1,n
	if(ne-2*j)260,261,261
  261 el=2.0
	go to 265
  260 if(ne+1-2*j)266,262,262
  262 el=1.0
	go to 265
  266 el=0.0
  265 continue
	d=d+s(i,n-j+1)**2*el
	if(ii-n)599,99,99
  599 pbo=pbo+s(i,n-j+1)*s(k,n-j+1)*el
   99 continue
	if(ii-n)799,699,699
  799 write(9,199)i,k,pbo,i,d
	go to 111
  699 write(9,204)i,d
  111 ii=ii+1
	if(ii-n)201,201,200
  200 continue
	write(9,8)
	write(9,8)
	if(ind-99)211,212,211
  212 go to 135
  211 stop
	end

	subroutine jacb
	common/jac01/n,a(40,40),s(40,40)
	indic=0
  151  vi=0.0
	do 106 i=1,n
	do 106 j=1,n
	if(i-j)107,206,107
  107  vi=vi+a(i,j)**2
	s(i,j)=0
	go to 106
  206  s(i,j)=1.0
  106  continue
	vi=sqrt(vi)
	vf=vi*0.1e-07
	an=n
  128  vi=vi/an
  137  iq=1
  124  iq=iq+1
	ip=0
  121  ip=ip+1
	if(a(ip,iq))108,120,109
  108  if(-a(ip,iq)-vi)120,112,112
  109  if(a(ip,iq)-vi)120,112,112
  112  indic=1
	alam=-a(ip,iq)
	amu=0.5*(a(ip,ip)-a(iq,iq))
	if(amu)113,114,114
  113  sgn=-1.0
	go to 115
  114  sgn=+1.0
  115 omega=sgn*alam/sqrt(alam**2+amu**2)
	stht=omega/sqrt(2.0+2.0*sqrt(1.0-omega**2))
	ctht=sqrt(1.0-stht**2)
	do 116 i=1,n
	if(i-ip)117,118,117
  117  if(i-iq)119,118,119
  119  aipi=a(ip,i)*ctht-a(iq,i)*stht                    
	aiqi=a(ip,i)*stht+a(iq,i)*ctht
	a(ip,i)=aipi
	a(iq,i)=aiqi
  118  aipi=s(i,ip)*ctht-s(i,iq)*stht
	aiqi=s(i,ip)*stht+s(i,iq)*ctht
	s(i,ip)=aipi
  116  s(i,iq)=aiqi
	aipi=a(ip,ip)*ctht**2+a(iq,iq)*stht**2-2.*a(ip,iq)*stht*ctht
	aiqi=a(ip,ip)*stht**2+a(iq,iq)*ctht**2
     1        +2.0*a(ip,iq)*stht*ctht
	aipiq=(a(ip,ip)-a(iq,iq))*ctht*stht
     1        +a(ip,iq)*(ctht**2-stht**2)
	a(ip,ip)=aipi
	a(iq,iq)=aiqi
	a(ip,iq)=aipiq
	a(iq,ip)=a(ip,iq)
	do 123 i=1,n
	a(i,ip)=a(ip,i)
  123  a(i,iq)=a(iq,i)
  120  if(ip-iq+1)121,122,122
  122  if(iq-n)124,125,125
  125  if(indic)126,127,126
  126  indic=0
	go to 137
  127  if(vi-vf)129,129,128
  129  return
	end