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),
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
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.
& 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
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
& 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
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) /
c call flush(iout)
next=.true.
+ npeak=0
+ ipeak=0
+ nhpb_peak=0
+
DO WHILE (next)
call card_concat(controlcard)
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_
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
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
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
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
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)