subroutine cartder 
***********************************************************************
* This subroutine calculates the derivatives of the consecutive virtual
* bond vectors and the SC vectors in the virtual-bond angles theta and
* virtual-torsional angles phi, as well as the derivatives of SC vectors
* in the angles alpha and omega, describing the location of a side chain
* in its local coordinate system.
*
* The derivatives are stored in the following arrays:
*
* DDCDV - the derivatives of virtual-bond vectors DC in theta and phi.
* The structure is as follows:
*
* dDC(x,2)/dT(3),...,dDC(z,2)/dT(3),0,             0,             0
* dDC(x,3)/dT(4),...,dDC(z,3)/dT(4),dDC(x,3)/dP(4),dDC(y,4)/dP(4),dDC(z,4)/dP(4)
*         . . . . . . . . . . . .  . . . . . .
* dDC(x,N-1)/dT(4),...,dDC(z,N-1)/dT(4),dDC(x,N-1)/dP(4),dDC(y,N-1)/dP(4),dDC(z,N-1)/dP(4)
*                          .
*                          .
*                          .
* dDC(x,N-1)/dT(N),...,dDC(z,N-1)/dT(N),dDC(x,N-1)/dP(N),dDC(y,N-1)/dP(N),dDC(z,N-1)/dP(N)
*
* DXDV - the derivatives of the side-chain vectors in theta and phi. 
* The structure is same as above.
*
* DCDS - the derivatives of the side chain vectors in the local spherical
* andgles alph and omega:
*
* dX(x,2)/dA(2),dX(y,2)/dA(2),dX(z,2)/dA(2),dX(x,2)/dO(2),dX(y,2)/dO(2),dX(z,2)/dO(2)
* dX(x,3)/dA(3),dX(y,3)/dA(3),dX(z,3)/dA(3),dX(x,3)/dO(3),dX(y,3)/dO(3),dX(z,3)/dO(3)
*                          .
*                          .
*                          .
* dX(x,N-1)/dA(N-1),dX(y,N-1)/dA(N-1),dX(z,N-1)/dA(N-1),dX(x,N-1)/dO(N-1),dX(y,N-1)/dO(N-1),dX(z,N-1)/dO(N-1)
*
* Version of March '95, based on an early version of November '91.
*
*********************************************************************** 
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      dimension drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),temp(3,3),
     &     fromto(3,3,maxdim),prordt(3,3,maxres),prodrt(3,3,maxres)
      dimension xx(3),xx1(3)
c      common /przechowalnia/ fromto
* get the position of the jth ijth fragment of the chain coordinate system      
* in the fromto array.
      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
*
* calculate the derivatives of transformation matrix elements in theta
*
      do i=1,nres-2
        rdt(1,1,i)=-rt(1,2,i)
        rdt(1,2,i)= rt(1,1,i)
        rdt(1,3,i)= 0.0d0
        rdt(2,1,i)=-rt(2,2,i)
        rdt(2,2,i)= rt(2,1,i)
        rdt(2,3,i)= 0.0d0
        rdt(3,1,i)=-rt(3,2,i)
        rdt(3,2,i)= rt(3,1,i)
        rdt(3,3,i)= 0.0d0
      enddo
*
* derivatives in phi
*
      do i=2,nres-2
        drt(1,1,i)= 0.0d0
        drt(1,2,i)= 0.0d0
        drt(1,3,i)= 0.0d0
        drt(2,1,i)= rt(3,1,i)
        drt(2,2,i)= rt(3,2,i)
        drt(2,3,i)= rt(3,3,i)
        drt(3,1,i)=-rt(2,1,i)
        drt(3,2,i)=-rt(2,2,i)
        drt(3,3,i)=-rt(2,3,i)
      enddo 
*
* generate the matrix products of type r(i)t(i)...r(j)t(j)
*
      do i=2,nres-2
        ind=indmat(i,i+1)
        do k=1,3
          do l=1,3
            temp(k,l)=rt(k,l,i)
          enddo
        enddo
        do k=1,3
          do l=1,3
            fromto(k,l,ind)=temp(k,l)
          enddo
        enddo  
        do j=i+1,nres-2
          ind=indmat(i,j+1)
          do k=1,3
            do l=1,3
              dpkl=0.0d0
              do m=1,3
                dpkl=dpkl+temp(k,m)*rt(m,l,j)
              enddo
              dp(k,l)=dpkl
              fromto(k,l,ind)=dpkl
            enddo
          enddo
          do k=1,3
            do l=1,3
              temp(k,l)=dp(k,l)
            enddo
          enddo
        enddo
      enddo
*
* Calculate derivatives.
*
      ind1=0
      do i=1,nres-2
	ind1=ind1+1
*
* Derivatives of DC(i+1) in theta(i+2)
*
        do j=1,3
          do k=1,2
            dpjk=0.0D0
            do l=1,3
              dpjk=dpjk+prod(j,l,i)*rdt(l,k,i)
            enddo
            dp(j,k)=dpjk
            prordt(j,k,i)=dp(j,k)
          enddo
          dp(j,3)=0.0D0
          dcdv(j,ind1)=vbld(i+1)*dp(j,1)       
        enddo
*
* Derivatives of SC(i+1) in theta(i+2)
* 
        xx1(1)=-0.5D0*xloc(2,i+1)
        xx1(2)= 0.5D0*xloc(1,i+1)
        do j=1,3
          xj=0.0D0
          do k=1,2
            xj=xj+r(j,k,i)*xx1(k)
          enddo
          xx(j)=xj
        enddo
        do j=1,3
          rj=0.0D0
          do k=1,3
            rj=rj+prod(j,k,i)*xx(k)
          enddo
          dxdv(j,ind1)=rj
        enddo
