subroutine molread(*) C C Read molecular data. C implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' 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 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 '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 c write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) call card_concat(controlcard,.true.) 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) 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 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) 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 do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) 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 dhpb(nhpb)=ddjk forcon(nhpb)=wfrag_(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) endif else nhpb=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 ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo endif enddo do i=1,npair_ if (wpair_(i).gt.0.0d0) then 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) 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 ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) enddo enddo endif enddo 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 nhpb=nhpb+1 if (ibecarb(i).gt.0) then ihpb(i)=ihpb(i)+nres jhpb(i)=jhpb(i)+nres endif if (dhpb(nhpb).eq.0.0d0) & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) endif 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 call flush(iout) return end c====------------------------------------------------------------------- subroutine read_constr_homology include 'DIMENSIONS' include 'DIMENSIONS.ZSCOPT' #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.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 logical lprn /.true./ logical unres_pdb c c FP - Nov. 2014 Temporary specifications for new vars c double precision rescore_tmp,x12,y12,z12 double precision, dimension (max_template,maxres) :: rescore 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 readi(controlcard,"HOMOL_NSET",homol_nset,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 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 c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d 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)) open(ipdbin,file=pdbfile,status='old',err=33) goto 34 33 write (iout,'(a)') 'Error opening PDB file.' 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" c tpl_k_sigma_odl="template"//kic2//".sigma_odl" c tpl_k_sigma_dih="template"//kic2//".sigma_dih" c tpl_k_sigma_theta="template"//kic2//".sigma_theta" c tpl_k_sigma_d="template"//kic2//".sigma_d" 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') do irec=1,maxdim ! loop for reading res sim if (irec.eq.1) then rescore(k,irec)=0.0d0 goto 1301 endif 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) 1301 continue enddo 1401 continue close (ientin) c open (ientin,file=tpl_k_sigma_odl,status='old') c do irec=1,maxdim ! loop for reading sigma_odl c read (ientin,*,end=1401) i, j, c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?) c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose? c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) c enddo c 1401 continue c close (ientin) if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 ! right? without parallel. do j=i+2,nct ! right? c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology c do j=i+2,nres ! ibid c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above) c do j=i+2,nct ! ibid ii=ii+1 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 c c Attempt to replace dist(i,j) by its definition in ... c 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) odl(k,ii)=distal c c odl(k,ii)=dist(i,j) c write (iout,*) "dist(",i,j,") =",dist(i,j) c write (iout,*) "distal = ",distal c write (iout,*) "odl(",k,ii,") =",odl(k,ii) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), c & "rescore(",k,j,") =",rescore(k,j) c c Calculation of sigma from res sim c c if (odl(k,ii).le.6.0d0) then c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) c Other functional forms possible depending on odl(k,ii), eg. c if (odl(k,ii).le.dist_cut) then sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) else #ifdef OLDSIGMA sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif c Following expr replaced by a positive exp argument c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)* c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) endif c sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii) c c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?) c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity enddo c read (ientin,*) sigma_odl(k,ii) ! 1st variant enddo c lim_odl=ii c if (constr_homology.gt.0) call homology_partition 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 ! right? without parallel. c do i=1,nres ! alternative for bounds acc to readpdb? c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel. 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) ! right expression ? c 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 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) if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right? c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file enddo 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. 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) ! right expression ? 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) if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right? 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 close (ientin) 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) goto 1 ! right? 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 ? 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? 1 continue enddo endif close(ientin) enddo if (waga_dist.ne.0.0d0) lim_odl=ii 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,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology) enddo write (iout,*) "Dihedral angle restraints from templates" do i=nnt+3,lim_dih write (iout,'(i5,10(2f8.2,4x))') 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,10(2f8.2,4x))') 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 write(iout,'(i5,10(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----------------------------------------------------------------------