subroutine read_dist_constr
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.CONTROL'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      include 'COMMON.SBRIDGE'
      integer ifrag_(2,100),ipair_(2,100)
      double precision wfrag_(100),wpair_(100)
      character*500 controlcard
      logical lprn /.true./
      logical normalize,next
      integer restr_type
      double precision xlink(4,0:4) /
c           a          b       c     sigma
     &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
     &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
     &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
     &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
     &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
      write (iout,*) "Calling read_dist_constr"
c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
c      call flush(iout)
      next=.true.

      DO WHILE (next)

      call card_concat(controlcard)
      next = index(controlcard,"NEXT").gt.0
      call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
      write (iout,*) "restr_type",restr_type
      call readi(controlcard,"NFRAG",nfrag_,0)
      call readi(controlcard,"NFRAG",nfrag_,0)
      call readi(controlcard,"NPAIR",npair_,0)
      call readi(controlcard,"NDIST",ndist_,0)
      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
      if (restr_type.eq.10) 
     &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
      call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
      call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
      call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
      call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
      normalize = index(controlcard,"NORMALIZE").gt.0
      write (iout,*) "WBOLTZD",wboltzd
c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
c      write (iout,*) "IFRAG"
c      do i=1,nfrag_
c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
c      enddo
c      write (iout,*) "IPAIR"
c      do i=1,npair_
c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
c      enddo
      if (nfrag_.gt.0) then
     &  write (iout,*) 
     &   "Distance restraints as generated from reference structure"
        read(inp,'(a)') pdbfile
        write (iout,'(2a,1h.)') 'PDB data will be read from file ',
     &    pdbfile(:ilen(pdbfile))
        open(ipdbin,file=pdbfile,status='old',err=33)
        goto 34
  33    write (iout,'(a)') 'Error opening PDB file.'
        return1
  34    continue
        do i=1,nres
          itype_pdb(i)=itype(i)
        enddo
        call readpdb(.true.)
        do i=1,nres
          iaux=itype_pdb(i)
          itype_pdb(i)=itype(i)
          itype(i)=iaux
        enddo
        close (ipdbin)
      endif
      do i=1,nfrag_
        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
        if (ifrag_(2,i).gt.nstart_sup+nsup-1)
     &    ifrag_(2,i)=nstart_sup+nsup-1
c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
c        call flush(iout)
        if (wfrag_(i).eq.0.0d0) cycle
        do j=ifrag_(1,i),ifrag_(2,i)-1
          do k=j+1,ifrag_(2,i)
c            write (iout,*) "j",j," k",k
            ddjk=dist(j,k)
            if (restr_type.eq.1) then
              nhpb=nhpb+1
              irestr_type(nhpb)=1
              ihpb(nhpb)=j
              jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
              forcon(nhpb)=wfrag_(i) 
            else if (constr_dist.eq.2) then
              if (ddjk.le.dist_cut) then
                nhpb=nhpb+1
                irestr_type(nhpb)=1
                ihpb(nhpb)=j
                jhpb(nhpb)=k
                dhpb(nhpb)=ddjk
                forcon(nhpb)=wfrag_(i) 
              endif
            else if (restr_type.eq.3) then
              nhpb=nhpb+1
              irestr_type(nhpb)=1
              ihpb(nhpb)=j
              jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
            endif
#ifdef MPI
            if (.not.out1file .or. me.eq.king) 
     &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
          enddo
        enddo
      enddo
      do i=1,npair_
        if (wpair_(i).eq.0.0d0) cycle
        ii = ipair_(1,i)
        jj = ipair_(2,i)
        if (ii.gt.jj) then
          itemp=ii
          ii=jj
          jj=itemp
        endif
        do j=ifrag_(1,ii),ifrag_(2,ii)
          do k=ifrag_(1,jj),ifrag_(2,jj)
            if (restr_type.eq.1) then
              nhpb=nhpb+1
              irestr_type(nhpb)=1
              ihpb(nhpb)=j
              jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
              forcon(nhpb)=wfrag_(i) 
            else if (constr_dist.eq.2) then
              if (ddjk.le.dist_cut) then
                nhpb=nhpb+1
                irestr_type(nhpb)=1
                ihpb(nhpb)=j
                jhpb(nhpb)=k
                dhpb(nhpb)=ddjk
                forcon(nhpb)=wfrag_(i) 
              endif
            else if (restr_type.eq.3) then
              nhpb=nhpb+1
              irestr_type(nhpb)=1
              ihpb(nhpb)=j
              jhpb(nhpb)=k
              dhpb(nhpb)=ddjk
              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
            endif
