double precision function qwolynes(ib,iref,iprot)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.HEADER'
      integer ib,iref,iprot
      integer i,j,jl,ilnres,jlnres,klnres,k,l,il,kl,nl,np,ip,kp
      integer nsep /3/
      double precision dist,qm
      double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
      double precision qcontrib
      external qcontrib
      logical lprn /.false./
      logical flag
      qq = 0.0d0
      nl=0 
      if (lprn) then
      call flush(iout)
      endif
      do il=nnt+nsep,nct
        if (itype(il).eq.ntyp1) cycle
        do jl=nnt,il-nsep
          if (itype(jl).eq.ntyp1) cycle
          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,iref,ib,iprot)
          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 
      qwolynes = qq/nl
      return 
      end
c--------------------------------------------------------------------------- 
      double precision function qcontrib(il,jl,il1,jl1,iref,ib,iprot)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.HEADER'
      include 'COMMON.CLASSES'
      integer i,j,k,il,jl,il1,jl1,nd,iref,ib,iprot
      double precision dist
      external dist
      double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
     &  fac1,ddave,ssij,ddqij
      logical lprn /.false./
      d0ij1=dsqrt(
     &           (cref(1,jl,iref,ib,iprot)-cref(1,il,iref,ib,iprot))**2+
     &           (cref(2,jl,iref,ib,iprot)-cref(2,il,iref,ib,iprot))**2+
     &           (cref(3,jl,iref,ib,iprot)-cref(3,il,iref,ib,iprot))**2)
      dij1=dist(il,jl)
      ddave=(dij1-d0ij1)**2
      nd=1
      if (jl1.ne.jl) then
        d0ij2=dsqrt((cref(1,jl1,iref,ib,iprot)-
     &          cref(1,il,iref,ib,iprot))**2+
     &         (cref(2,jl1,iref,ib,iprot)-cref(2,il,iref,ib,iprot))**2+
     &         (cref(3,jl1,iref,ib,iprot)-cref(3,il,iref,ib,iprot))**2)
        dij2=dist(il,jl1)
        ddave=ddave+(dij2-d0ij2)**2
        nd=nd+1
      endif
      if (il1.ne.il) then
        d0ij3=dsqrt((cref(1,jl,iref,ib,iprot)-
     &          cref(1,il1,iref,ib,iprot))**2+
     &         (cref(2,jl,iref,ib,iprot)-cref(2,il1,iref,ib,iprot))**2+
     &         (cref(3,jl,iref,ib,iprot)-cref(3,il1,iref,ib,iprot))**2)
        dij3=dist(il1,jl)
        ddave=ddave+(dij3-d0ij3)**2
        nd=nd+1
      endif
      if (il1.ne.il .and. jl1.ne.jl) then
        d0ij4=dsqrt((cref(1,jl1,iref,ib,iprot)-
     &          cref(1,il1,iref,ib,iprot))**2+
     &         (cref(2,jl1,iref,ib,iprot)-cref(2,il1,iref,ib,iprot))**2+
     &         (cref(3,jl1,iref,ib,iprot)-cref(3,il1,iref,ib,iprot))**2)
        dij4=dist(il1,jl1)
        ddave=ddave+(dij4-d0ij4)**2
        nd=nd+1
      endif
      ddave=ddave/nd
      if (lprn) then
        write (iout,*) "il",il," jl",jl,
     &  " itype",itype(il),itype(jl)," nd",nd
        write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
     &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
        call flush(iout)
      endif
c      ssij = (0.25d0*d0ij1)**2
      if (il.ne.il1 .and. jl.ne.jl1) then
c        ssij = 16.0d0/(d0ij1*d0ij4)
        ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij4)
      else
c        ssij = 16.0d0/(d0ij1*d0ij1)
c        ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij1)
        ssij = 1.0d0/(sigma2(iprot)*sigma2(iprot))
      endif
      qcontrib = dexp(-0.5d0*ddave*ssij)
      return
      end