      subroutine bofill(icycle,n,hessian,grad_1,coord_1,grad_0,
     &                  coord_0,delta_g,delta_x,wkmat1,wkmat2)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)  

      dimension hessian(n,*),grad_1(*),coord_1(*),grad_0(*),coord_0(*),
     &          delta_g(*),delta_x(*),wkmat1(n,*),wkmat2(n,*)

      data zero,one,two,pi/0.0d0,1.0d0,2.0d0,3.1415926535898d0/

      rewind 21

      call fastrd(21,hessian,n*n)
      call fastrd(21,grad_0,n)
      call fastrd(21,coord_0,n)

      do 1001 i=1,n
 1001 delta_x(i)=coord_1(i)-coord_0(i)

      do 1002 i=1,n
 1002 delta_g(i)=grad_1(i)-grad_0(i)

      do 1003 i=1,n
      if(delta_x(i).lt.(-pi)) delta_x(i)=delta_x(i)+two*pi
 1003 if(delta_x(i).gt.(+pi)) delta_x(i)=delta_x(i)-two*pi

      do 1004 j=1,n
      do 1004 i=1,n
 1004 delta_g(i)=delta_g(i)-hessian(i,j)*delta_x(j)

      fac1=zero
      fac2=zero
      fac3=zero

      do 1005 i=1,n
 1005 fac1=fac1+delta_x(i)*delta_x(i)

      do 1006 i=1,n
 1006 fac2=fac2+delta_g(i)*delta_g(i)

      do 1007 i=1,n
 1007 fac3=fac3+delta_g(i)*delta_x(i)

      factor=fac3/fac1

      do 1008 j=1,n
      do 1008 i=1,n
 1008 wkmat1(i,j)=zero

      do 1009 j=1,n
      do 1009 i=1,n
 1009 wkmat1(i,j)=wkmat1(i,j)+delta_g(i)*delta_x(j)
     &                       +delta_x(i)*delta_g(j)
     &                       -delta_x(i)*delta_x(j)*factor

      do 1010 j=1,n
      do 1010 i=1,n
 1010 wkmat1(i,j)=wkmat1(i,j)/fac1

      do 1011 j=1,n
      do 1011 i=1,n
 1011 wkmat2(i,j)=zero

      do 1012 j=1,n
      do 1012 i=1,n
 1012 wkmat2(i,j)=wkmat2(i,j)+delta_g(i)*delta_g(j)

      do 1013 j=1,n
      do 1013 i=1,n
 1013 wkmat2(i,j)=wkmat2(i,j)/fac3

      phi=one-((fac3*fac3)/(fac1*fac2))

      do 1014 j=1,n
      do 1014 i=1,n
 1014 hessian(i,j)=hessian(i,j)+phi*wkmat1(i,j)+(one-phi)*wkmat2(i,j)

      return
      end
