subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.MD'
      include 'COMMON.IOUNITS'
      external ufparm
      integer uiparm(1)
      double precision urparm(1)
      dimension x(maxvar),g(maxvar)
c
c This subroutine calculates total internal coordinate gradient.
c Depending on the number of function evaluations, either whole energy 
c is evaluated beforehand, Cartesian coordinates and their derivatives in 
c internal coordinates are reevaluated or only the cartesian-in-internal
c coordinate derivatives are evaluated. The subroutine was designed to work
c with SUMSL.
c 
c
      icg=mod(nf,2)+1

cd      print *,'grad',nf,icg
      if (nf-nfl+1) 20,30,40
   20 call func(n,x,nf,f,uiparm,urparm,ufparm)
c     write (iout,*) 'grad 20'
      if (nf.eq.0) return
      goto 40
   30 call var_to_geom(n,x)
      call chainbuild 
c     write (iout,*) 'grad 30'
C
C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
C
   40 call cartder
c     write (iout,*) 'grad 40'
c     print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
C
C Convert the Cartesian gradient into internal-coordinate gradient.
C
      ind=0
      ind1=0
      do i=1,nres-2
	gthetai=0.0D0
	gphii=0.0D0
	do j=i+1,nres-1
          ind=ind+1
c         ind=indmat(i,j)
c         print *,'GRAD: i=',i,' jc=',j,' ind=',ind
	  do k=1,3
            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
          enddo
	  do k=1,3
	    gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
          enddo
        enddo
	do j=i+1,nres-1
          ind1=ind1+1
c         ind1=indmat(i,j)
c         print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
	  do k=1,3
	    gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
	    gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
          enddo
        enddo
	if (i.gt.1) g(i-1)=gphii
	if (n.gt.nphi) g(nphi+i)=gthetai
      enddo
      if (n.le.nphi+ntheta) goto 10
      do i=2,nres-1
	if (itype(i).ne.10) then
          galphai=0.0D0
	  gomegai=0.0D0
	  do k=1,3
	    galphai=galphai+dxds(k,i)*gradx(k,i,icg)
          enddo
	  do k=1,3
	    gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
          enddo
          g(ialph(i,1))=galphai
	  g(ialph(i,1)+nside)=gomegai
        endif
      enddo
C
C Add the components corresponding to local energy terms.
C
   10 continue
      do i=1,nvar
cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
        g(i)=g(i)+gloc(i,icg)
      enddo
C Uncomment following three lines for diagnostics.
cd    call intout
cd    call briefout(0,0.0d0)
cd    write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
      return
      end
C-------------------------------------------------------------------------
      subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.IOUNITS'
      external ufparm
      integer uiparm(1)
      double precision urparm(1)
      dimension x(maxvar),g(maxvar)

      icg=mod(nf,2)+1
      if (nf-nfl+1) 20,30,40
   20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
c     write (iout,*) 'grad 20'
      if (nf.eq.0) return
      goto 40
   30 continue
#ifdef OSF
c     Intercept NaNs in the coordinates
c      write(iout,*) (var(i),i=1,nvar)
      x_sum=0.D0
      do i=1,n
        x_sum=x_sum+x(i)
      enddo
      if (x_sum.ne.x_sum) then
        write(iout,*)" *** grad_restr : Found NaN in coordinates"
        call flush(iout)
        print *," *** grad_restr : Found NaN in coordinates"
        return
      endif
#endif
      call var_to_geom_restr(n,x)
      call chainbuild 
C
C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
C
   40 call cartder
C
C Convert the Cartesian gradient into internal-coordinate gradient.
C

      ig=0
      ind=nres-2                                                                    
      do i=2,nres-2                
       IF (mask_phi(i+2).eq.1) THEN                                             
        gphii=0.0D0                                                             
        do j=i+1,nres-1                                                         
          ind=ind+1                                 
          do k=1,3                                                              
            gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
            gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
          enddo                                                                 
        enddo                                                                   
        ig=ig+1
        g(ig)=gphii
       ELSE
        ind=ind+nres-1-i
       ENDIF
      enddo                                        


      ind=0
      do i=1,nres-2
       IF (mask_theta(i+2).eq.1) THEN
        ig=ig+1
	gthetai=0.0D0
	do j=i+1,nres-1
          ind=ind+1
	  do k=1,3
            gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
            gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
          enddo
        enddo
        g(ig)=gthetai
       ELSE
        ind=ind+nres-1-i
       ENDIF
      enddo

      do i=2,nres-1
	if (itype(i).ne.10) then
         IF (mask_side(i).eq.1) THEN
          ig=ig+1
          galphai=0.0D0
	  do k=1,3
	    galphai=galphai+dxds(k,i)*gradx(k,i,icg)
          enddo
          g(ig)=galphai
         ENDIF
        endif
      enddo

      
      do i=2,nres-1
        if (itype(i).ne.10) then
         IF (mask_side(i).eq.1) THEN
          ig=ig+1
	  gomegai=0.0D0
	  do k=1,3
	    gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
          enddo
	  g(ig)=gomegai
         ENDIF
        endif
      enddo

