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'
      include 'COMMON.MD'
      integer i,j,jl,k,l,il,kl,nl,np,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
      qq = 0.0d0
      nl=0 
      do i=0,nres
        do j=1,3
          dqwol(j,i)=0.0d0
          dxqwol(j,i)=0.0d0
        enddo
      enddo 
      if (lprn) then
      write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4,
     & " flag",flag
      call flush(iout)
      endif
      if (flag) then
        do il=seg1+nsep,seg2
          do jl=seg1,il-nsep
            nl=nl+1
            if (itype(il).ne.10) then
              ilnres=il+nres
            else
              ilnres=il
            endif
            if (itype(jl).ne.10) then
              jlnres=jl+nres
            else
              jlnres=jl
            endif
            qqijCM = qcontrib(il,jl,ilnres,jlnres)
            qq = qq+qqijCM
            if (lprn) then
              write (iout,*) "qqijCM",qqijCM
              call flush(iout)
            endif
          enddo
        enddo
        if (lprn) then
          write (iout,*) "nl",nl," qq",qq
          call flush(iout)
        endif 
      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
            if (itype(il).ne.10) then
              ilnres=il+nres
            else
              ilnres=il
            endif
            if (itype(jl).ne.10) then
              jlnres=jl+nres
            else
              jlnres=jl
            endif
            qqijCM = qcontrib(il,jl,ilnres,jlnres)
            qq = qq+qqijCM
            if (lprn) then
              write (iout,*) "qqijCM",qqijCM
              call flush(iout)
            endif
          enddo
        enddo
      endif
      qq = qq/nl
      qwolynes=1.0d0-qq
      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'
      include 'COMMON.MD'
      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-7/
      write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4
      write(iout,*) "dQ/dc backbone "
       do i=0,nres
        write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3)
      enddo      
      write(iout,*) "dQ/dX side chain "
      do i=1,nres
            write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3)
      enddo
      do i=1,nres
        do j=1,3
          cdummy(j,i)=c(j,i)
          c(j,i)=c(j,i)-delta
          q1=qwolynes(seg1,seg2,flag,seg3,seg4)
          c(j,i)=cdummy(j,i)+delta
          q2=qwolynes(seg1,seg2,flag,seg3,seg4)
          qwolan(j,i)=0.5d0*(q2-q1)/delta
          c(j,i)=cdummy(j,i)
c          write (iout,*) "i",i," j",j," q1",q1," a2",q2
        enddo
      enddo
      do i=1,nres
        do j=1,3
          cdummy(j,i+nres)=c(j,i+nres)
          c(j,i+nres)=c(j,i+nres)-delta
          q1=qwolynes(seg1,seg2,flag,seg3,seg4)
          c(j,i+nres)=cdummy(j,i+nres)+delta
          q2=qwolynes(seg1,seg2,flag,seg3,seg4)
          qwolxan(j,i)=0.5d0*(q2-q1)/delta
          c(j,i+nres)=cdummy(j,i+nres)
        enddo
      enddo  
      write(iout,*) "Numerical Q cartesian gradients backbone: "
      do i=0,nres
        write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
      enddo
      write(iout,*) "Numerical Q cartesian gradients side-chain: "
      do i=0,nres
        write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
      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 Calculating the derivatives of Q with respect to cartesian coordinates
         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
c      write (iout,*) "Calling qwol_num"
c      call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.,idummy,idummy)
      enddo
c      stop
      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 Calculating dQ/dXi
         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/dc backbone "
