cluster_wham Adam's new constr_dist single chain
authorCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 22 Jun 2018 22:14:43 +0000 (00:14 +0200)
committerCezary Czaplewski <czarek@chem.univ.gda.pl>
Fri, 22 Jun 2018 22:14:43 +0000 (00:14 +0200)
source/cluster/wham/src/energy_p_new.F
source/cluster/wham/src/include_unres/COMMON.SBRIDGE
source/cluster/wham/src/readrtns.F

index 8006a41..a56af26 100644 (file)
       include 'COMMON.DERIV'
       include 'COMMON.VAR'
       include 'COMMON.INTERACT'
-      include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
+      include 'COMMON.IOUNITS'
       dimension ggg(3)
       ehpb=0.0D0
+      ggg=0.0d0
+C      write (iout,*) ,"link_end",link_end,constr_dist
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
-cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
+c      write(iout,*)'link_start=',link_start,' link_end=',link_end,
+c     &  " constr_dist",constr_dist
       if (link_end.eq.0) return
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
@@ -2899,101 +2902,81 @@ c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
 c     &    dhpb(i),dhpb1(i),forcon(i)
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
+C        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C     & iabs(itype(jjj)).eq.1) then
+cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
         if (.not.dyn_ss .and. i.le.nss) then
 C 15/02/13 CC dynamic SSbond - additional check
-        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
-          call ssbond_ene(iii,jjj,eij)
-          ehpb=ehpb+2*eij
+          if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+     &        iabs(itype(jjj)).eq.1) then
+           call ssbond_ene(iii,jjj,eij)
+           ehpb=ehpb+2*eij
+         endif
 cd          write (iout,*) "eij",eij
-        endif
-        else if (ii.gt.nres .and. jj.gt.nres) then
-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
-         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 !end dhpb1(i).gt.0
-         endif !end const_dist=11
-          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
-        else
+cd   &   ' waga=',waga,' fac=',fac
+!        else if (ii.gt.nres .and. jj.gt.nres) then
+        else 
 C Calculate the distance between the two points and its difference from the
 C target distance.
           dd=dist(ii,jj)
-C          write(iout,*) "after",dd
-          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
-C            ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i))
-C            fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd
-C            print *,ehpb,"tu?"
-C            write(iout,*) ehpb,"btu?",
-C     & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i)
-C          write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj,
-C     &    ehpb,fordepth(i),dd
-           else   
-          if (dhpb1(i).gt.0.0d0) then
-            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          if (irestr_type(i).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
+c            if (energy_dec) write (iout,'(a6,2i5,6f10.3,i5)')
+c     &        "edisL",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+c     &        ehpb,irestr_type(i)
+          else if (irestr_type(i).eq.10) then
+c AL 6//19/2018 cross-link restraints
+            xdis = 0.5d0*(dd/forcon(i))**2
+            expdis = dexp(-xdis)
+c            aux=(dhpb(i)+dhpb1(i)*xdis)*expdis+fordepth(i)
+            aux=(dhpb(i)+dhpb1(i)*xdis*xdis)*expdis+fordepth(i)
+c            write (iout,*)"HERE: xdis",xdis," expdis",expdis," aux",aux,
+c     &          " wboltzd",wboltzd
+            ehpb=ehpb-wboltzd*xlscore(i)*dlog(aux)
+c            fac=-wboltzd*(dhpb1(i)*(1.0d0-xdis)-dhpb(i))
+            fac=-wboltzd*xlscore(i)*(dhpb1(i)*(2.0d0-xdis)*xdis-dhpb(i))
+     &           *expdis/(aux*forcon(i)**2)
+c            if (energy_dec) write(iout,'(a6,2i5,6f10.3,i5)') 
+c     &        "edisX",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),fordepth(i),
+c     &        -wboltzd*xlscore(i)*dlog(aux),irestr_type(i)
+          else if (irestr_type(i).eq.2) then
+c Quartic restraints
+            ehpb=ehpb+forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+c            if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') 
+c     &      "edisQ",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+c     &      forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)),irestr_type(i)
             fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
-c            write (iout,*) "alph nmr",
-c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
           else
+c Quadratic restraints
             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,*) "alpha reg",dd,waga*rdis*rdis
+            ehpb=ehpb+0.5d0*waga*rdis*rdis
+c            if (energy_dec) write(iout,'(a6,2i5,5f10.3,i5)') 
+c     &      "edisS",ii,jj,dd,dhpb(i),dhpb1(i),forcon(i),
+c     &       0.5d0*waga*rdis*rdis,irestr_type(i)
 C
 C Evaluate gradient.
 C
             fac=waga*rdis/dd
           endif
-          endif
-cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-cd   &   ' waga=',waga,' fac=',fac
-            do j=1,3
-              ggg(j)=fac*(c(j,jj)-c(j,ii))
-            enddo
+c Calculate Cartesian gradient
+          do j=1,3
+            ggg(j)=fac*(c(j,jj)-c(j,ii))
+          enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
           if (iii.lt.ii) then
-          do j=1,3
-            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
-            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
-          enddo
+            do j=1,3
+              ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+              ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+            enddo
           endif
           do k=1,3
             ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
@@ -3001,7 +2984,6 @@ C Cartesian gradient in the SC vectors (ghpbx).
           enddo
         endif
       enddo
-      if (constr_dist.ne.11) ehpb=0.5D0*ehpb
       return
       end
 C--------------------------------------------------------------------------
index b68d0e3..d02241d 100644 (file)
@@ -2,17 +2,20 @@
       integer ns,nss,nfree,iss
       common /sbridge/ ss_depth,ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,
      & ns,nss,nfree,iss(maxss)
