c---------------------------------------------------------------------------------------------------------
      double precision function qwolynes(ilevel,jfrag,ib,iprot)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.PEPTCONT'
      include 'COMMON.CONTACTS1'
      include 'COMMON.HEADER'
      include 'COMMON.CLASSES'
      integer ilevel,jfrag,ib,iprot
      integer i,j,jl,k,l,il,kl,nl,np,ip,kp
      integer nsep /3/
      double precision dist
      double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
      logical lprn /.false./
      double precision sigm,x
      sigm(x)=0.25d0*x
      if (lprn) write (iout,*) "QWolyes: iprot",iprot," ib",ib,
     &      " jfrag",jfrag," ilevel",ilevel
      qq = 0.0d0
      if (ilevel.eq.0) then
        if (lprn) write (iout,*) "Q computed for whole molecule"
        nl=0
        do il=nnt+nsep,nct
          do jl=nnt,il-nsep
            dij=0.0d0
            dijCM=0.0d0
            d0ij=0.0d0
            d0ijCM=0.0d0
            qqij=0.0d0
            qqijCM=0.0d0
            nl=nl+1
            d0ij=dsqrt((cref(1,jl,iprot)-cref(1,il,iprot))**2+
     &                 (cref(2,jl,iprot)-cref(2,il,iprot))**2+
     &                 (cref(3,jl,iprot)-cref(3,il,iprot))**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,iprot)-cref(1,il+nres,iprot))**2+
     &               (cref(2,jl+nres,iprot)-cref(2,il+nres,iprot))**2+
     &               (cref(3,jl+nres,iprot)-cref(3,il+nres,iprot))**2)
              dijCM=dist(il+nres,jl+nres)
              qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
            endif
            qq = qq+qqij+qqijCM
c            if (lprn) then
c              write (iout,*) "il",il," jl",jl,
c     &         " itype",itype(il),itype(jl)
c              write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
c     &         " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
c            endif
          enddo
        enddo
        qq = qq/nl
        if (lprn) write (iout,*) "nl",nl," qq",qq
      else if (ilevel.eq.1) then
        if (lprn) then
          write (iout,*) "Level",ilevel," fragment",jfrag
          write (iout,*) "nlist_frag",nlist_frag(jfrag,ib,iprot),
     &      " list_frag",(list_frag(j,jfrag,ib,iprot),
     &      j=1,nlist_frag(jfrag,ib,iprot))
        endif
        nl=0
        do i=2,nlist_frag(jfrag,ib,iprot)
          do j=1,i-1
            il=list_frag(i,jfrag,ib,iprot)
            jl=list_frag(j,jfrag,ib,iprot)
            if (iabs(il-jl).gt.nsep) then
              dij=0.0d0
              dijCM=0.0d0
              d0ij=0.0d0
              d0ijCM=0.0d0
              qqij=0.0d0
              qqijCM=0.0d0
              nl=nl+1
              d0ij=dsqrt((cref(1,jl,iprot)-cref(1,il,iprot))**2+
     &                   (cref(2,jl,iprot)-cref(2,il,iprot))**2+
     &                   (cref(3,jl,iprot)-cref(3,il,iprot))**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,iprot)-cref(1,il+nres,iprot))**2+
     &                 (cref(2,jl+nres,iprot)-cref(2,il+nres,iprot))**2+
     &                 (cref(3,jl+nres,iprot)-cref(3,il+nres,iprot))**2)
                dijCM=dist(il+nres,jl+nres)
               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
              endif
              qq = qq+qqij+qqijCM
c              if (lprn) then
c                write (iout,*) "i",i," j",j," il",il," jl",jl,
c     &           " itype",itype(il),itype(jl)
c                write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
c     &           " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
c              endif
            endif
          enddo
        enddo
        qq = qq/nl
        if (lprn) write (iout,*) "nl",nl," qq",qq
      else if (ilevel.eq.2) then
        np=npiece(jfrag,ilevel,ib,iprot)
        nl=0
        if (lprn) then
          write (iout,*) "npiece",npiece(jfrag,ilevel,ib,iprot),
     &     " ipiece",(ipiece(i,jfrag,ilevel,ib,iprot),
     &       i=1,npiece(jfrag,ilevel,ib,iprot))
        endif
        do i=2,np
          ip=ipiece(i,jfrag,ilevel,ib,iprot)
          do j=1,nlist_frag(ip,ib,iprot) 
            il=list_frag(j,ip,ib,iprot)
            do k=1,i-1 
              kp=ipiece(k,jfrag,ilevel,ib,iprot)
              do l=1,nlist_frag(kp,ib,iprot)
                kl=list_frag(l,kp,ib,iprot)
                if (iabs(kl-il).gt.nsep) then 
                  nl=nl+1
                  dij=0.0d0
                  dijCM=0.0d0
                  d0ij=0.0d0
                  d0ijCM=0.0d0
                  qqij=0.0d0
                  qqijCM=0.0d0
                  d0ij=dsqrt((cref(1,kl,iprot)-cref(1,il,iprot))**2+
     &                       (cref(2,kl,iprot)-cref(2,il,iprot))**2+
     &                       (cref(3,kl,iprot)-cref(3,il,iprot))**2)
                  dij=dist(il,kl)
                  qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
                  if (itype(il).ne.10 .or. itype(kl).ne.10) then
                    nl=nl+1
                    d0ijCM=dsqrt(
     &                 (cref(1,kl+nres,iprot)-cref(1,il+nres,iprot))**2+
     &                 (cref(2,kl+nres,iprot)-cref(2,il+nres,iprot))**2+
     &                 (cref(3,kl+nres,iprot)-cref(3,il+nres,iprot))**2)
                    dijCM=dist(il+nres,kl+nres)
                    qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/
     &                (sigm(d0ijCM)))**2)
                  endif
                  qq = qq+qqij+qqijCM
c                  if (lprn) then
c                    write (iout,*) "i",i," j",j," k",k," l",l," il",il,
c     &                " kl",kl," itype",itype(il),itype(kl)
c                    write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",
c     &              d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
c                  endif
                endif
              enddo  ! l
            enddo    ! k
          enddo      ! j
        enddo        ! i
        qq = qq/nl
        if (lprn) write (iout,*) "nl",nl," qq",qq
      else
        write (iout,*)"Error: Q can be computed only for level 1 and 2."
      endif
      qwolynes=1.0d0-qq
      return
      end