subroutine cart2intgrad(n,g) 
***********************************************************************
* This subroutine thransforms the gradient in virtual-bond vectors to
* that in the backbone and side-chain angular variables.
* Adapted from the cartder subroutine.
*
* 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'
      integer n
      double precision g(n)
      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),aux(6)
      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.
c      integer indmat
c      indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
c      call chainbuild_extconf
c      call cartprint
c      call intout
      g=0.0d0
* 3/13/20 Adam: Skip calculating backbone derivatives if SC only
* requested.
      if (sideonly) goto 10
*
* 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 backbone derivatives.
* This code invlves N^2 effort and should be parallelized, to be done
* later.
      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)
        if (n.gt.nphi) then

        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+2)*dp(j,1)       
          g(nphi+i)=g(nphi+i)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
        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
c          dxdv(j,ind1)=rj
c          write (iout,*) "1:i",i," j",i+1,"ind1",ind1," dxdthet",rj,
c     &      " gradx",gradx(j,i+1,icg)
          g(nphi+i)=g(nphi+i)+rj*gradx(j,i+1,icg)
        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.
*
        if (i.lt.nres-2) then
        do j=1,3
          dxoiij=0.0D0
          do k=1,3
            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
          enddo
c          dxdv(j,ind1+1)=dxoiij
c          write (iout,*) "2:i",i," j",i+1,"ind1",ind1+1,
c     &      " dxdthet",dxoiij," gradx",gradx(j,i+2,icg)
          g(nphi+i)=g(nphi+i)+dxoiij*gradx(j,i+2,icg)
        enddo
        endif
c        write(iout,*)ind1+1,(dxdv(j,ind1+1),j=1,3)

        endif
*
* Derivatives of DC(i+1) in phi(i+2)
*
        if (i.gt.1) then
        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+2)*dp(j,1)
          g(i-1)=g(i-1)+vbld(i+2)*dp(j,1)*gradc(j,i+1,icg)
        enddo
        endif
*
* 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)
        if (i.gt.1) then
        do j=1,3
          rj=0.0D0
          do k=2,3
            rj=rj+prod(j,k,i)*xx(k)
          enddo
c          dxdv(j+3,ind1)=-rj
c          write (iout,*) "1:i",i," j",i+1,"ind1",ind1," dxdphi",-rj,
c     &      " gradx",gradx(j,i+1,icg)
          g(i-1)=g(i-1)-rj*gradx(j,i+1,icg)
        enddo
        endif
*
* Derivatives of SC(i+1) in phi(i+3).
*
        if (i.gt.1) then
        do j=1,3
          dxoiij=0.0D0
          do k=1,3
            dxoiij=dxoiij+dp(j,k)*xrot(k,i+2)
          enddo
c          dxdv(j+3,ind1+1)=dxoiij
          g(i-1)=g(i-1)+dxoiij*gradx(j,i+2,icg)
c          write (iout,*) "2:i",i," j",i+2," ind1",ind1+1,
c     &      " dxdphi",dxoiij," gradx",gradx(j,i+2,icg)
        enddo
        endif
*
* 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
c          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)
          if (n.gt.nphi) then
* Derivatives of virtual-bond vectors in theta
          do k=1,3
c            dcdv(k,ind1)=vbld(j+2)*temp(k,1)
            g(nphi+i)=g(nphi+i)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
          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
c            dxdv(k,ind1+1)=dxoijk
c            write (iout,*) "3:i",i+1," j",j+2,"ind1",ind1+1,
c     &      " dxdthet",dxoijk," gradx",gradx(k,j+2,icg)
            g(nphi+i)=g(nphi+i)+dxoijk*gradx(k,j+2,icg)
          enddo
c          write(iout,'(7htheta  3f10.5)')(dxdv(k,ind1),k=1,3)
          endif
*
*--- 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)
              enddo
              temp(k,l)=tempkl
            enddo
          enddo
          if (i.gt.1) then
          do k=1,3
c            dcdv(k+3,ind1)=vbld(j+2)*temp(k,1)
            g(i-1)=g(i-1)+vbld(j+2)*temp(k,1)*gradc(k,j+1,icg)
          enddo
          do k=1,3
            dxoijk=0.0D0
            do l=1,3
              dxoijk=dxoijk+temp(k,l)*xrot(l,j+2)
            enddo
c            dxdv(k+3,ind1+1)=dxoijk
            g(i-1)=g(i-1)+dxoijk*gradx(k,j+2,icg)
c            write (iout,*) "3:i",i," j",j+2," ind1",ind1+1,
c     &      " dxdphi",dxoijk," gradx",gradx(k,j+2,icg)
          enddo
          endif
        enddo
      enddo

      if (nvar.le.nphi+ntheta) return

   10 continue
*
* Derivatives in alpha and omega:
*
      do i=2,nres-1
        if (iabs(itype(i)).eq.10 .or. itype(i).eq.ntyp1!) cycle
     &    .or. mask_side(i).eq.0 ) cycle
        ii=ialph(i,1)
        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
c            dxds(jjj+k,i)=dj
            aux(jjj+k)=dj
          enddo
          jjj=jjj+3
        enddo
        do k=1,3
          g(ii)=g(ii)+aux(k)*gradx(k,i,icg)
          g(ii+nside)=g(ii+nside)+aux(k+3)*gradx(k,i,icg)
        enddo
      enddo
      return
      end
c-----------------------------------------------------------------------------
      subroutine build_fromto(i,j,fromto)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      integer i,j,jj,k,l,m
      double precision fromto(3,3),temp(3,3),dp(3,3)
      double precision dpkl
      save temp
*
* generate the matrix products of type r(i)t(i)...r(j)t(j) on the fly
*
c      write (iout,*) "temp on entry"
c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
c      do i=2,nres-2
c        ind=indmat(i,i+1)
      if (j.eq.i+1) then
        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)=temp(k,l)
          enddo
        enddo  
      else
c        do j=i+1,nres-2
c          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-1)
              enddo
              dp(k,l)=dpkl
              fromto(k,l)=dpkl
            enddo
          enddo
          do k=1,3
            do l=1,3
              temp(k,l)=dp(k,l)
            enddo
          enddo
      endif
c      write (iout,*) "temp upon exit"
c      write (iout,'(3f10.5)') ((temp(k,l),l=1,3),k=1,3)
c        enddo
c      enddo
      return
      end