c       do ii=0,nres
c        write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
c      enddo      
c      write(iout,*) "dU/dX side chain "
c      do ii=1,nres
c            write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
c      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,(dudxconst(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
         do j=1,3
           duxcartan(j,i)=0.0d0
         enddo
         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
      write(iout,*) "Numerical dUconst/ddx side-chain "
      do ii=1,nres
         write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
      enddo 
      return
      end
c--------------------------------------------------------------------------- 
      double precision function qcontrib(il,jl,il1,jl1)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.MD'
      include 'COMMON.LOCAL'
      integer i,j,k,il,jl,il1,jl1,nd,itl,jtl
      double precision dist
      external dist
      double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120
     &  ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12
      double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3),
     &  aux1,aux2
      double precision scalar
      external scalar
      logical lprn /.false./
      if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1
      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)
      dij1=dist(il1,jl1)
      do i=1,3
        er(i)=(c(i,jl1)-c(i,il1))/dij1
      enddo
      do i=1,3
        er0(i)=cref(i,jl1)-cref(i,il1)
      enddo
      d0ij1=dsqrt(scalar(er0,er0))
      do i=1,3
        er0(i)=er0(i)/d0ij1
      enddo
      if (il.ne.il1 .or. jl.ne.jl1) then
        ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2)
        nd=2
      else
        ddave=(dij-d0ij)**2
        nd=1
      endif
      if (il.ne.il1) then
        do i=1,3
          u(i)=cref(i,il1)-cref(i,il)
        enddo
        d0ii1=dsqrt(scalar(u,u))
        do i=1,3
          u(i)=u(i)/d0ii1
        enddo
        if (lprn) then
        write (iout,*) "u",(u(i),i=1,3)
        write (iout,*) "er0",(er0(i),i=1,3)
        om10=scalar(er0,u)
        om1=scalar(er,dc_norm(1,il1))
        write (iout,*) "om10",om10," om1",om1
        endif
      else
        om1=0.0d0
        om10=0.0d0
      endif
      if (jl.ne.jl1) then
        do i=1,3
          v(i)=cref(i,jl1)-cref(i,jl)
        enddo
        d0jj1=dsqrt(scalar(v,v))
        do i=1,3
          v(i)=v(i)/d0jj1
        enddo
        if (lprn) then
        write (iout,*) "v",(v(i),i=1,3)
        write (iout,*) "er0",(er0(i),i=1,3)
        om20=scalar(er,v)
        om2=scalar(er,dc_norm(1,jl1))
        write (iout,*) "om20",om20," om2",om2
        endif
      else
        om2=0.0d0
        om20=0.0d0
      endif
      if (il.ne.il1 .and. jl.ne.jl1) then
        om120=scalar(u,v)
        om12=scalar(dc_norm(1,il1),dc_norm(1,jl1))
      else
        om12=0.0d0
        om120=0.0d0
      endif
      if (lprn) then
        write (iout,*) "il",il," jl",jl,itype(il),itype(jl)
        write (iout,*)"d0ij",d0ij," om10",om10," om20",om20,
     &   " om120",om120,
     &  " dij",dij," om1",om1," om2",om2," om12",om12
        call flush(iout)
      endif
      ssij = 16.0d0/(d0ij*d0ij)
      qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2
     &                       +(om2-om20)**2+(om12-om120)**2)))
      if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib
c      qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2)
c      qcontrib = dexp(-0.5d0*(ddave*ssij))
c Compute gradient - radial component
      fac1 = qcontrib*ssij/nd
      fac = fac1*(dij-d0ij)/dij
      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 (il1.ne.il .or. jl1.ne.jl) then
        fac = fac1*(dij1-d0ij1)/dij1
        do k=1,3
          ddqij = (c(k,il1)-c(k,jl1))*fac
          if (il1.ne.il) then
            dxqwol(k,il)=dxqwol(k,il)+ddqij
          else
            dqwol(k,il)=dqwol(k,il)+ddqij
          endif
          if (jl1.ne.jl) then
            dxqwol(k,jl)=dxqwol(k,jl)-ddqij
          else
            dqwol(k,jl)=dqwol(k,jl)-ddqij
          endif
        enddo
      endif
c      return
c Orientational contributions
      rij=1.0d0/dij1
      eom1=qcontrib*(om1-om10)
      eom2=qcontrib*(om2-om20)
      eom12=qcontrib*(om12-om120)
      do k=1,3
        dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k))
        dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k))
      enddo
      do k=1,3
        ddqij=eom1*dcosom1(k)+eom2*dcosom2(k)
        aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
     &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
        aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
     &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
        dqwol(k,il)=dqwol(k,il)-ddqij-aux1
        dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2
        dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1
c     &            +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
c     &            +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
        dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2
c     &            +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
c     &            +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
      enddo
      return
      end