      subroutine get_bmatrix(natoms,nmodes,nbonds,nangles,ndihedrals,
     &                       nimpropers,i_bond,i_angle,i_dihedral,
     &                       i_improper,coord,bmatrix)

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 i_bond(2,*),i_angle(3,*),i_dihedral(4,*),i_improper(4,*)

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

      dimension vec12(3),vec31(3),vec32(3),vec23(3),vec34(3),vec43(3),
     &          vec21(3),vec41(3),vec42(3),veca(3),vecb(3),vecc(3)

      data zero,one/0.0d0,1.0d0/

      do 1001 i=1,3*natoms
      do 1001 j=1,nmodes
 1001 bmatrix(j,i)=zero

      do 1002 i=1,nbonds

      natom1=i_bond(1,i)
      natom2=i_bond(2,i)

      iloc1=3*(natom1-1)
      iloc2=3*(natom2-1)

      call unitvec(vec12,coord(1,natom1),coord(1,natom2))

      do 1003 k=1,3
      bmatrix(i,iloc1+k)=-vec12(k)
 1003 bmatrix(i,iloc2+k)=+vec12(k)

 1002 continue

      do 1004 i=1,nangles

      natom1=i_angle(1,i)
      natom2=i_angle(2,i)
      natom3=i_angle(3,i)

      iloc1=3*(natom1-1)
      iloc2=3*(natom2-1)
      iloc3=3*(natom3-1)

      call unitvec(vec31,coord(1,natom2),coord(1,natom1))
      call unitvec(vec32,coord(1,natom2),coord(1,natom3))

      dist31=dist(coord(1,natom2),coord(1,natom1))
      dist32=dist(coord(1,natom2),coord(1,natom3))

      cos_phi=dot(3,vec31,vec32)
      sin_phi=sqrt(one-cos_phi*cos_phi)

      do 1005 k=1,3
      bmatrix(nbonds+i,iloc1+k)=+(cos_phi*vec31(k)-vec32(k))
     &                          /(dist31*sin_phi)
      bmatrix(nbonds+i,iloc2+k)=+((dist31-dist32*cos_phi)*vec31(k)
     &                           +(dist32-dist31*cos_phi)*vec32(k))
     &                          /(dist31*dist32*sin_phi)
 1005 bmatrix(nbonds+i,iloc3+k)=+(cos_phi*vec32(k)-vec31(k))
     &                          /(dist32*sin_phi)

 1004 continue

      do 1006 i=1,ndihedrals

      natom1=i_dihedral(1,i)
      natom2=i_dihedral(2,i)
      natom3=i_dihedral(3,i)
      natom4=i_dihedral(4,i)

      iloc1=3*(natom1-1)
      iloc2=3*(natom2-1)
      iloc3=3*(natom3-1)
      iloc4=3*(natom4-1)

      call unitvec(vec12,coord(1,natom1),coord(1,natom2))
      call unitvec(vec23,coord(1,natom2),coord(1,natom3))
      call unitvec(vec34,coord(1,natom3),coord(1,natom4))

      call unitvec(vec43,coord(1,natom4),coord(1,natom3))
      call unitvec(vec32,coord(1,natom3),coord(1,natom2))
      call unitvec(vec21,coord(1,natom2),coord(1,natom1))

      cos_phi2=dot(3,vec21,vec23)
      cos_phi3=dot(3,vec32,vec34)

      sin_phi2=sqrt(one-cos_phi2*cos_phi2)
      sin_phi3=sqrt(one-cos_phi3*cos_phi3)

      call crossprod(veca,vec12,vec23)
      call crossprod(vecb,vec43,vec32)

      dist12=dist(coord(1,natom1),coord(1,natom2))
      dist23=dist(coord(1,natom2),coord(1,natom3))
      dist34=dist(coord(1,natom3),coord(1,natom4))

      factor1=(dist23-dist12*cos_phi2)/(dist23*dist12*sin_phi2*sin_phi2)
      factor2=cos_phi3/(dist23*sin_phi3*sin_phi3)
      factor3=(dist23-dist34*cos_phi3)/(dist23*dist34*sin_phi3*sin_phi3)
      factor4=cos_phi2/(dist23*sin_phi2*sin_phi2)

      do 1007 k=1,3
      bmatrix(nbonds+nangles+i,iloc1+k)=-veca(k)
     &                                  /(dist12*sin_phi2*sin_phi2)
      bmatrix(nbonds+nangles+i,iloc2+k)=factor1*veca(k)+factor2*vecb(k)
      bmatrix(nbonds+nangles+i,iloc3+k)=factor3*vecb(k)+factor4*veca(k)
 1007 bmatrix(nbonds+nangles+i,iloc4+k)=-vecb(k)
     &                                  /(dist34*sin_phi3*sin_phi3)

 1006 continue

      do 1008 i=1,nimpropers

      natom1=i_improper(1,i)
      natom2=i_improper(2,i)
      natom3=i_improper(3,i)
      natom4=i_improper(4,i)

      iloc1=3*(natom1-1)
      iloc2=3*(natom2-1)
      iloc3=3*(natom3-1)
      iloc4=3*(natom4-1)

      phi_1=dang(coord(1,natom2),coord(1,natom4),coord(1,natom3))

      theta=dimp(coord(1,natom1),coord(1,natom2),
     &           coord(1,natom3),coord(1,natom4))

      den=cos(theta)*sin(phi_1)

      fact1=tan(theta)
      fact2=tan(theta)/(sin(phi_1)**2)
      fact3=cos(phi_1)

      d41=one/dist(coord(1,natom4),coord(1,natom1))
      d42=one/dist(coord(1,natom4),coord(1,natom2))
      d43=one/dist(coord(1,natom4),coord(1,natom3))

      call unitvec(vec41,coord(1,natom4),coord(1,natom1))
      call unitvec(vec42,coord(1,natom4),coord(1,natom2))
      call unitvec(vec43,coord(1,natom4),coord(1,natom3))

      call crossprod(veca,vec42,vec43)
      call crossprod(vecb,vec43,vec41)
      call crossprod(vecc,vec41,vec42)

      ii=nbonds+nangles+ndihedrals+i

      do 1009 k=1,3
      bmatrix(ii,iloc1+k)=+d41*((veca(k)/den)-fact1*vec41(k))
      bmatrix(ii,iloc2+k)=+d42*((vecb(k)/den)
     &                    -fact2*(vec42(k)-fact3*vec43(k)))
      bmatrix(ii,iloc3+k)=+d43*((vecc(k)/den)
     &                    -fact2*(vec43(k)-fact3*vec42(k)))
 1009 bmatrix(ii,iloc4+k)=-d41*((veca(k)/den)-fact1*vec41(k))
     &                    -d42*((vecb(k)/den)
     &                    -fact2*(vec42(k)-fact3*vec43(k)))
     &                    -d43*((vecc(k)/den)
     &                    -fact2*(vec43(k)-fact3*vec42(k)))

 1008 continue

      return
      end
