subroutine dEconstrQ_num
c Calculating numerical dUconst/ddc and dUconst/ddx      
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.NAMES'
      include 'COMMON.TIME1'
      double precision uzap1,uzap2
      double precision dUcartan(3,0:MAXRES)
     & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
      integer kstart,kend,lstart,lend,idummy
      double precision delta /1.0d-7/
c     For the backbone
      do i=0,nres-1
         do j=1,3
            dUcartan(j,i)=0.0d0
            cdummy(j,i)=dc(j,i)
            dc(j,i)=dc(j,i)+delta
            call chainbuild_cart
	    uzap2=0.0d0
            do ii=1,nfrag
             qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
     &         ,idummy,idummy)
               uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
     &          qinfrag(ii,iset))
            enddo
            do ii=1,npair
               kstart=ifrag(1,ipair(1,ii,iset),iset)
               kend=ifrag(2,ipair(1,ii,iset),iset)
               lstart=ifrag(1,ipair(2,ii,iset),iset)
               lend=ifrag(2,ipair(2,ii,iset),iset)
               qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
               uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
     &           qinpair(ii,iset))
            enddo
            dc(j,i)=cdummy(j,i)
            call chainbuild_cart
            uzap1=0.0d0
             do ii=1,nfrag
             qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
     &         ,idummy,idummy)
               uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
     &          qinfrag(ii,iset))
            enddo
            do ii=1,npair
               kstart=ifrag(1,ipair(1,ii,iset),iset)
               kend=ifrag(2,ipair(1,ii,iset),iset)
               lstart=ifrag(1,ipair(2,ii,iset),iset)
               lend=ifrag(2,ipair(2,ii,iset),iset)
               qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
               uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
     &          qinpair(ii,iset))
            enddo
            ducartan(j,i)=(uzap2-uzap1)/(delta)	    
         enddo
      enddo
c Calculating numerical gradients for dU/ddx
      do i=0,nres-1
         duxcartan(j,i)=0.0d0
         do j=1,3
            cdummy(j,i)=dc(j,i+nres)
            dc(j,i+nres)=dc(j,i+nres)+delta
            call chainbuild_cart
	    uzap2=0.0d0
            do ii=1,nfrag
             qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
     &         ,idummy,idummy)
               uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
     &          qinfrag(ii,iset))
            enddo
            do ii=1,npair
               kstart=ifrag(1,ipair(1,ii,iset),iset)
               kend=ifrag(2,ipair(1,ii,iset),iset)
               lstart=ifrag(1,ipair(2,ii,iset),iset)
               lend=ifrag(2,ipair(2,ii,iset),iset)
               qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
               uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
     &          qinpair(ii,iset))
            enddo
            dc(j,i+nres)=cdummy(j,i)
            call chainbuild_cart
            uzap1=0.0d0
             do ii=1,nfrag
               qfrag(ii)=qwolynes(ifrag(1,ii,iset),
     &          ifrag(2,ii,iset),.true.,idummy,idummy)
               uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
     &          qinfrag(ii,iset))
            enddo
            do ii=1,npair
               kstart=ifrag(1,ipair(1,ii,iset),iset)
               kend=ifrag(2,ipair(1,ii,iset),iset)
               lstart=ifrag(1,ipair(2,ii,iset),iset)
               lend=ifrag(2,ipair(2,ii,iset),iset)
               qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
               uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
     &          qinpair(ii,iset))
            enddo
            duxcartan(j,i)=(uzap2-uzap1)/(delta)	    
         enddo
      enddo    
      write(iout,*) "Numerical dUconst/ddc backbone "
      do ii=0,nres
        write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
      enddo
c      write(iout,*) "Numerical dUconst/ddx side-chain "
c      do ii=1,nres
c         write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
c      enddo 
      return
      end
c---------------------------------------------------------------------------