X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fmolread_zs.F;h=c1a6469e0b61c214c20a8d8216b256f821e2d831;hb=ce876d8e106475d6cf7684c287eda53aeb15ca2f;hp=b11e3fbd90456edaee145f0b7b8a0661ac14845e;hpb=42251e287a50083c0337e4d78b357b370f832666;p=unres.git diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F index b11e3fb..c1a6469 100644 --- a/source/wham/src/molread_zs.F +++ b/source/wham/src/molread_zs.F @@ -375,81 +375,100 @@ 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 + 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,.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) - if (.not.refstr .and. nfrag_.gt.0) then + 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,*) - & "ERROR: no reference structure to compute distance restraints" - write (iout,*) - & "Restraints must be specified explicitly (NDIST=number)" - stop + & "Distance restraints will be constructed from structure ",pdbfile + open(ipdbin,file=pdbfile,status='old',err=11) + call readpdb(.true.) + nres=nres0 + close(ipdbin) 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) 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 - 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 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) endif - 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 @@ -459,40 +478,153 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) 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_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), - & ibecarb(i),forcon(nhpb+1) - 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 + 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 +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) + if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) 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 - do i=1,nhpb - write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ", - & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i) - 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 - - - c====------------------------------------------------------------------- subroutine read_constr_homology @@ -508,6 +640,7 @@ c====------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.INTERACT' + include 'COMMON.NAMES' include 'COMMON.HOMRESTR' c c For new homol impl @@ -528,7 +661,7 @@ c & sigma_odl_temp(maxres,maxres,max_template) logical lprn /.true./ integer ilen external ilen - logical unres_pdb + logical unres_pdb,liiflag c c FP - Nov. 2014 Temporary specifications for new vars c @@ -551,7 +684,7 @@ c Alternative: reading from input call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) - call readi(controlcard,"HOMOL_SET",homol_nset,1) + call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) call readi(controlcard,"IHSET",ihset,1) if (homol_nset.gt.1)then @@ -606,6 +739,10 @@ c521 continue enddo c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d + if (read_homol_frag) then + call read_klapaucjusz + else + do k=1,constr_homology read(inp,'(a)') pdbfile @@ -663,7 +800,7 @@ c c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore open (ientin,file=tpl_k_rescore,status='old') if (nnt.gt.1) rescore(k,1)=0.0d0 - do irec=nnt,maxdim ! loop for reading res sim + do irec=nnt,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, & idomain_tmp @@ -671,6 +808,9 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore idomain(k,i_tmp)=idomain_tmp rescore(k,i_tmp)=rescore_tmp rescore2(k,i_tmp)=rescore2_tmp + write(iout,'(a7,i5,2f10.5,i5)') "rescore", + & i_tmp,rescore2_tmp,rescore_tmp, + & idomain_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -770,7 +910,8 @@ c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2)*rescore(k,i-3) ! right expression ? c Instead of res sim other local measure of b/b str reliability possible - sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + if (sigma_dih(k,i).ne.0) + & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) enddo lim_dih=nct-nnt-2 @@ -802,14 +943,14 @@ c read (ientin,*) sigma_theta(k,i) ! 1st variant sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2))/3.0 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 - sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + if (sigma_theta(k,i).ne.0) + & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* c rescore(k,i-2) ! right expression ? c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) enddo endif - lim_theta=nct-nnt-1 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') @@ -836,14 +977,13 @@ c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i) sigma_d(k,i)=rescore(k,i) ! right expression ? - sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + if (sigma_d(k,i).ne.0) + & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) c read (ientin,*) sigma_d(k,i) ! 1st variant - if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? enddo - lim_xx=nct-nnt+1 endif enddo c @@ -852,26 +992,41 @@ c shift data in all arrays c if (waga_dist.ne.0.0d0) then ii=0 + liiflag=.true. do i=nnt,nct-2 do j=i+2,nct ii=ii+1 - if (ii_in_use(ii).eq.0) then - do ki=ii,lim_odl-1 - ires_homo(ki)=ires_homo(ki+1) - jres_homo(ki)=jres_homo(ki+1) - ii_in_use(ki)=ii_in_use(ki+1) + if (ii_in_use(ii).eq.0.and.liiflag) then + liiflag=.false. + iistart=ii + endif + if (ii_in_use(ii).ne.0.and..not.liiflag.or. + & .not.liiflag.and.ii.eq.lim_odl) then + if (ii.eq.lim_odl) then + iishift=ii-iistart+1 + else + iishift=ii-iistart + endif + liiflag=.true. + do ki=iistart,lim_odl-iishift + ires_homo(ki)=ires_homo(ki+iishift) + jres_homo(ki)=jres_homo(ki+iishift) + ii_in_use(ki)=ii_in_use(ki+iishift) do k=1,constr_homology - odl(k,ki)=odl(k,ki+1) - sigma_odl(k,ki)=sigma_odl(k,ki+1) - l_homo(k,ki)=l_homo(k,ki+1) + odl(k,ki)=odl(k,ki+iishift) + sigma_odl(k,ki)=sigma_odl(k,ki+iishift) + l_homo(k,ki)=l_homo(k,ki+iishift) enddo enddo - ii=ii-1 - lim_odl=lim_odl-1 + ii=ii-iishift + lim_odl=lim_odl-iishift endif enddo enddo endif + + endif ! .not. klapaucjusz + if (constr_homology.gt.0) call homology_partition if (constr_homology.gt.0) call init_int_table cd write (iout,*) "homology_partition: lim_theta= ",lim_theta, @@ -892,17 +1047,19 @@ cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" - do i=nnt+3,lim_dih - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*dih(ki,i), + do i=nnt+3,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + & (rad2deg*dih(ki,i), & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) enddo write (iout,*) "Virtual-bond angle restraints from templates" - do i=nnt+2,lim_theta - write (iout,'(i5,100(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + do i=nnt+2,nct + write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)), + & (rad2deg*thetatpl(ki,i), & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) enddo write (iout,*) "SC restraints from templates" - do i=nnt,lim_xx + do i=nnt,nct write(iout,'(i5,100(4f8.2,4x))') i, & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) @@ -912,3 +1069,264 @@ c ----------------------------------------------------------------- return end c---------------------------------------------------------------------- + subroutine read_klapaucjusz + + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.FREE' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + include 'COMMON.NAMES' + include 'COMMON.HOMRESTR' + character*256 fragfile + integer ninclust(maxclust),inclust(max_template,maxclust), + & nresclust(maxclust),iresclust(maxres,maxclust) + + character*2 kic2 + character*24 model_ki_dist, model_ki_angle + character*500 controlcard + integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp + integer idomain(max_template,maxres) + logical lprn /.true./ + integer ilen + external ilen + logical unres_pdb,liiflag +c +c + double precision rescore_tmp,x12,y12,z12,rescore2_tmp + double precision, dimension (max_template,maxres) :: rescore + double precision, dimension (max_template,maxres) :: rescore2 + character*24 tpl_k_rescore + +c +c For new homol impl +c + include 'COMMON.VAR' +c + double precision chomo(3,maxres2+2,max_template) + call getenv("FRAGFILE",fragfile) + open(ientin,file=fragfile,status="old",err=10) + read(ientin,*) constr_homology,nclust + l_homo = .false. + sigma_theta=0.0 + sigma_d=0.0 + sigma_dih=0.0 +c Read pdb files + do k=1,constr_homology + read(ientin,'(a)') pdbfile + write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', + & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) + stop + 34 continue + unres_pdb=.false. + call readpdb + do i=1,2*nres + do j=1,3 + chomo(j,i,k)=c(j,i) + enddo + enddo + do i=1,nres + rescore(k,i)=0.2d0 + rescore2(k,i)=1.0d0 + enddo + enddo +c Read clusters + do i=1,nclust + read(ientin,*) ninclust(i),nresclust(i) + read(ientin,*) (inclust(k,i),k=1,ninclust(i)) + read(ientin,*) (iresclust(k,i),k=1,nresclust(i)) + enddo +c +c Loop over clusters +c + do l=1,nclust + do ll = 1,ninclust(l) + + k = inclust(ll,l) + do i=1,nres + idomain(k,i)=0 + enddo + do i=1,nresclust(l) + if (nnt.gt.1) then + idomain(k,iresclust(i,l)+1) = 1 + else + idomain(k,iresclust(i,l)) = 1 + endif + enddo +c +c Distance restraints +c +c ... --> odl(k,ii) +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,k) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + if (waga_dist.ne.0.0d0) then + ii=0 + do i = nnt,nct-2 + do j=i+2,nct + + x12=c(1,i)-c(1,j) + y12=c(2,i)-c(2,j) + z12=c(3,i)-c(3,j) + distal=dsqrt(x12*x12+y12*y12+z12*z12) +c write (iout,*) k,i,j,distal,dist2_cut + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 + & .and. distal.le.dist2_cut ) then + + ii=ii+1 + ii_in_use(ii)=1 + l_homo(k,ii)=.true. + +c write (iout,*) "k",k +c write (iout,*) "i",i," j",j," constr_homology", +c & constr_homology + ires_homo(ii)=i + jres_homo(ii)=j + odl(k,ii)=distal + if (read2sigma) then + sigma_odl(k,ii)=0 + do ik=i,j + sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) + enddo + sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) + if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = + & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) + else + if (odl(k,ii).le.dist_cut) then + sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) + else +#ifdef OLDSIGMA + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) +#else + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) +#endif + endif + endif + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) + else + ii=ii+1 +c l_homo(k,ii)=.false. + endif + enddo + enddo + lim_odl=ii + endif +c +c Theta, dihedral and SC retraints +c + if (waga_angle.gt.0.0d0) then + do i = nnt+3,nct + if (idomain(k,i).eq.0) then +c sigma_dih(k,i)=0.0 + cycle + endif + dih(k,i)=phiref(i) + sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2)+rescore(k,i-3))/4.0 +c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i), +c & " sigma_dihed",sigma_dih(k,i) + if (sigma_dih(k,i).ne.0) + & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) + enddo + lim_dih=nct-nnt-2 + endif + + if (waga_theta.gt.0.0d0) then + do i = nnt+2,nct + if (idomain(k,i).eq.0) then +c sigma_theta(k,i)=0.0 + cycle + endif + thetatpl(k,i)=thetaref(i) + sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2))/3.0 + if (sigma_theta(k,i).ne.0) + & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + enddo + endif + + if (waga_d.gt.0.0d0) then + do i = nnt,nct + if (itype(i).eq.10) cycle + if (idomain(k,i).eq.0 ) then +c sigma_d(k,i)=0.0 + cycle + endif + xxtpl(k,i)=xxref(i) + yytpl(k,i)=yyref(i) + zztpl(k,i)=zzref(i) + sigma_d(k,i)=rescore(k,i) + if (sigma_d(k,i).ne.0) + & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 + enddo + endif + enddo ! l + enddo ! ll +c +c remove distance restraints not used in any model from the list +c shift data in all arrays +c + if (waga_dist.ne.0.0d0) then + ii=0 + liiflag=.true. + do i=nnt,nct-2 + do j=i+2,nct + ii=ii+1 + if (ii_in_use(ii).eq.0.and.liiflag) then + liiflag=.false. + iistart=ii + endif + if (ii_in_use(ii).ne.0.and..not.liiflag.or. + & .not.liiflag.and.ii.eq.lim_odl) then + if (ii.eq.lim_odl) then + iishift=ii-iistart+1 + else + iishift=ii-iistart + endif + liiflag=.true. + do ki=iistart,lim_odl-iishift + ires_homo(ki)=ires_homo(ki+iishift) + jres_homo(ki)=jres_homo(ki+iishift) + ii_in_use(ki)=ii_in_use(ki+iishift) + do k=1,constr_homology + odl(k,ki)=odl(k,ki+iishift) + sigma_odl(k,ki)=sigma_odl(k,ki+iishift) + l_homo(k,ki)=l_homo(k,ki+iishift) + enddo + enddo + ii=ii-iishift + lim_odl=lim_odl-iishift + endif + enddo + enddo + endif + + return + 10 stop "Error infragment file" + end +c----------------------------------------------------------------------