Adam's new constr_dist single chain
[unres.git] / source / unres / src_MD / readrtns.F
index 6b4dcd0..38f8e89 100644 (file)
@@ -41,13 +41,18 @@ C Print restraint information
 #ifdef MPI
       if (.not. out1file .or. me.eq.king) then
 #endif
-      if (nhpb.gt.nss) 
-     &write (iout,'(a,i5,a)') "The following",nhpb-nss,
-     & " distance constraints have been imposed"
-      do i=nss+1,nhpb
-        write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
-     &     ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
-      enddo
+      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
 #ifdef MPI
       endif
 #endif
@@ -2575,19 +2580,40 @@ 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 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
+c      print *, "WCHODZE" 
+      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_
@@ -2597,46 +2623,38 @@ c      write (iout,*) "IPAIR"
 c      do i=1,npair_
 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 c      enddo
-      if (.not.refstr .and. nfrag.gt.0) then
-        write (iout,*) 
-     &  "ERROR: no reference structure to compute distance restraints"
-        write (iout,*)
-     &  "Restraints must be specified explicitly (NDIST=number)"
-        stop 
-      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 
-      call flush(iout)
+      if (nfrag_.gt.0) write (iout,*) 
+     &   "Distance restraints as generated from reference structure"
       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)
-        call flush(iout)
-        if (wfrag_(i).gt.0.0d0) then
+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 (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
@@ -2644,18 +2662,17 @@ c            write (iout,*) "j",j," k",k
             endif
 #ifdef MPI
             if (.not.out1file .or. me.eq.king) 
-     &      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)
 #else
-            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)
 #endif
           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
@@ -2665,59 +2682,161 @@ c            write (iout,*) "j",j," k",k
         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)
+            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.constr ",
+     &      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.constr ",
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           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
+          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
-          if (ibecarb(i).gt.0) then
-            ihpb(i)=ihpb(i)+nres
-            jhpb(i)=jhpb(i)+nres
+          irestr_type(nhpb)=10
+          if (ibecarb(nhpb).gt.0) then
+            ihpb(nhpb)=ihpb(nhpb)+nres
+            jhpb(nhpb)=jhpb(nhpb)+nres
           endif
-          if (dhpb(nhpb).eq.0.0d0) 
+#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
-      enddo
 #ifdef MPI
-      if (.not.out1file .or. me.eq.king) then
+          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
-      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
+        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
-      call flush(iout)
-#ifdef MPI
+
+      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
-#endif
       return
       end
 c-------------------------------------------------------------------------------