NMR log-exp from Adam

diff --git a/source/unres/src_MD-M/COMMON.SBRIDGE b/source/unres/src_MD-M/COMMON.SBRIDGE
index 469add2..0df27ec 100644 (file)
--- a/source/unres/src_MD-M/COMMON.SBRIDGE
+++ b/source/unres/src_MD-M/COMMON.SBRIDGE
@@ -2,16 +2,24 @@
       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,xlscore,wboltzd
-      integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,irestr_type
+      double precision dhpb,dhpb1,forcon,fordepth,xlscore,wboltzd,
+     & dhpb_peak,dhpb1_peak,forcon_peak,fordepth_peak,scal_peak
+      integer ihpb,jhpb,nhpb,idssb,jdssb,ibecarb,ibecarb_peak,npeak,
+     & ipeak,
+     & irestr_type,irestr_type_peak,ihpb_peak,jhpb_peak,nhpb_peak
       common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
      & fordepth(maxdim),xlscore(maxdim),wboltzd, 
      & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),irestr_type(maxdim),
      & nhpb
+      common /NMRpeaks/ dhpb_peak(maxdim),dhpb1_peak(maxdim),
+     & forcon_peak(maxdim),fordepth_peak(maxdim),scal_peak,
+     & ihpb_peak(maxdim),jhpb_peak(maxdim),ibecarb_peak(maxdim),
+     & irestr_type_peak(maxdim),ipeak(2,maxdim),npeak,nhpb_peak
       double precision weidis
       common /restraints/ weidis
-      integer link_start,link_end
-      common /links_split/ link_start,link_end
+      integer link_start,link_end,link_start_peak,link_end_peak
+      common /links_split/ link_start,link_end,link_start_peak,
+     & link_end_peak
       double precision Ht,dyn_ssbond_ij,dtriss,atriss,btriss,ctriss
       logical dyn_ss,dyn_ss_mask
       common /dyn_ssbond/ dyn_ssbond_ij(maxres,maxres),

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index badacd2..332f1ee 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -5192,7 +5192,7 @@ C
       include 'COMMON.INTERACT'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
-      dimension ggg(3)
+      dimension ggg(3),ggg_peak(3,20)
       ehpb=0.0D0
       do i=1,3
        ggg(i)=0.0d0
@@ -5200,8 +5200,67 @@ C
 C      write (iout,*) ,"link_end",link_end,constr_dist
 cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 c      write(iout,*)'link_start=',link_start,' link_end=',link_end,
-c     &  " constr_dist",constr_dist
-      if (link_end.eq.0) return
+c     &  " constr_dist",constr_dist," link_start_peak",link_start_peak,
+c     &  " link_end_peak",link_end_peak
+      if (link_end.eq.0.and.link_end_peak.eq.0) return
+      do i=link_start_peak,link_end_peak
+        ehpb_peak=0.0d0
+c        print *,"i",i," link_end_peak",link_end_peak," ipeak",
+c     &   ipeak(1,i),ipeak(2,i)
+        do ip=ipeak(1,i),ipeak(2,i)
+          ii=ihpb_peak(ip)
+          jj=jhpb_peak(ip)
+          dd=dist(ii,jj)
+          iip=ip-ipeak(1,i)+1
+C iii and jjj point to the residues for which the distance is assigned.
+          if (ii.gt.nres) then
+            iii=ii-nres
+            jjj=jj-nres 
+          else
+            iii=ii
+            jjj=jj
+          endif
+          aux=rlornmr1(dd,dhpb_peak(ip),dhpb1_peak(ip),forcon_peak(ip))
+          aux=dexp(-scal_peak*aux)
+          ehpb_peak=ehpb_peak+aux
+          fac=rlornmr1prim(dd,dhpb_peak(ip),dhpb1_peak(ip),
+     &      forcon_peak(ip))*aux/dd
+          do j=1,3
+            ggg_peak(j,iip)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          if (energy_dec) write (iout,'(a6,3i5,6f10.3,i5)')
+     &      "edisL",i,ii,jj,dd,dhpb_peak(ip),dhpb1_peak(ip),
+     &      forcon_peak(ip),fordepth_peak(ip),ehpb_peak
+        enddo
+c        write (iout,*) "ehpb_peak",ehpb_peak," scal_peak",scal_peak
+        ehpb=ehpb-fordepth_peak(ipeak(1,i))*dlog(ehpb_peak)/scal_peak
+        do ip=ipeak(1,i),ipeak(2,i)
+          iip=ip-ipeak(1,i)+1
+          do j=1,3
+            ggg(j)=ggg_peak(j,iip)/ehpb_peak
+          enddo
+          ii=ihpb_peak(ip)
+          jj=jhpb_peak(ip)
+C iii and jjj point to the residues for which the distance is assigned.
+          if (ii.gt.nres) then
+            iii=ii-nres
+            jjj=jj-nres 
+          else
+            iii=ii
+            jjj=jj
+          endif
+          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
+          endif
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
+        enddo
+      enddo
       do i=link_start,link_end
 C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
 C CA-CA distance used in regularization of structure.

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index d789d82..11a8267 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -1526,3 +1526,25 @@ cd     &   " lim_dih",lim_dih
 #endif
       return
       end
+c------------------------------------------------------------------------------
+      subroutine NMRpeak_partition
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.SETUP'
+#ifdef MPI
+      call int_bounds(npeak,link_start_peak,link_end_peak)
+      write (iout,*) 'Processor',fg_rank,' CG group',kolor,
+     &  ' absolute rank',MyRank,
+     &  ' npeak',npeak,' link_start_peak=',link_start_peak,
+     &  ' link_end_peak',link_end_peak
+#else
+      link_start_peak=1
+      link_end_peak=npeak
+#endif
+      return
+      end

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index add54af..a94bb53 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -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)