subroutine molread(*) C C Read molecular data. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' c include 'include_unres/COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.TORCNSTR' include 'COMMON.CONTROL' character*4 sequence(maxres) integer rescode double precision x(maxvar) character*320 controlcard,ucase dimension itype_pdb(maxres) logical seq_comp call card_concat(controlcard,.true.) call reada(controlcard,'SCAL14',scal14,0.4d0) call reada(controlcard,'SCALSCP',scalscp,1.0d0) call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0) call reada(controlcard,'DELT_CORR',delt_corr,0.5d0) C Bartek call reada(controlcard,'WDFAD',wdfa_dist,0.0d0) call reada(controlcard,'WDFAT',wdfa_tor,0.0d0) call reada(controlcard,'WDFAN',wdfa_nei,0.0d0) call reada(controlcard,'WDFAB',wdfa_beta,0.0d0) write (iout,*) "wdfa_dist",wdfa_dist," wdfa_tor",wdfa_tor, & " wdfa_nei",wdfa_nei," wdfa_beta",wdfa_beta r0_corr=cutoff_corr-delt_corr call readi(controlcard,"NRES",nres,0) iscode=index(controlcard,"ONE_LETTER") if (nres.le.0) then write (iout,*) "Error: no residues in molecule" return1 endif if (nres.gt.maxres) then write (iout,*) "Error: too many residues",nres,maxres endif write(iout,*) 'nres=',nres C Read sequence of the protein if (iscode.gt.0) then read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres) else read (inp,'(20(1x,a3))') (sequence(i),i=1,nres) endif C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo if (itype(2).eq.10.and.itype(1).eq.ntyp1) then write (iout,*) & "Glycine is the first full residue, initial dummy deleted" do i=1,nres itype(i)=itype(i+1) enddo nres=nres-1 endif if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then write (iout,*) & "Glycine is the last full residue, terminal dummy deleted" nres=nres-1 endif write (iout,*) "Numeric code:" write (iout,'(20i4)') (itype(i),i=1,nres) do i=1,nres-1 #ifdef PROCOR if (itype(i).eq.21 .or. itype(i+1).eq.21) then #else if (itype(i).eq.21) then #endif itel(i)=0 #ifdef PROCOR else if (itype(i+1).ne.20) then #else else if (itype(i).ne.20) then #endif itel(i)=1 else itel(i)=2 endif enddo call read_bridge if (with_dihed_constr) then read (inp,*) ndih_constr if (ndih_constr.gt.0) then read (inp,*) ftors write (iout,*) 'FTORS',ftors read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) enddo do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo endif endif nnt=1 nct=nres if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 write(iout,*) 'NNT=',NNT,' NCT=',NCT C Juyong:READ init_vars C Initialize variables! C Juyong:READ read_info C READ fragment information!! C both routines should be in dfa.F file!! if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and. & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then write (iout,*) "Calling init_dfa_vars" call flush(iout) call init_dfa_vars write (iout,*) 'init_dfa_vars finished!' call flush(iout) call read_dfa_info write (iout,*) 'read_dfa_info finished!' call flush(iout) endif c Read distance restraints if (constr_dist.gt.0) then if (refstr) call read_ref_structure(*11) call read_dist_constr call hpb_partition endif if (constr_homology.gt.0) then c write (iout,*) "About to call read_constr_homology" c call flush(iout) call read_constr_homology c write (iout,*) "Exit read_constr_homology" c call flush(iout) if (indpdb.gt.0 .or. pdbref) then do i=1,2*nres do j=1,3 c(j,i)=crefjlee(j,i) cref(j,i)=crefjlee(j,i) enddo enddo endif #ifdef DEBUG write (iout,*) "Array C" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3), & (c(j,i+nres),j=1,3) enddo write (iout,*) "Array Cref" do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3), & (cref(j,i+nres),j=1,3) enddo #endif #ifdef DEBUG call int_from_cart1(.false.) call sc_loc_geom(.false.) do i=1,nres thetaref(i)=theta(i) phiref(i)=phi(i) write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i) enddo do i=1,nres-1 do j=1,3 dc(j,i)=c(j,i+1)-c(j,i) dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo #endif else homol_nset=0 endif call setup_var call init_int_table if (ns.gt.0) then write (iout,'(/a,i3,a)') 'The chain contains',ns, & ' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')', & dhpb(i),ebr,forcon(i) enddo endif write (iout,'(a)') return 11 stop "Error reading reference structure" end c----------------------------------------------------------------------------- logical function seq_comp(itypea,itypeb,length) implicit none integer length,itypea(length),itypeb(length) integer i do i=1,length if (itypea(i).ne.itypeb(i)) then seq_comp=.false. return endif enddo seq_comp=.true. return end c----------------------------------------------------------------------------- subroutine read_bridge C Read information about disulfide bridges. implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.INTERACT' include 'COMMON.NAMES' include 'COMMON.CHAIN' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) print *,'ns=',ns write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then write (iout,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(iss(i)),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') & 'Do you REALLY think that the residue ',restyp(iss(i)),i, & ' can form a disulfide bridge?!!!' stop endif enddo C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of C bridging residues. do i=1,nss do j=1,i-1 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then write (iout,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' write (*,'(a,i3,a)') 'Disulfide pair',i, & ' contains residues present in other pairs.' stop endif enddo do j=1,ns if (ihpb(i).eq.iss(j)) goto 10 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 10 continue do j=1,ns if (jhpb(i).eq.iss(j)) goto 20 enddo write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.' 20 continue dhpb(i)=dbr forcon(i)=fbr enddo do i=1,nss ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres enddo endif endif if (ns.gt.0.and.dyn_ss) then C /06/28/2013 Adasko:ns is number of Cysteins bonded also called half of C the bond do i=nss+1,nhpb C /06/28/2013 Adasko: nss number of full SS bonds ihpb(i-nss)=ihpb(i) jhpb(i-nss)=jhpb(i) forcon(i-nss)=forcon(i) dhpb(i-nss)=dhpb(i) enddo nhpb=nhpb-nss nss=0 call hpb_partition do i=1,ns dyn_ss_mask(iss(i))=.true. C /06/28/2013 Adasko: dyn_ss_mask which Cysteins can form disulfidebond c write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU" enddo endif return end c------------------------------------------------------------------------------ subroutine read_angles(kanal,iscor,energ,iprot,*) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' include 'COMMON.GEO' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.IOUNITS' character*80 lineh read(kanal,'(a80)',end=10,err=10) lineh read(lineh(:5),*,err=8) ic read(lineh(6:),*,err=8) energ goto 9 8 ic=1 print *,'error, assuming e=1d10',lineh energ=1d10 nss=0 9 continue read(lineh(18:),*,end=10,err=10) nss IF (NSS.LT.9) THEN read (lineh(20:),*,end=10,err=10) & (IHPB(I),JHPB(I),I=1,NSS),iscor ELSE read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8) read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I), & I=9,NSS),iscor ENDIF c print *,"energy",energ," iscor",iscor read (kanal,*,err=10,end=10) (theta(i),i=3,nres) read (kanal,*,err=10,end=10) (phi(i),i=4,nres) read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1) read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1) do i=1,nres theta(i)=deg2rad*theta(i) phi(i)=deg2rad*phi(i) alph(i)=deg2rad*alph(i) omeg(i)=deg2rad*omeg(i) enddo return 10 return1 end c------------------------------------------------------------------------------- subroutine read_dist_constr implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' include 'DIMENSIONS.FREE' include 'COMMON.CONTROL' 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 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) 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) 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 (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,1pe12.2)') "+dist.restr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo enddo do i=1,npair_ if (wpair_(i).eq.0.0d0) cycle ii = ipair_(1,i) jj = ipair_(2,i) if (ii.gt.jj) then itemp=ii ii=jj jj=itemp endif do j=ifrag_(1,ii),ifrag_(2,ii) do k=ifrag_(1,jj),ifrag_(2,jj) 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 enddo c print *,ndist_ write (iout,*) "Distance restraints as read from input" do i=1,ndist_ 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 (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 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 c====------------------------------------------------------------------- subroutine read_constr_homology 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' c c For new homol impl c include 'COMMON.VAR' c include 'include_unres/COMMON.VAR' c c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, c & dist_cut c common /przechowalnia/ odl_temp(maxres,maxres,max_template), c & sigma_odl_temp(maxres,maxres,max_template) 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 FP - Nov. 2014 Temporary specifications for new vars 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 Reading multiple PDB ref structures and calculation of retraints c not using pre-computed ones stored in files model_ki_{dist,angle} c FP (Nov., 2014) c ----------------------------------------------------------------- c c c Alternative: reading from input call card_concat(controlcard,.true.) call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new 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_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) call readi(controlcard,"IHSET",ihset,1) if (homol_nset.gt.1)then call card_concat(controlcard,.true.) read(controlcard,*) (waga_homology(i),i=1,homol_nset) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "iset homology_weight " c do i=1,homol_nset c write(iout,*) i,waga_homology(i) c enddo endif iset=mod(kolor,homol_nset)+1 else iset=1 waga_homology(1)=1.0 endif c write(iout,*) "waga_homology(",iset,")",waga_homology(iset) cd write (iout,*) "nnt",nnt," nct",nct cd call flush(iout) lim_odl=0 lim_dih=0 c c New c lim_theta=0 lim_xx=0 c c Reading HM global scores (prob not required) c do i = nnt,nct do k=1,constr_homology idomain(k,i)=0 enddo enddo c open (4,file="HMscore") c do k=1,constr_homology c read (4,*,end=521) hmscore_tmp c hmscore(k)=hmscore_tmp ! Another transformation can be used c write(*,*) "Model", k, ":", hmscore(k) c enddo c521 continue ii=0 do i = nnt,nct-2 do j=i+2,nct ii=ii+1 ii_in_use(ii)=0 enddo 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 c Next stament causes error upon compilation (?) c if(me.eq.king.or. .not. out1file) c write (iout,'(2a)') 'PDB data will be read from file ', c & pdbfile(:ilen(pdbfile)) write (iout,'(a,5x,a)') 'HOMOL: 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 c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c write(kic2,'(bz,i2.2)') k tpl_k_rescore="template"//kic2//".sco" unres_pdb=.false. call readpdb do i=1,2*nres do j=1,3 crefjlee(j,i)=c(j,i) enddo enddo #ifdef DEBUG do i=1,nres write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), & (crefjlee(j,i+nres),j=1,3) enddo #endif write (iout,*) "read_constr_homology: after reading pdb file" call flush(iout) 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)=cref(j,i) c write (iout,*) "c(",j,i,") =",c(j,i) enddo enddo c c From read_dist_constr (commented out 25/11/2014 <-> res sim) 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,nct ! loop for reading res sim if (read2sigma) then read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp, & idomain_tmp i_tmp=i_tmp+nnt-1 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 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue close (ientin) 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 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 c open (ientin,file=tpl_k_sigma_dih,status='old') c do irec=1,maxres-3 ! loop for reading sigma_dih c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_dih(k,i+nnt-1) c enddo c1402 continue c close (ientin) do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2), c & "rescore(",k,i-3,") =",rescore(k,i-3) sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ & rescore(k,i-2)+rescore(k,i-3))/4.0 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 c write (iout,*) "Raw sigmas for dihedral angle restraints" 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 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 endif if (waga_theta.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_theta,status='old') c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_theta(k,i+nnt-1) c enddo c1403 continue c close (ientin) do i = nnt+2,nct ! right? without parallel. c do i = i=1,nres ! alternative for bounds acc to readpdb? c do i=ithet_start,ithet_end ! with FG parallel. if (idomain(k,i).eq.0) then sigma_theta(k,i)=0.0 cycle endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,i-1,") =",rescore(k,i-1), c & "rescore(",k,i-2,") =",rescore(k,i-2) 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 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 if (waga_d.gt.0.0d0) then c open (ientin,file=tpl_k_sigma_d,status='old') c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity c & sigma_d(k,i+nnt-1) c enddo c1404 continue do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? c do i=loc_start,loc_end ! with FG parallel. if (itype(i).eq.10) cycle if (idomain(k,i).eq.0 ) then sigma_d(k,i)=0.0 cycle endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) 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 ? 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 enddo endif enddo 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 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, cd & "lim_xx=",lim_xx c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end c c Print restraints c if (.not.lprn) return cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write (iout,*) "Distance restraints from templates" do ii=1,lim_odl write(iout,'(3i5,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii), & ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" 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,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,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) enddo endif 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----------------------------------------------------------------------