      subroutine getdibydc(nint,ncentr,izmat1,izmat2,izmat3,dibydc,
     &                     coord)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

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

      dimension izmat1(*),izmat2(*),izmat3(*)

      dimension dibydc(nint,*),coord(3,*)

      dimension vec1(3),vec2(3),vec3(3),vec4(3),vec5(3)

      data zero,one/0.0d0,1.0d0/

      do 1001 i=1,3*ncentr
      do 1001 j=1,nint
 1001 dibydc(j,i)=zero

      iint=0

      do 1002 natom1=1,ncentr

      if(natom1.eq.1) ndegfreedom=0
      if(natom1.eq.2) ndegfreedom=1
      if(natom1.eq.3) ndegfreedom=2
      if(natom1.ge.4) ndegfreedom=3

      if(ndegfreedom.eq.0) goto 1002

      do 1003 n=1,ndegfreedom

      iint=iint+1

      natom2=izmat1(natom1)
      natom3=izmat2(natom1)
      natom4=izmat3(natom1)

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

      goto (2001,3001,4001) n

 2001 continue

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

      do 2002 k=1,3
      dibydc(iint,iloc1+k)=-vec1(k)
 2002 dibydc(iint,iloc2+k)=+vec1(k)

      goto 1003

 3001 continue

      call unitvec(vec1,coord(1,natom2),coord(1,natom1))
      call unitvec(vec2,coord(1,natom2),coord(1,natom3))

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

      factor1=dot(3,vec1,vec2)
      factor2=sqrt(one-factor1*factor1)

      do 3002 k=1,3
      dibydc(iint,iloc1+k)= (factor1*vec1(k)-vec2(k))/(dist12*factor2)
      dibydc(iint,iloc2+k)=-(factor1*vec1(k)-vec2(k))/(dist12*factor2)
     &                     -(factor1*vec2(k)-vec1(k))/(dist23*factor2)
 3002 dibydc(iint,iloc3+k)= (factor1*vec2(k)-vec1(k))/(dist23*factor2)

      goto 1003

 4001 continue

      call unitvec(vec1,coord(1,natom1),coord(1,natom2))
      call unitvec(vec2,coord(1,natom2),coord(1,natom3))
      call unitvec(vec3,coord(1,natom3),coord(1,natom4))

      factor1=-dot(3,vec1,vec2)
      factor2=-dot(3,vec2,vec3)

      factor3=sqrt(one-factor1*factor1)
      factor4=sqrt(one-factor2*factor2)

      call crossprod(vec4,vec1,vec2)
      call crossprod(vec5,vec2,vec3)

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

      factor5=(dist23-dist12*factor1)/(dist23*dist12*factor3*factor3)
      factor6=(dist23-dist34*factor2)/(dist23*dist34*factor4*factor4)

      do 4002 k=1,3
      dibydc(iint,iloc1+k)=-vec4(k)/(dist12*factor3*factor3)
      dibydc(iint,iloc2+k)=+vec4(k)*factor5-vec5(k)*factor2/
     &                                      (dist23*factor4*factor4)
      dibydc(iint,iloc3+k)=-vec5(k)*factor6+vec4(k)*factor1/
     &                                      (dist23*factor3*factor3)
 4002 dibydc(iint,iloc4+k)=+vec5(k)/(dist34*factor4*factor4)

 1003 continue

 1002 continue

      return
      end
