c Restraints from contact prediction
          dd=dist(ii,jj)
          if (constr_dist.eq.11) then
            ehpb=ehpb+fordepth(i)!**4.0d0
     &          *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i))
            fac=fordepth(i)!**4.0d0
     &          *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd
            if(energy_dec) write (iout,'(a6,2i5,6f10.3)') "edisl",ii,jj,
     &        dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),ehpb
          else if (constr_dist.eq.10) then
c AL 6//19/2018 cross-link restraints
            xdis = (dd/forcon(i))**2
            expon = dexp(-xx)
            aux=(dhpb(i)+dhpb1(i)*xx)*expon+fordepth(i)
            ehpb=ehpb-wboltzd*dlog(aux)
            fac=-wboltzd*dhpb1(i)*(1.0d0-xx)*dxpon/(aux*forcon(i)**2)
            if (energy_dec)write (iout,'(a6,2i5,6f10.3)')"edisl",ii,jj,
     &      dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),ehpb
          else
            if (dhpb1(i).gt.0.0d0) then
              ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
              fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
c             write (iout,*) "beta nmr",
c     &       dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
            else
              dd=dist(ii,jj)
              rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
              waga=forcon(i)
C Calculate the contribution to energy.
              ehpb=ehpb+waga*rdis*rdis
c             write (iout,*) "beta reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
              fac=waga*rdis/dd
            endif
          endif
          do j=1,3
            ggg(j)=fac*(c(j,jj)-c(j,ii))
          enddo
          do j=1,3
            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
          enddo
          do k=1,3
            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
          enddo