NMR log-exp from Adam

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

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.