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.
*
* 03/11/20 Adam. Array fromto eliminated, computed on the fly
*     Fixed the problem with vbld indices, which caused errors in
*     derivatives when the backbone virtual bond lengths were not equal.
*********************************************************************** 
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      double precision drt(3,3,maxres),rdt(3,3,maxres),dp(3,3),
     &temp(3,3),prordt(3,3,maxres),prodrt(3,3,maxres)
      double precision xx(3),xx1(3),alphi,omegi,xj,dpjk,yp,xp,xxp,yyp
      double precision cosalphi,sinalphi,cosomegi,sinomegi,theta2,
     & cost2,sint2,rj,dxoiij,tempkl,dxoijk,dsci,zzp,dj,dpkl
      double precision fromto(3,3)
      integer i,ii,j,jjj,k,l,m,indi,ind,ind1
* get the position of the jth ijth fragment of the chain coordinate system      
* in the fromto array.
      integer indmat
      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
      call chainbuild_extconf
      call cartprint
      call intout
*
* 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 
*
* Calculate derivatives.
*
      ind1=0
      do i=1,nres-2
        ind1=ind1+1
*
* Derivatives of DC(i+1) in theta(i+2)
*
c        write (iout,*) "theta i",i
c        write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
c        write(iout,'(7hrdt    9f10.5)')((rdt(k,l,i),l=1,3),k=1,3)
c        write(iout,*) "vbld",vbld(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
c          dcdv(j,ind1)=vbld(i+1)*dp(j,1)       
          dcdv(j,ind1)=vbld(i+2)*dp(j,1)       
        enddo
c        write(iout,'(7hdcdv   3f10.5)')(dcdv(k,ind1),k=1,3)
*
* 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
c        write (iout,*) "dxdv",(dxdv(j,ind1),j=1,3)
*
* 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
c        write(iout,*)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 
c          dcdv(j+3,ind1)=vbld(i+1)*dp(j,1)
          dcdv(j+3,ind1)=vbld(i+2)*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)
c          write(iout,*)'i=',i,' j=',j,' ind=',ind,' ind1=',ind1
          call build_fromto(i+1,j+1,fromto)
c          write(iout,'(7hfromto 9f10.5)')((fromto(k,l),l=1,3),k=1,3)
          do k=1,3
            do l=1,3
              tempkl=0.0D0
              do m=1,2
                tempkl=tempkl+prordt(k,m,i)*fromto(m,l)
              enddo
              temp(k,l)=tempkl
            enddo
          enddo  
c          write(iout,'(7hfromto 9f10.5)')((fromto(k,l,ind),l=1,3),k=1,3)
c          write(iout,'(7hprod   9f10.5)')((prod(k,l,i),l=1,3),k=1,3)
c          write(iout,'(7htemp   9f10.5)')((temp(k,l),l=1,3),k=1,3)
* Derivatives of virtual-bond vectors in theta
          do k=1,3
c            dcdv(k,ind1)=vbld(i+1)*temp(k,1)
            dcdv(k,ind1)=vbld(j+2)*temp(k,1)
          enddo
c          write(iout,'(7hdcdv   3f10.5)')(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
c          write(iout,'(7htheta  3f10.5)')(dxdv(k,ind1),k=1,3)
*
*--- Calculate the derivatives in phi
*
#ifdef FIVEDIAG
          do k=1,3
            do l=1,3
              tempkl=0.0D0
              do m=1,3
                tempkl=tempkl+prodrt(k,m,i)*fromto(m,l)
              enddo
              temp(k,l)=tempkl
            enddo
          enddo
#else
          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
#endif
          do k=1,3
c            dcdv(k+3,ind1)=vbld(i+1)*temp(k,1)
            dcdv(k+3,ind1)=vbld(j+2)*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
#ifdef DEBUG
      write (iout,*)
      write (iout,'(a)') '****************** ddc/dtheta'
      write (iout,*)
      do i=1,nres-2
        do j=i+1,nres-1
	  ii = indmat(i,j)
          write (iout,'(2i4,3e14.6)') i,j,(dcdv(k,ii),k=1,3)
        enddo
      enddo    
      write (iout,*) 
      write (iout,'(a)') '******************* ddc/dphi'
      write (iout,*)
      do i=1,nres-3
        do j=i+2,nres-1
	  ii = indmat(i+1,j)
          write (iout,'(2i4,3e14.6)') i,j,(dcdv(k+3,ii),k=1,3)
	  write (iout,'(a)')
        enddo
      enddo   
      write (iout,'(a)')
      write (iout,'(a)') '**************** dx/dtheta'
      write (iout,'(a)')
      do i=3,nres
        do j=i-1,nres-1
	  ii = indmat(i-2,j)
          write (iout,'(2i4,3e14.6)') i,j,(dxdv(k,ii),k=1,3)
        enddo
      enddo
      write (iout,'(a)')
      write (iout,'(a)') '***************** dx/dphi'
      write (iout,'(a)')
      do i=4,nres
        do j=i-1,nres-1
	  ii = indmat(i-2,j)
          write (iout,'(2i4,3e14.6)') i,j,(dxdv(k+3,ii),k=1,3)
          write(iout,'(a)')
        enddo
      enddo
#endif
*
* 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