#ifdef MPI
            if (.not.out1file .or. me.eq.king)
     &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#else
            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
     &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
#endif
          enddo
        enddo
      enddo 

c      print *,ndist_
      write (iout,*) "Distance restraints as read from input"
      do i=1,ndist_
        if (restr_type.eq.11) then
          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
     &     dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
          if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
          nhpb=nhpb+1
          irestr_type(nhpb)=11
#ifdef MPI
          if (.not.out1file .or. me.eq.king)
     &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
     &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
#else
          write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
     &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
#endif
          if (ibecarb(nhpb).gt.0) then
            ihpb(nhpb)=ihpb(nhpb)+nres
            jhpb(nhpb)=jhpb(nhpb)+nres
          endif
        else if (constr_dist.eq.10) then
c Cross-lonk Markov-like potential
          call card_concat(controlcard)
          call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
          call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
          ibecarb(nhpb+1)=0
          if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
          if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
          if (index(controlcard,"ZL").gt.0) then
            link_type=1
          else if (index(controlcard,"ADH").gt.0) then
            link_type=2
          else if (index(controlcard,"PDH").gt.0) then
            link_type=3
          else if (index(controlcard,"DSS").gt.0) then
            link_type=4
          else
            link_type=0
          endif
          call reada(controlcard,"AXLINK",dhpb(nhpb+1),
     &       xlink(1,link_type))
          call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
     &       xlink(2,link_type))
          call reada(controlcard,"CXLINK",fordepth(nhpb+1),
     &       xlink(3,link_type))
          call reada(controlcard,"SIGMA",forcon(nhpb+1),
     &       xlink(4,link_type))
          call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
          if (forcon(nhpb+1).le.0.0d0 .or. 
     &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
          nhpb=nhpb+1
          irestr_type(nhpb)=10
          if (ibecarb(nhpb).gt.0) then
            ihpb(nhpb)=ihpb(nhpb)+nres
            jhpb(nhpb)=jhpb(nhpb)+nres
          endif
#ifdef MPI
          if (.not.out1file .or. me.eq.king)
     &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
     &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
     &     irestr_type(nhpb)
#else
          write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
     &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
     &     irestr_type(nhpb)
#endif
        else
C        print *,"in else"
          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
     &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
          if (forcon(nhpb+1).gt.0.0d0) then
          nhpb=nhpb+1
          if (dhpb1(nhpb).eq.0.0d0) then
            irestr_type(nhpb)=1
          else
            irestr_type(nhpb)=2
          endif
          if (ibecarb(nhpb).gt.0) then
            ihpb(nhpb)=ihpb(nhpb)+nres
            jhpb(nhpb)=jhpb(nhpb)+nres
          endif
          if (dhpb(nhpb).eq.0.0d0)
     &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
        endif
#ifdef MPI
          if (.not.out1file .or. me.eq.king)
     &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#else
          write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
#endif
        endif
C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
C        if (forcon(nhpb+1).gt.0.0d0) then
C          nhpb=nhpb+1
C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
      enddo

      ENDDO ! next

      fordepthmax=0.0d0
      if (normalize) then
        do i=nss+1,nhpb
          if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax) 
     &      fordepthmax=fordepth(i)
        enddo
        do i=nss+1,nhpb
          if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
        enddo
      endif
      if (nhpb.gt.nss)  then
        write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
     &  "The following",nhpb-nss,
     &  " distance restraints have been imposed:",
     &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
     &  "  score"," type"
        do i=nss+1,nhpb
          write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
     &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
     &  irestr_type(i)
        enddo
      endif
      write (iout,*) "Calling HPB_PARTINION"
      call hpb_partition
      call flush(iout)
      return
      end