      subroutine hguess(nmodes,nbonds,nangles,ndihedrals,nimpropers,
     &                  natomtype,i_bond,i_angle,i_dihedral,i_improper,
     &                  coord,hessian)

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 radii(0:54),radius(0:54)

      dimension natomtype(*),i_bond(2,*),i_angle(3,*),i_dihedral(4,*),
     &          i_improper(4,*)

      dimension coord(3,*),hessian(nmodes,*)

      data zero,bohr/0.0d0,0.529177d0/

      data radii/1.00d0,0.32d0,0.60d0,1.20d0,1.05d0,0.81d0,0.77d0,
     &                  0.74d0,0.74d0,0.72d0,0.72d0,1.50d0,1.40d0,
     &                  1.30d0,1.17d0,1.10d0,1.04d0,0.99d0,0.99d0,
     &                  1.80d0,1.60d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.30d0,1.20d0,1.20d0,1.10d0,1.10d0,
     &                  1.80d0,1.60d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,1.40d0,
     &                  1.40d0,1.30d0,1.20d0,1.20d0,1.10d0,1.10d0/

      data a_b,b_b/0.3601d0,1.944d0/
      data a_a,b_a,c_a,d_a/0.089d0,0.11d0,0.44d0,-0.42d0/
      data a_d,b_d,c_d,d_d,e_d/0.0015d0,14.0d0,2.85d0,0.57d0,4.00d0/
      data a_i,b_i,c_i,d_i,e_i/0.0025d0,0.0061d0,3.00d0,4.00d0,0.80d0/

      do 1001 i=0,54
 1001 radius(i)=radii(i)/bohr

      do 1002 j=1,nmodes
      do 1002 i=1,nmodes
 1002 hessian(i,j)=zero

      do 1003 i=1,nbonds
      r_ab=dist(coord(1,i_bond(1,i)),coord(1,i_bond(2,i)))
      r_ab_cov=radius(natomtype(i_bond(1,i)))
     &        +radius(natomtype(i_bond(2,i)))
 1003 hessian(i,i)=a_b*exp(-b_b*(r_ab-r_ab_cov))

      do 1004 i=1,nangles
      r_ab=dist(coord(1,i_angle(2,i)),coord(1,i_angle(1,i)))
      r_ac=dist(coord(1,i_angle(2,i)),coord(1,i_angle(3,i)))
      r_ab_cov=radius(natomtype(i_angle(2,i)))
     &        +radius(natomtype(i_angle(1,i)))
      r_ac_cov=radius(natomtype(i_angle(2,i)))
     &        +radius(natomtype(i_angle(3,i)))
      factor1=c_a*(r_ab+r_ac-r_ab_cov-r_ac_cov)
      factor2=(r_ab_cov*r_ac_cov)**(d_a)
 1004 hessian(nbonds+i,nbonds+i)=a_a+(b_a/factor2)*exp(-factor1)

      do 1005 i=1,ndihedrals

      l=0

      do 1006 j=1,nbonds
      if(i_bond(1,j).eq.i_dihedral(2,i)) l=l+1
      if(i_bond(1,j).eq.i_dihedral(3,i)) l=l+1
      if(i_bond(2,j).eq.i_dihedral(2,i)) l=l+1
      if(i_bond(2,j).eq.i_dihedral(3,i)) l=l+1
 1006 continue

      r_ab=dist(coord(1,i_dihedral(2,i)),coord(1,i_dihedral(3,i)))
      r_ab_cov=radius(natomtype(i_dihedral(2,i)))
     &        +radius(natomtype(i_dihedral(3,i)))

      factor1=exp(-c_d*(r_ab-r_ab_cov))
      factor2=(b_d*(real64(l)**d_d))/((r_ab*r_ab_cov)**e_d)

 1005 hessian(nbonds+nangles+i,nbonds+nangles+i)=a_d+factor1*factor2

      do 1007 i=1,nimpropers

      theta=dimp(coord(1,i_improper(1,i)),coord(1,i_improper(2,i)),
     &           coord(1,i_improper(3,i)),coord(1,i_improper(4,i)))

      r_ax=dist(coord(1,i_improper(4,i)),coord(1,i_improper(1,i)))

      r_ab_cov=radius(natomtype(i_improper(4,i)))
     &        +radius(natomtype(i_improper(2,i)))
      r_ac_cov=radius(natomtype(i_improper(4,i)))
     &        +radius(natomtype(i_improper(3,i)))
      r_ax_cov=radius(natomtype(i_improper(4,i)))
     &        +radius(natomtype(i_improper(1,i)))

      factor1=exp(-c_i*(r_ax-r_ax_cov))
      factor2=b_i*((r_ab_cov*r_ac_cov)**e_i)*((cos(theta))**d_i)

      ii=nbonds+nangles+ndihedrals+i

 1007 hessian(ii,ii)=a_i+factor1*factor2

      return
      end