*
* Derivatives of SC(i+1) in theta(i+3). The have to be handled differently
* than the other off-diagonal derivatives.
*
        do j=1,3
          dxoiij=0.0D0
          do k=1,3
            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
          enddo
          dxdv(j,ind1+1)=dxoiij
        enddo
cd      print *,ind1+1,(dxdv(j,ind1+1),j=1,3)
*
* Derivatives of DC(i+1) in phi(i+2)
*
        do j=1,3
          do k=1,3
            dpjk=0.0
            do l=2,3
              dpjk=dpjk+prod(j,l,i)*drt(l,k,i)
            enddo
            dp(j,k)=dpjk
            prodrt(j,k,i)=dp(j,k)
          enddo 
          dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
        enddo
*
* Derivatives of SC(i+1) in phi(i+2)
*
        xx(1)= 0.0D0 
        xx(3)= xloc(2,i+1)*r(2,2,i)+xloc(3,i+1)*r(2,3,i)
        xx(2)=-xloc(2,i+1)*r(3,2,i)-xloc(3,i+1)*r(3,3,i)
        do j=1,3
          rj=0.0D0
          do k=2,3
            rj=rj+prod(j,k,i)*xx(k)
          enddo
          dxdv(j+3,ind1)=-rj
        enddo
*
* Derivatives of SC(i+1) in phi(i+3).
*
        do j=1,3
          dxoiij=0.0D0
          do k=1,3
            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
          enddo
          dxdv(j+3,ind1+1)=dxoiij
        enddo
*
* Calculate the derivatives of DC(i+1) and SC(i+1) in theta(i+3) thru 
* theta(nres) and phi(i+3) thru phi(nres).
*
        do j=i+1,nres-2
	  ind1=ind1+1
	  ind=indmat(i+1,j+1)
cd        print *,'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
          do k=1,3
            do l=1,3
              tempkl=0.0D0
              do m=1,2
                tempkl=tempkl+prordt(k,m,i)*fromto(m,l,ind)
              enddo
              temp(k,l)=tempkl
            enddo
          enddo  
cd        print '(9f8.3)',((fromto(k,l,ind),l=1,3),k=1,3)
cd        print '(9f8.3)',((prod(k,l,i),l=1,3),k=1,3)
cd        print '(9f8.3)',((temp(k,l),l=1,3),k=1,3)
* Derivatives of virtual-bond vectors in theta
          do k=1,3
            dcdv(k,ind1)=vbld(i+1)*temp(k,1)
          enddo
cd        print '(3f8.3)',(dcdv(k,ind1),k=1,3)
* Derivatives of SC vectors in theta
          do k=1,3
            dxoijk=0.0D0
            do l=1,3
              dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
            enddo
            dxdv(k,ind1+1)=dxoijk
          enddo
*
*--- Calculate the derivatives in phi
*
          do k=1,3
            do l=1,3
              tempkl=0.0D0
              do m=1,3
                tempkl=tempkl+prodrt(k,m,i)*fromto(m,l,ind)
              enddo
              temp(k,l)=tempkl
            enddo
          enddo
          do k=1,3
            dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
	  enddo
          do k=1,3
            dxoijk=0.0D0
            do l=1,3
              dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
            enddo
            dxdv(k+3,ind1+1)=dxoijk
          enddo
        enddo
      enddo
*
* Derivatives in alpha and omega:
*
      do i=2,nres-1
c       dsci=dsc(itype(i))
        dsci=vbld(i+nres)
#ifdef OSF
        alphi=alph(i)
        omegi=omeg(i)
        if(alphi.ne.alphi) alphi=100.0 
        if(omegi.ne.omegi) omegi=-100.0
#else
	alphi=alph(i)
	omegi=omeg(i)
#endif
cd      print *,'i=',i,' dsci=',dsci,' alphi=',alphi,' omegi=',omegi
	cosalphi=dcos(alphi)
	sinalphi=dsin(alphi)
	cosomegi=dcos(omegi)
	sinomegi=dsin(omegi)
	temp(1,1)=-dsci*sinalphi
	temp(2,1)= dsci*cosalphi*cosomegi
	temp(3,1)=-dsci*cosalphi*sinomegi
	temp(1,2)=0.0D0
	temp(2,2)=-dsci*sinalphi*sinomegi
	temp(3,2)=-dsci*sinalphi*cosomegi
	theta2=pi-0.5D0*theta(i+1)
	cost2=dcos(theta2)
	sint2=dsin(theta2)
	jjj=0
cd      print *,((temp(l,k),l=1,3),k=1,2)
        do j=1,2
	  xp=temp(1,j)
	  yp=temp(2,j)
	  xxp= xp*cost2+yp*sint2
	  yyp=-xp*sint2+yp*cost2
	  zzp=temp(3,j)
	  xx(1)=xxp
	  xx(2)=yyp*r(2,2,i-1)+zzp*r(2,3,i-1)
	  xx(3)=yyp*r(3,2,i-1)+zzp*r(3,3,i-1)
	  do k=1,3
	    dj=0.0D0
	    do l=1,3
	      dj=dj+prod(k,l,i-1)*xx(l)
            enddo
	    dxds(jjj+k,i)=dj
          enddo
	  jjj=jjj+3
	enddo
      enddo
      return
      end