+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