C
C Add the components corresponding to local energy terms.
C

      ig=0
      igall=0
      do i=4,nres
        igall=igall+1
        if (mask_phi(i).eq.1) then
          ig=ig+1
          g(ig)=g(ig)+gloc(igall,icg)
        endif
      enddo

      do i=3,nres
        igall=igall+1
        if (mask_theta(i).eq.1) then
          ig=ig+1
          g(ig)=g(ig)+gloc(igall,icg)
        endif
      enddo
     
      do ij=1,2
      do i=2,nres-1
        if (itype(i).ne.10) then
          igall=igall+1
          if (mask_side(i).eq.1) then
            ig=ig+1
            g(ig)=g(ig)+gloc(igall,icg)
          endif
        endif
      enddo
      enddo

cd      do i=1,ig
cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
cd      enddo
      return
      end
C-------------------------------------------------------------------------
      subroutine cartgrad
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.VAR'
      include 'COMMON.INTERACT'
      include 'COMMON.FFIELD'
      include 'COMMON.MD'
      include 'COMMON.IOUNITS'
      include 'COMMON.TIME1'
c
c This subrouting calculates total Cartesian coordinate gradient. 
c The subroutine chainbuild_cart and energy MUST be called beforehand.
c
#ifdef TIMING
      time00=MPI_Wtime()
#endif
      icg=1
      call sum_gradient
#ifdef TIMING
#endif
#ifdef DEBUG
      write (iout,*) "After sum_gradient"
      do i=1,nres-1
        write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
        write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
      enddo
#endif
c If performing constraint dynamics, add the gradients of the constraint energy
      if(usampl.and.totT.gt.eq_time) then
         do i=1,nct
           do j=1,3
             gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i)
             gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i)
           enddo
         enddo
         do i=1,nres-3
           gloc(i,icg)=gloc(i,icg)+dugamma(i)
         enddo
         do i=1,nres-2
           gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
         enddo
      endif 
#ifdef TIMING
      time01=MPI_Wtime()
#endif
      call intcartderiv
#ifdef TIMING
      time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
#endif
cd      call checkintcartgrad
cd      write(iout,*) 'calling int_to_cart'
#ifdef DEBUG
      write (iout,*) "gcart, gxcart, gloc before int_to_cart"
#endif
      do i=1,nct
        do j=1,3
          gcart(j,i)=gradc(j,i,icg)
          gxcart(j,i)=gradx(j,i,icg)
        enddo
#ifdef DEBUG
        write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),
     &    (gxcart(j,i),j=1,3),gloc(i,icg)
#endif
      enddo
#ifdef TIMING
      time01=MPI_Wtime()
#endif
      call int_to_cart
#ifdef TIMING
      time_inttocart=time_inttocart+MPI_Wtime()-time01
#endif
#ifdef DEBUG
      write (iout,*) "gcart and gxcart after int_to_cart"
      do i=0,nres-1
        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
     &      (gxcart(j,i),j=1,3)
      enddo
#endif
#ifdef TIMING
      time_cartgrad=time_cartgrad+MPI_Wtime()-time00
#endif
      return
      end
C-------------------------------------------------------------------------
      subroutine zerograd
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.DERIV'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.MD'
C
C Initialize Cartesian-coordinate gradient
C
      do i=1,nres
	do j=1,3
	  gvdwx(j,i)=0.0D0
          gradx_scp(j,i)=0.0D0
	  gvdwc(j,i)=0.0D0
          gvdwc_scp(j,i)=0.0D0
          gvdwc_scpp(j,i)=0.0d0
	  gelc (j,i)=0.0D0
	  gelc_long(j,i)=0.0D0
          gradb(j,i)=0.0d0
          gradbx(j,i)=0.0d0
          gvdwpp(j,i)=0.0d0
          gel_loc(j,i)=0.0d0
          gel_loc_long(j,i)=0.0d0
	  ghpbc(j,i)=0.0D0
	  ghpbx(j,i)=0.0D0
          gcorr3_turn(j,i)=0.0d0
          gcorr4_turn(j,i)=0.0d0
          gradcorr(j,i)=0.0d0
          gradcorr_long(j,i)=0.0d0
          gradcorr5_long(j,i)=0.0d0
          gradcorr6_long(j,i)=0.0d0
          gcorr6_turn_long(j,i)=0.0d0
          gradcorr5(j,i)=0.0d0
          gradcorr6(j,i)=0.0d0
          gcorr6_turn(j,i)=0.0d0
          gsccorc(j,i)=0.0d0
          gsccorx(j,i)=0.0d0
          gradc(j,i,icg)=0.0d0
          gradx(j,i,icg)=0.0d0
          gscloc(j,i)=0.0d0
          gsclocx(j,i)=0.0d0
        enddo
      enddo
C
C Initialize the gradient of local energy terms.
C
      do i=1,4*nres
        gloc(i,icg)=0.0D0
      enddo
      do i=1,nres
        gel_loc_loc(i)=0.0d0
        gcorr_loc(i)=0.0d0
        g_corr5_loc(i)=0.0d0
        g_corr6_loc(i)=0.0d0
        gel_loc_turn3(i)=0.0d0
        gel_loc_turn4(i)=0.0d0
        gel_loc_turn6(i)=0.0d0
        gsccor_loc(i)=0.0d0
      enddo
c initialize gcart and gxcart
      do i=0,nres
        do j=1,3
          gcart(j,i)=0.0d0
          gxcart(j,i)=0.0d0
        enddo
      enddo
      return
      end
c-------------------------------------------------------------------------
      double precision function fdum()
      fdum=0.0D0
      return
      end