subroutine EconstrQ
c     MD with umbrella_sampling using Wolyne's distance measure as a constraint
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
      include 'COMMON.QRESTR'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
#ifdef FIVEDIAG
      include 'COMMON.LANGEVIN.lang0.5diag'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
#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,hm1,hm2,hmnum
      double precision ucdelan,dUcartan(3,0:MAXRES)
     & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
     &  duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
      integer kstart,kend,lstart,lend,idummy
      double precision delta /1.0d-7/
      do i=0,nres
         do j=1,3
            duconst(j,i)=0.0d0
            dudconst(j,i)=0.0d0	    
            duxconst(j,i)=0.0d0
            dudxconst(j,i)=0.0d0	     
         enddo
      enddo
      Uconst=0.0d0
#ifdef DEBUG
      write (iout,*) "econstrq: nfrag",nfrag," iset",iset,
     &   " ifrag",(ifrag(1,i,iset),ifrag(2,i,iset),i=1,nfrag)
      do i=1,nfrag
         qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
     &    ,idummy,idummy)
         write (iout,*) "fragment",i," qfrag",qfrag(i)
#endif
         Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
c Calculating the derivatives of Constraint energy with respect to Q
         Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
     &     qinfrag(i,iset))
c         hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
c      	 hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
c         hmnum=(hm2-hm1)/delta	 	 
c         write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
c     &   qinfrag(i,iset))
c         write(iout,*) "harmonicnum frag", hmnum	  	 
c Calculating the derivatives of Q with respect to cartesian coordinates
         call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
     &   ,idummy,idummy)
c         write(iout,*) "dqwol "
c         do ii=1,nres
c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
c         enddo
c         write(iout,*) "dxqwol "
c         do ii=1,nres
c           write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
c         enddo
c  The gradients of Uconst in Cs
         do ii=0,nres
            do j=1,3
               duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
               dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
            enddo
         enddo
      enddo	
      do i=1,npair
         kstart=ifrag(1,ipair(1,i,iset),iset)
         kend=ifrag(2,ipair(1,i,iset),iset)
         lstart=ifrag(1,ipair(2,i,iset),iset)
         lend=ifrag(2,ipair(2,i,iset),iset)
         qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
         Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
c  Calculating dU/dQ
         Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
c         hm1=harmonic(qpair(i),qinpair(i,iset))
c      	 hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
c         hmnum=(hm2-hm1)/delta	 	 
c         write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
c     &   qinpair(i,iset))
c         write(iout,*) "harmonicnum pair ", hmnum	 
c Calculating dQ/dXi
         call qwolynes_prim(kstart,kend,.false.
     &   ,lstart,lend)
c         write(iout,*) "dqwol "
c         do ii=1,nres
c          write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
c         enddo
c         write(iout,*) "dxqwol "
c         do ii=1,nres
c          write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
c        enddo
c The gradients of Uconst in Cs
         do ii=0,nres
            do j=1,3
               duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
               dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
            enddo
         enddo
      enddo
c      write(iout,*) "Uconst inside subroutine ", Uconst
c Transforming the gradients from Cs to dCs for the backbone
      do i=0,nres
         do j=i+1,nres
           do k=1,3
             dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
           enddo
         enddo
      enddo
c  Transforming the gradients from Cs to dCs for the side chains      
      do i=1,nres
         do j=1,3
           dudxconst(j,i)=duxconst(j,i)
         enddo
      enddo	 	   	 
c      write(iout,*) "dU/ddc backbone "
c       do ii=0,nres
c        write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
c      enddo      
c      write(iout,*) "dU/ddX side chain "
c      do ii=1,nres
c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
c      enddo
      return
      end