      subroutine bfgs_r(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

      den=zero

      do 1004 i=1,n
 1004 den=den+delta_x(i)*delta_g(i)

      do 1005 j=1,n
      do 1006 i=1,n
 1006 wkmat1(i,j)=   -(delta_x(i)*delta_g(j))/den
 1005 wkmat1(j,j)=one-(delta_x(j)*delta_g(j))/den

      do 1007 j=1,n
      do 1008 i=1,n
 1008 wkmat2(i,j)=zero
      do 1007 k=1,n
      do 1007 i=1,n
 1007 wkmat2(i,j)=wkmat2(i,j)+hessian(i,k)*wkmat1(j,k)

      do 1009 j=1,n
      do 1010 i=1,n
 1010 hessian(i,j)=zero
      do 1009 k=1,n
      do 1009 i=1,n
 1009 hessian(i,j)=hessian(i,j)+wkmat1(i,k)*wkmat2(k,j)

      do 1011 j=1,n
      do 1011 i=1,n
 1011 hessian(i,j)=hessian(i,j)+(delta_x(i)*delta_x(j))/den

      return
      end
