NMR log-exp from Adam
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index add54af..a94bb53 100644 (file)
@@ -53,6 +53,24 @@ C Print restraint information
      &  irestr_type(i)
         enddo
       endif
+      if (npeak.gt.0) then
+        write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
+     &  "The following",npeak,
+     &  " NMR peak restraints have been imposed:",
+     &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
+     &  "  score"," type"," ipeak"
+        do i=1,npeak
+          do j=ipeak(1,i),ipeak(2,i)
+            write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
+     &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
+     &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
+          enddo
+        enddo
+        write (iout,*) "The ipeak array"
+        do i=1,npeak
+          write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
+        enddo
+      endif
 #ifdef MPI
       endif
 #endif
@@ -998,11 +1016,14 @@ czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
           enddo
           call contact(.true.,ncont_ref,icont_ref,co)
         endif
-c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+        endif
+c        print *, "A TU"
+        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         call flush(iout)
         if (constr_dist.gt.0) call read_dist_constr
-c        write (iout,*) "After read_dist_constr nhpb",nhpb
+        write (iout,*) "After read_dist_constr nhpb",nhpb
         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
+        call NMRpeak_partition
         if(me.eq.king.or..not.out1file)
      &   write (iout,*) 'Contact order:',co
         if (pdbref) then
@@ -1018,7 +1039,6 @@ c        write (iout,*) "After read_dist_constr nhpb",nhpb
      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
         enddo
         endif
-      endif
       write (iout,*) "calling read_saxs_consrtr",nsaxs
       if (nsaxs.gt.0) call read_saxs_constr
  
@@ -2548,9 +2568,9 @@ c-------------------------------------------------------------------------------
       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
+      integer ifrag_(2,100),ipair_(2,1000)
+      double precision wfrag_(100),wpair_(1000)
+      character*5000 controlcard
       logical normalize,next
       integer restr_type
       double precision xlink(4,0:4) /
@@ -2566,6 +2586,10 @@ c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 c      call flush(iout)
       next=.true.
 
+      npeak=0
+      ipeak=0
+      nhpb_peak=0
       DO WHILE (next)
 
       call card_concat(controlcard)
@@ -2579,21 +2603,24 @@ c      call flush(iout)
       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
       if (restr_type.eq.10) 
      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
+      if (restr_type.eq.12) 
+     &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
       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
+      write (iout,*) "SCAL_PEAK",scal_peak
+      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
       if (nfrag_.gt.0) write (iout,*) 
      &   "Distance restraints as generated from reference structure"
       do i=1,nfrag_
@@ -2653,13 +2680,14 @@ c            write (iout,*) "j",j," k",k
         endif
         do j=ifrag_(1,ii),ifrag_(2,ii)
           do k=ifrag_(1,jj),ifrag_(2,jj)
+            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) 
+              forcon(nhpb)=wpair_(i) 
             else if (constr_dist.eq.2) then
               if (ddjk.le.dist_cut) then
                 nhpb=nhpb+1
@@ -2667,7 +2695,7 @@ c            write (iout,*) "j",j," k",k
                 ihpb(nhpb)=j
                 jhpb(nhpb)=k
                 dhpb(nhpb)=ddjk
-                forcon(nhpb)=wfrag_(i) 
+                forcon(nhpb)=wpair_(i) 
               endif
             else if (restr_type.eq.3) then
               nhpb=nhpb+1
@@ -2675,14 +2703,14 @@ c            write (iout,*) "j",j," k",k
               ihpb(nhpb)=j
               jhpb(nhpb)=k
               dhpb(nhpb)=ddjk
-              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+              forcon(nhpb)=wpair_(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 ",
+     &      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,f10.1)') "+dist.restr ",
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 #endif
           enddo
@@ -2692,9 +2720,42 @@ c            write (iout,*) "j",j," k",k
 c      print *,ndist_
       write (iout,*) "Distance restraints as read from input"
       do i=1,ndist_
-        if (restr_type.eq.11) then
+        if (restr_type.eq.12) then
+          read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
+     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
+     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
+     &    fordepth_peak(nhpb_peak+1),npeak
+c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
+c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
+c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
+c     &    fordepth_peak(nhpb_peak+1),npeak
+          if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
+     &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
+          nhpb_peak=nhpb_peak+1
+          irestr_type_peak(nhpb_peak)=12
+          if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i 
+          ipeak(2,npeak)=i
+#ifdef MPI
+          if (.not.out1file .or. me.eq.king)
+     &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
+     &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
+     &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
+     &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
+#else
+          write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+     &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
+     &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
+     &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
+     &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
+#endif
+          if (ibecarb_peak(nhpb_peak).gt.0) then
+            ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
+            jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
+          endif
+        else 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)
+     &     dhpb1(nhpb+1),ibecarb(nhpb+1),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
@@ -2713,7 +2774,7 @@ c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
             ihpb(nhpb)=ihpb(nhpb)+nres
             jhpb(nhpb)=jhpb(nhpb)+nres
           endif
-        else if (constr_dist.eq.10) then
+        else if (restr_type.eq.10) then
 c Cross-lonk Markov-like potential
           call card_concat(controlcard)
           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)