-      double precision dhpb,dhpb1,forcon,fordepth
-      integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb
+      double precision dhpb,dhpb1,forcon,fordepth,xlscore,wboltzd
+      integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,irestr_type
       common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
-     & fordepth(maxdim),
-     & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
+     & fordepth(maxdim),xlscore(maxdim),wboltzd,
+     & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),irestr_type(maxdim),
+     & nhpb
       double precision weidis
       common /restraints/ weidis
       integer link_start,link_end
       common /links_split/ link_start,link_end
-      double precision Ht,dyn_ssbond_ij
+      double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss
       logical dyn_ss,dyn_ss_mask
-      common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres),
-     &  idssb(maxdim),jdssb(maxdim),
-     &  Ht,dyn_ss,dyn_ss_mask(maxres)
+      common /dyn_ssbond/ dtriss,atriss,btriss,ctriss,Ht,
+     &  dyn_ssbond_ij(maxres,maxres),
+     &  idssb(maxdim),jdssb(maxdim)
+      common /dyn_ss_logic/
+     &  dyn_ss,dyn_ss_mask(maxres)
index 47ff41d..9a98fe1 100644 (file)
@@ -769,6 +769,9 @@ c-------------------------------------------------------------------------------
       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'
@@ -776,90 +779,101 @@ c-------------------------------------------------------------------------------
       integer ifrag_(2,100),ipair_(2,100)
       double precision wfrag_(100),wpair_(100)
       character*500 controlcard
-c      write (iout,*) "Calling read_dist_constr"
+      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
-      call card_concat(controlcard)
 c      call flush(iout)
-c      write (iout,'(a)') controlcard
+      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)
-      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
-      write (iout,*) "IFRAG"
-      do i=1,nfrag_
-        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
-      enddo
-      write (iout,*) "IPAIR"
-      do i=1,npair_
-        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
-      enddo
-#ifdef AIX
-      call flush_(iout)
-#else
-      call flush(iout)
-#endif
-      if (.not.refstr .and. nfrag_.gt.0) then
+      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
+        nres0=nres
+        read(inp,'(a)') pdbfile
         write (iout,*) 
-     &  "ERROR: no reference structure to compute distance restraints"
-        write (iout,*)
-     &  "Restraints must be specified explicitly (NDIST=number)"
-        stop 
+     & "Distance restraints will be constructed from structure ",pdbfile
+        open(ipdbin,file=pdbfile,status='old',err=11)
+        call readpdb(.true.)
+        nres=nres0
+        close(ipdbin)
       endif
-      if (nfrag_.lt.2 .and. npair_.gt.0) then 
-        write (iout,*) "ERROR: Less than 2 fragments specified",
-     &   " but distance restraints between pairs requested"
-        stop 
-      endif 
-#ifdef AIX
-      call flush_(iout)
-#else
-      call flush(iout)
-#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).gt.0.0d0) then
+        if (wfrag_(i).eq.0.0d0) cycle
         do j=ifrag_(1,i),ifrag_(2,i)-1
           do k=j+1,ifrag_(2,i)
-            write (iout,*) "j",j," k",k
+c            write (iout,*) "j",j," k",k
             ddjk=dist(j,k)
-            if (constr_dist.eq.1) then
-            nhpb=nhpb+1
-            ihpb(nhpb)=j
-            jhpb(nhpb)=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) 
+              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
+            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
-            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
           enddo
         enddo
-        endif
       enddo
       do i=1,npair_
-        if (wpair_(i).gt.0.0d0) then
+        if (wpair_(i).eq.0.0d0) cycle
         ii = ipair_(1,i)
         jj = ipair_(2,i)
         if (ii.gt.jj) then
@@ -869,54 +883,152 @@ c        call flush(iout)
         endif
         do j=ifrag_(1,ii),ifrag_(2,ii)
           do k=ifrag_(1,jj),ifrag_(2,jj)
-            nhpb=nhpb+1
-            ihpb(nhpb)=j
-            jhpb(nhpb)=k
-            forcon(nhpb)=wpair_(i)
-            dhpb(nhpb)=dist(j,k)
-            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+            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
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
           enddo
         enddo
-        endif
       enddo 
+
+c      print *,ndist_
+      write (iout,*) "Distance restraints as read from input"
       do i=1,ndist_
-       if (constr_dist.eq.11) then
-         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
-     &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
-         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
-       else
-        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
-     &     ibecarb(i),forcon(nhpb+1)
-       endif
-        if (forcon(nhpb+1).gt.0.0d0) then
+        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
-          if (ibecarb(i).gt.0) then
-            ihpb(i)=ihpb(i)+nres
-            jhpb(i)=jhpb(i)+nres
+          irestr_type(nhpb)=11
+          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)
+          if (ibecarb(nhpb).gt.0) then
+            ihpb(nhpb)=ihpb(nhpb)+nres
+            jhpb(nhpb)=jhpb(nhpb)+nres
           endif
-          if (dhpb(nhpb).eq.0.0d0) 
+        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
+          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
+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
+        write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+     &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
+        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
-      do i=1,nhpb
-       if (constr_dist.eq.11) then
-          write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
-     &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
-     &     fordepth(i)
-       else
-          write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
-     &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
-       endif
-      enddo
-#ifdef AIX
-      call flush_(iout)
-#else
+
+      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
+      call hpb_partition
       call flush(iout)
-#endif
       return
+   11 write (iout,*)"read_dist_restr: error reading reference structure"
+      stop
       end
-
 c====-------------------------------------------------------------------
       subroutine read_constr_homology(lprn)