implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.FREE'
-#ifdef MPI
- include 'mpif.h'
-#endif
include 'COMMON.CONTROL'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
integer ifrag_(2,100),ipair_(2,100)
double precision wfrag_(100),wpair_(100)
character*500 controlcard
- logical lprn /.true./
+ logical normalize,next
+ integer restr_type
+ double precision xlink(4,0:4) /
+c a b c sigma
+ & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
+ & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
+ & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
+ & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
+ & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
write (iout,*) "Calling read_dist_constr"
-C write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
-C call flush(iout)
- write(iout,*) "TU sie wywalam?"
- call card_concat(controlcard,.false.)
- write (iout,*) controlcard
- call flush(iout)
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
+ next=.true.
+
+ DO WHILE (next)
+
+ call card_concat(controlcard,.true.)
+ next = index(controlcard,"NEXT").gt.0
+ call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
+ write (iout,*) "restr_type",restr_type
+ call readi(controlcard,"NFRAG",nfrag_,0)
call readi(controlcard,"NFRAG",nfrag_,0)
call readi(controlcard,"NPAIR",npair_,0)
call readi(controlcard,"NDIST",ndist_,0)
call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ if (restr_type.eq.10)
+ & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
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)
- 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
- call flush(iout)
+ 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
+ if (nfrag_.gt.0) then
+ nres0=nres
+ read(inp,'(a)') pdbfile
+ write (iout,*)
+ & "Distance restraints will be constructed from structure ",pdbfile
+ open(ipdbin,file=pdbfile,status='old',err=11)
+ call readpdb(.true.)
+ nres=nres0
+ close(ipdbin)
+ endif
do i=1,nfrag_
if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
if (ifrag_(2,i).gt.nstart_sup+nsup-1)
& ifrag_(2,i)=nstart_sup+nsup-1
c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
- call flush(iout)
- if (wfrag_(i).gt.0.0d0) then
+c call flush(iout)
+ if (wfrag_(i).eq.0.0d0) cycle
do j=ifrag_(1,i),ifrag_(2,i)-1
do k=j+1,ifrag_(2,i)
- write (iout,*) "j",j," k",k
+c write (iout,*) "j",j," k",k
ddjk=dist(j,k)
- if (constr_dist.eq.1) then
+ if (restr_type.eq.1) then
nhpb=nhpb+1
+ irestr_type(nhpb)=1
ihpb(nhpb)=j
jhpb(nhpb)=k
dhpb(nhpb)=ddjk
else if (constr_dist.eq.2) then
if (ddjk.le.dist_cut) then
nhpb=nhpb+1
+ irestr_type(nhpb)=1
ihpb(nhpb)=j
jhpb(nhpb)=k
dhpb(nhpb)=ddjk
forcon(nhpb)=wfrag_(i)
endif
- else
+ else if (restr_type.eq.3) then
nhpb=nhpb+1
+ irestr_type(nhpb)=1
ihpb(nhpb)=j
jhpb(nhpb)=k
dhpb(nhpb)=ddjk
forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
endif
- if (lprn)
- & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
enddo
enddo
- endif
enddo
do i=1,npair_
- if (wpair_(i).gt.0.0d0) then
+ if (wpair_(i).eq.0.0d0) cycle
ii = ipair_(1,i)
jj = ipair_(2,i)
if (ii.gt.jj) then
endif
do j=ifrag_(1,ii),ifrag_(2,ii)
do k=ifrag_(1,jj),ifrag_(2,jj)
- nhpb=nhpb+1
- ihpb(nhpb)=j
- jhpb(nhpb)=k
- forcon(nhpb)=wpair_(i)
- dhpb(nhpb)=dist(j,k)
- write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ 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)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else if (restr_type.eq.3) then
+ nhpb=nhpb+1
+ irestr_type(nhpb)=1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
& nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
enddo
enddo
- endif
enddo
+
+c print *,ndist_
+ write (iout,*) "Distance restraints as read from input"
do i=1,ndist_
- if (constr_dist.eq.11) then
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
- & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
- fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
-C write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
-C & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ 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)
+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
+ irestr_type(nhpb)=11
+ write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ else if (constr_dist.eq.10) then
+c Cross-lonk Markov-like potential
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
+ call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
+ ibecarb(nhpb+1)=0
+ if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
+ if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
+ if (index(controlcard,"ZL").gt.0) then
+ link_type=1
+ else if (index(controlcard,"ADH").gt.0) then
+ link_type=2
+ else if (index(controlcard,"PDH").gt.0) then
+ link_type=3
+ else if (index(controlcard,"DSS").gt.0) then
+ link_type=4
+ else
+ link_type=0
+ endif
+ call reada(controlcard,"AXLINK",dhpb(nhpb+1),
+ & xlink(1,link_type))
+ call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
+ & xlink(2,link_type))
+ call reada(controlcard,"CXLINK",fordepth(nhpb+1),
+ & xlink(3,link_type))
+ call reada(controlcard,"SIGMA",forcon(nhpb+1),
+ & xlink(4,link_type))
+ call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
+c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
+c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
+ if (forcon(nhpb+1).le.0.0d0 .or.
+ & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
+ nhpb=nhpb+1
+ irestr_type(nhpb)=10
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
+ endif
+ write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
+ & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
+ & irestr_type(nhpb)
else
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
- endif
- if (forcon(nhpb+1).gt.0.0d0) then
+C print *,"in else"
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
+ & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
nhpb=nhpb+1
- if (ibecarb(i).gt.0) then
- ihpb(i)=ihpb(i)+nres
- jhpb(i)=jhpb(i)+nres
+ if (dhpb1(nhpb).eq.0.0d0) then
+ irestr_type(nhpb)=1
+ else
+ irestr_type(nhpb)=2
+ endif
+ if (ibecarb(nhpb).gt.0) then
+ ihpb(nhpb)=ihpb(nhpb)+nres
+ jhpb(nhpb)=jhpb(nhpb)+nres
endif
if (dhpb(nhpb).eq.0.0d0)
& dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
-C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
- write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
- & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
endif
-C endif
+ write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
+ endif
+C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+C if (forcon(nhpb+1).gt.0.0d0) then
+C nhpb=nhpb+1
+C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
enddo
+
+ ENDDO ! next
+
+ fordepthmax=0.0d0
+ if (normalize) then
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
+ & fordepthmax=fordepth(i)
+ enddo
+ do i=nss+1,nhpb
+ if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
+ enddo
+ endif
+ if (nhpb.gt.nss) then
+ write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
+ & "The following",nhpb-nss,
+ & " distance restraints have been imposed:",
+ & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
+ & " score"," type"
+ do i=nss+1,nhpb
+ write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
+ & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
+ & irestr_type(i)
+ enddo
+ endif
+ write (iout,*)
call hpb_partition
call flush(iout)
return
+ 11 write (iout,*)"read_dist_restr: error reading reference structure"
+ stop
end
projects / unres.git / blobdiff