X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=38f8e8958c99b2d2c9882583518377c715214589;hb=8568fc505b3e0d3f39901e7cf9de933603c71011;hp=6b4dcd0acabf84cd9ff05adc848244b50368966a;hpb=94b8e262bc9a6abf9c4a779950b4aaff41acd251;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 6b4dcd0..38f8e89 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -41,13 +41,18 @@ C Print restraint information #ifdef MPI if (.not. out1file .or. me.eq.king) then #endif - if (nhpb.gt.nss) - &write (iout,'(a,i5,a)') "The following",nhpb-nss, - & " distance constraints have been imposed" - do i=nss+1,nhpb - write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i), - & ibecarb(i),dhpb(i),dhpb1(i),forcon(i) - enddo + 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 #ifdef MPI endif #endif @@ -2575,19 +2580,40 @@ c------------------------------------------------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard -c write (iout,*) "Calling read_dist_constr" + 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 +c print *, "WCHODZE" + write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) + next=.true. + + DO WHILE (next) + call card_concat(controlcard) + 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) + 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_ @@ -2597,46 +2623,38 @@ c write (iout,*) "IPAIR" c do i=1,npair_ c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) c enddo - if (.not.refstr .and. nfrag.gt.0) then - write (iout,*) - & "ERROR: no reference structure to compute distance restraints" - write (iout,*) - & "Restraints must be specified explicitly (NDIST=number)" - stop - endif - if (nfrag.lt.2 .and. npair.gt.0) then - write (iout,*) "ERROR: Less than 2 fragments specified", - & " but distance restraints between pairs requested" - stop - endif - call flush(iout) + if (nfrag_.gt.0) write (iout,*) + & "Distance restraints as generated from reference structure" 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) c write (iout,*) "j",j," k",k ddjk=dist(j,k) - if (constr_dist.eq.1) then - nhpb=nhpb+1 - ihpb(nhpb)=j - jhpb(nhpb)=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)=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 + else if (restr_type.eq.3) then nhpb=nhpb+1 + irestr_type(nhpb)=1 ihpb(nhpb)=j jhpb(nhpb)=k dhpb(nhpb)=ddjk @@ -2644,18 +2662,17 @@ c write (iout,*) "j",j," k",k endif #ifdef MPI if (.not.out1file .or. me.eq.king) - & 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) #else - 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) #endif 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 @@ -2665,59 +2682,161 @@ c write (iout,*) "j",j," k",k 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) + 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 #ifdef MPI if (.not.out1file .or. me.eq.king) - & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #else - write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif 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) - else - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), - & ibecarb(i),forcon(nhpb+1) - endif - if (forcon(nhpb+1).gt.0.0d0) then + 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 +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & 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) +#else + 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) +#endif + 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) + 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 - if (ibecarb(i).gt.0) then - ihpb(i)=ihpb(i)+nres - jhpb(i)=jhpb(i)+nres + irestr_type(nhpb)=10 + if (ibecarb(nhpb).gt.0) then + ihpb(nhpb)=ihpb(nhpb)+nres + jhpb(nhpb)=jhpb(nhpb)+nres endif - if (dhpb(nhpb).eq.0.0d0) +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & 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 + 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) +#endif + else +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 (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)) endif - enddo #ifdef MPI - if (.not.out1file .or. me.eq.king) then + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb) #endif - do i=1,nhpb - if (constr_dist.eq.11) then - write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ", - & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i), - & fordepth(i) - else - write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", - & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) - endif + 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 - call flush(iout) -#ifdef MPI + + 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 -#endif return end c-------------------------------------------------------------------------------