double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      integer i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,
     & secseg
      integer nsep /3/
      double precision dist,qm
      double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
      logical lprn /.false./
      logical flag
      double precision sigm,x
      sigm(x)=0.25d0*x
      qq = 0.0d0
      nl=0 
       if(flag) then
        do il=seg1+nsep,seg2
          do jl=seg1,il-nsep
            nl=nl+1
            d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
     &                 (cref(2,jl)-cref(2,il))**2+
     &                 (cref(3,jl)-cref(3,il))**2)
            dij=dist(il,jl)
            qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
            if (itype(il).ne.10 .or. itype(jl).ne.10) then
              nl=nl+1
              d0ijCM=dsqrt(
     &               (cref(1,jl+nres)-cref(1,il+nres))**2+
     &               (cref(2,jl+nres)-cref(2,il+nres))**2+
     &               (cref(3,jl+nres)-cref(3,il+nres))**2)
              dijCM=dist(il+nres,jl+nres)
              qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
            endif
            qq = qq+qqij+qqijCM
          enddo
        enddo 	
        qq = qq/nl
      else
      do il=seg1,seg2
        if((seg3-il).lt.3) then
             secseg=il+3
        else
             secseg=seg3
        endif 
          do jl=secseg,seg4
            nl=nl+1
            d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
     &                 (cref(2,jl)-cref(2,il))**2+
     &                 (cref(3,jl)-cref(3,il))**2)
            dij=dist(il,jl)
            qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
            if (itype(il).ne.10 .or. itype(jl).ne.10) then
              nl=nl+1
              d0ijCM=dsqrt(
     &               (cref(1,jl+nres)-cref(1,il+nres))**2+
     &               (cref(2,jl+nres)-cref(2,il+nres))**2+
     &               (cref(3,jl+nres)-cref(3,il+nres))**2)
              dijCM=dist(il+nres,jl+nres)
              qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
            endif
            qq = qq+qqij+qqijCM
          enddo
        enddo
      qq = qq/nl
      endif
      qwolynes=1.0d0-qq
      return 
      end
c-------------------------------------------------------------------
      subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.MD'
      integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
     & secseg
      integer nsep /3/
      double precision dist
      double precision dij,d0ij,dijCM,d0ijCM
      logical lprn /.false./
      logical flag
      double precision sigm,x,sim,dd0,fac,ddqij
      sigm(x)=0.25d0*x
      
      do i=0,nres
        do j=1,3
          dqwol(j,i)=0.0d0
          dxqwol(j,i)=0.0d0	  
        enddo
      enddo
      nl=0 
       if(flag) then
        do il=seg1+nsep,seg2
          do jl=seg1,il-nsep
            nl=nl+1
            d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
     &                 (cref(2,jl)-cref(2,il))**2+
     &                 (cref(3,jl)-cref(3,il))**2)
            dij=dist(il,jl)
            sim = 1.0d0/sigm(d0ij)
            sim = sim*sim
            dd0 = dij-d0ij
            fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
	    do k=1,3
              ddqij = (c(k,il)-c(k,jl))*fac
              dqwol(k,il)=dqwol(k,il)+ddqij
              dqwol(k,jl)=dqwol(k,jl)-ddqij
            enddo
            	     
            if (itype(il).ne.10 .or. itype(jl).ne.10) then
              nl=nl+1
              d0ijCM=dsqrt(
     &               (cref(1,jl+nres)-cref(1,il+nres))**2+
     &               (cref(2,jl+nres)-cref(2,il+nres))**2+
     &               (cref(3,jl+nres)-cref(3,il+nres))**2)
              dijCM=dist(il+nres,jl+nres)
              sim = 1.0d0/sigm(d0ijCM)
              sim = sim*sim
              dd0=dijCM-d0ijCM
              fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
              do k=1,3
                ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
                dxqwol(k,il)=dxqwol(k,il)+ddqij
                dxqwol(k,jl)=dxqwol(k,jl)-ddqij
              enddo
            endif 	    
          enddo
        enddo 	
       else
        do il=seg1,seg2
        if((seg3-il).lt.3) then
             secseg=il+3
        else
             secseg=seg3
        endif 
          do jl=secseg,seg4
            nl=nl+1
            d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
     &                 (cref(2,jl)-cref(2,il))**2+
     &                 (cref(3,jl)-cref(3,il))**2)
            dij=dist(il,jl)
            sim = 1.0d0/sigm(d0ij)
            sim = sim*sim
            dd0 = dij-d0ij
            fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
            do k=1,3
              ddqij = (c(k,il)-c(k,jl))*fac
              dqwol(k,il)=dqwol(k,il)+ddqij
              dqwol(k,jl)=dqwol(k,jl)-ddqij
            enddo
            if (itype(il).ne.10 .or. itype(jl).ne.10) then
              nl=nl+1
              d0ijCM=dsqrt(
     &               (cref(1,jl+nres)-cref(1,il+nres))**2+
     &               (cref(2,jl+nres)-cref(2,il+nres))**2+
     &               (cref(3,jl+nres)-cref(3,il+nres))**2)
              dijCM=dist(il+nres,jl+nres)
              sim = 1.0d0/sigm(d0ijCM)
              sim=sim*sim
              dd0 = dijCM-d0ijCM
              fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
              do k=1,3
               ddqij = (c(k,il+nres)-c(k,jl+nres))*fac             
               dxqwol(k,il)=dxqwol(k,il)+ddqij
               dxqwol(k,jl)=dxqwol(k,jl)-ddqij  
              enddo
            endif 
          enddo
        enddo	  	     
      endif
       do i=0,nres
         do j=1,3
           dqwol(j,i)=dqwol(j,i)/nl
           dxqwol(j,i)=dxqwol(j,i)/nl
         enddo
       enddo   	  	  	  	 	    	    	  	 
      return 
      end
c-------------------------------------------------------------------
      subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      integer seg1,seg2,seg3,seg4
      logical flag
      double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
     & qwolxan(3,0:maxres),q1,q2
      double precision delta /1.0d-10/
      do i=0,nres
        do j=1,3
          q1=qwolynes(seg1,seg2,flag,seg3,seg4)
          cdummy(j,i)=c(j,i)
          c(j,i)=c(j,i)+delta
          q2=qwolynes(seg1,seg2,flag,seg3,seg4)
          qwolan(j,i)=(q2-q1)/delta
          c(j,i)=cdummy(j,i)
        enddo
      enddo
      do i=0,nres
        do j=1,3
          q1=qwolynes(seg1,seg2,flag,seg3,seg4)
          cdummy(j,i+nres)=c(j,i+nres)
          c(j,i+nres)=c(j,i+nres)+delta
          q2=qwolynes(seg1,seg2,flag,seg3,seg4)
          qwolxan(j,i)=(q2-q1)/delta
          c(j,i+nres)=cdummy(j,i+nres)
        enddo
      enddo  
c      write(iout,*) "Numerical Q carteisan gradients backbone: "
c      do i=0,nct
c        write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
c      enddo
c      write(iout,*) "Numerical Q carteisan gradients side-chain: "
c      do i=0,nct
c        write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
c      enddo 
      return
      end
c------------------------------------------------------------------------  
      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'
#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,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
      do i=1,nfrag
         qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
     &    ,idummy,idummy)
         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 Calculating numerical gradients of dU/dQi and dQi/dxi
c        call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
c     &  ,idummy,idummy)
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 Calculating numerical gradients
c        call qwol_num(kstart,kend,.false.
c     &  ,lstart,lend)
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
c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
c      call dEconstrQ_num      
      return
      end
c-----------------------------------------------------------------------
      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---------------------------------------------------------------------------