subroutine Econstr_back_qlike
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 sigmaang/0.1d0/,sigmadih /0.1d0/,sigmasc /0.1d0/
c      double precision sigmaang/0.2d0/,sigmadih /0.4d0/,sigmasc /0.5d0/
      double precision auxvec(maxres),auxtab(3,maxres),
     & auxtab1(3,maxres),auxtabx(3,maxres)
      Uconst_back=0.0d0
      do i=1,nres
        dutheta(i)=0.0d0
        dugamma(i)=0.0d0
        do j=1,3
          duscdiff(j,i)=0.0d0
          duscdiffx(j,i)=0.0d0
        enddo
      enddo
c      write (iout,*) "Econstr_back_qlike",nfrag_back," iset",iset
      do i=1,nfrag_back
        if (wfrag_back(1,i,iset).gt.0.0d0) then
        ii = ifrag_back(2,i,iset)-ifrag_back(1,i,iset)
c        write (iout,*) "i",i," ifrag_back",ifrag_back(1,i,iset),
c     &     ifrag_back(2,i,iset)," ii",ii
c
c Deviations from theta angles
c
        utheta_i=0.0d0
        do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset)
          dtheta_i=theta(j)-thetaref(j)
          expthet=dexp(-0.5d0*dtheta_i*dtheta_i/sigmaang)
c          expthet=0.5d0*dtheta_i*dtheta_i
          utheta_i=utheta_i+expthet
          auxvec(j-2)=expthet*dtheta_i/sigmaang
c          auxvec(j-2)=dtheta_i
c          write (iout,*) "j",j," theta",theta(j)," thetaref",thetaref(j)
c          write (iout,*) "expthet",expthet
        enddo
        qloc(1,i)=1.0d0-utheta_i/(ii-1)
c        qloc(1,i,iset)=utheta_i/(ii-1)
c        utheta(i)=(qloc(1,i,iset)-qin_back(1,i,iset))**2
        utheta(i)=(qloc(1,i)-qin_back(1,i,iset))**2
c        utheta(i)=qloc(1,i,iset)
c        write (iout,*) "utheta",utheta(i)
        do j=ifrag_back(1,i,iset)+2,ifrag_back(2,i,iset)
c          dutheta(j-2)=dutheta(j-2)+wfrag_back(1,i,iset)*auxvec(j-2)
c     &                  /(ii-1)
          dutheta(j-2)=dutheta(j-2)+2*wfrag_back(1,i,iset)*auxvec(j-2)
     &                  *(qloc(1,i)-qin_back(1,i,iset))/(ii-1)
c          write (iout,*) i,j," dutheta",dutheta(j-2)
        enddo
        endif
c
c Deviations from gamma angles
c
        if (wfrag_back(2,i,iset).gt.0.0d0) then
        ugamma_i=0.0d0
        do j=ifrag_back(1,i,iset)+3,ifrag_back(2,i,iset)
          dgamma_i=pinorm(phi(j)-phiref(j))
c          write (iout,*) j,phi(j),phi(j)-phiref(j)
          expgam=dexp(-0.5d0*dgamma_i*dgamma_i/sigmadih) 
          ugamma_i=ugamma_i+expgam
          auxvec(j-3)=expgam*dgamma_i/sigmadih
        enddo
        qloc(2,i)=1.0d0-ugamma_i/(ii-2)
        ugamma(i)=(qloc(2,i)-qin_back(2,i,iset))**2
        do j=ifrag_back(1,i,iset)+3,ifrag_back(2,i,iset)
          dugamma(j-3)=dugamma(j-3)+2*wfrag_back(2,i,iset)*auxvec(j-3)*
     &        (qloc(2,i)-qin_back(2,i,iset))/(ii-2)
c          write (iout,*) i,j," dugamma",dugamma(j-3)
        enddo
        endif
c
c Deviations from local SC geometry
c
        if (wfrag_back(3,i,iset).gt.0.0d0) then
        usc_i=0.0d0
        do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1
c          write (iout,*) "uscdif j=",j
          if (itype(j).ne.10) then
          dxx=xxtab(j)-xxref(j)
          dyy=yytab(j)-yyref(j)
          dzz=zztab(j)-zzref(j)
          expsc=dexp(-0.5d0*(dxx*dxx+dyy*dyy+dzz*dzz)/sigmasc)
          usc_i=usc_i+expsc
          expsc=expsc/sigmasc
          do k=1,3
            auxtab1(k,j-1)=expsc*
     &       (dXX_C1tab(k,j)*dxx+dYY_C1tab(k,j)*dyy+dZZ_C1tab(k,j)*dzz)
            auxtab(k,j)=expsc*
     &       (dXX_Ctab(k,j)*dxx+dYY_Ctab(k,j)*dyy+dZZ_Ctab(k,j)*dzz)
            auxtabx(k,j)=expsc*
     &     (dXX_XYZtab(k,j)*dxx+dYY_XYZtab(k,j)*dyy+dZZ_XYZtab(k,j)*dzz)
          enddo
c          write (iout,'(i5,6f10.5)') j,xxtab(j),yytab(j),zztab(j),
c     &      xxref(j),yyref(j),zzref(j)
          endif
        enddo
        qloc(3,i)=1.0d0-usc_i/(ii-1)
        uscdiff(i)=(qloc(3,i)-qin_back(3,i,iset))**2
        do j=ifrag_back(1,i,iset)+1,ifrag_back(2,i,iset)-1
          if (itype(j).ne.10) then
          do k=1,3
            duscdiff(k,j-1)=duscdiff(k,j-1)+2*wfrag_back(3,i,iset)
     &       *auxtab1(k,j-1)*(qloc(3,i)-qin_back(3,i,iset))/(ii-1)
            duscdiff(k,j)=duscdiff(k,j)+2*wfrag_back(3,i,iset)
     &       *auxtab(k,j)*(qloc(3,i)-qin_back(3,i,iset))/(ii-1)
            duscdiffx(k,j)=duscdiffx(k,j)+2*wfrag_back(3,i,iset)
     &       *auxtabx(k,j)*(qloc(3,i)-qin_back(3,i,iset))/(ii-1)
          enddo
c          write (iout,*) i,j," duscdiff",(duscdiff(k,j),k=1,3)
          endif
        enddo
c        write (iout,*) i," uscdiff",uscdiff(i)
        endif
c
c Put together deviations from local geometry
c
c        Uconst_back=Uconst_back+
c     &    wfrag_back(3,i,iset)*uscdiff(i)
        Uconst_back=Uconst_back+wfrag_back(1,i,iset)*utheta(i)+
     &    wfrag_back(2,i,iset)*ugamma(i)+wfrag_back(3,i,iset)*uscdiff(i)
c        write(iout,*) "i",i," utheta",utheta(i)," ugamma",ugamma(i),
c     &   " uscdiff",uscdiff(i)," uconst_back",uconst_back
        utheta(i)=qloc(1,i)
        ugamma(i)=qloc(2,i)
        uscdiff(i)=qloc(3,i)
      enddo
      return
      end