X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc%2Freadrtns.F;h=fac2508c7a638464a9ad2f473198e1f7c4cfd998;hb=2c3a2cc1f2c409fa34495825d2c9febd53eb630e;hp=8e63ff8964d67fdfbfd34bdecfc9859bfd59a3e0;hpb=b0b9ebc44c9760d964bf33cc8c3a45b2e477be58;p=unres.git diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index 8e63ff8..fac2508 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -24,18 +24,21 @@ C call card_concat(controlcard) call readi(controlcard,'NRES',nres,0) + write (iout,*) "NRES",NRES call readi(controlcard,'RESCALE',rescale_mode,2) call readi(controlcard,'PDBOUT',outpdb,0) call readi(controlcard,'MOL2OUT',outmol2,0) refstr=(index(controlcard,'REFSTR').gt.0) write (iout,*) "REFSTR",refstr pdbref=(index(controlcard,'PDBREF').gt.0) + dyn_ss=(index(controlcard,'DYN_SS').gt.0) iscode=index(controlcard,'ONE_LETTER') tree=(index(controlcard,'MAKE_TREE').gt.0) min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) - call readi(controlcard,'NCUT',ncut,1) + call readi(controlcard,'NCUT',ncut,0) + call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) @@ -44,7 +47,8 @@ C lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) - call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) + if (ncut.gt.0) + & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 @@ -63,7 +67,21 @@ C call readi(controlcard,'CONSTR_DIST',constr_dist,0) write (iout,*) "with_dihed_constr ",with_dihed_constr, & " CONSTR_DIST",constr_dist + + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) + write (iout,*) "with_homology_constr ",with_dihed_constr, + & " CONSTR_HOMOLOGY",constr_homology + print_homology_restraints= + & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0 + print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0 + print_homology_models= + & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0 + +#ifdef AIX + call flush_(iout) +#else call flush(iout) +#endif if (min_var) iopt=1 return end @@ -98,6 +116,7 @@ C integer itype_pdb(maxres) logical seq_comp integer i,j + write (iout,*) " MOLREAD: NRES",NRES C C Body C @@ -126,9 +145,48 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'WSCCOR',wsccor,1.0D0) + call reada(weightcard,"D0CM",d0cm,3.78d0) + call reada(weightcard,"AKCM",akcm,15.1d0) + call reada(weightcard,"AKTH",akth,11.0d0) + call reada(weightcard,"AKCT",akct,12.0d0) + call reada(weightcard,"V1SS",v1ss,-1.08d0) + call reada(weightcard,"V2SS",v2ss,7.61d0) + call reada(weightcard,"V3SS",v3ss,13.7d0) + call reada(weightcard,"EBR",ebr,-5.50D0) +C Bartek + call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) + call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) + call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) + call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + endif + write (iout,'(/a)') "Disulfide bridge parameters:" + write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr + write (iout,'(a,f10.2)') 'S-S depth: ',ss_depth + write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm + write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct + write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss, + & ' v3ss:',v3ss + write (iout,'(2(a,f10.2))') 'ht:',ht,' eps:', eps(1,1) if (wcorr4.gt.0.0d0) wcorr=wcorr4 weights(1)=wsc weights(2)=wscp @@ -148,9 +206,13 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(16)=wvdwpp weights(17)=wbond weights(18)=scal14 + weights(22)=wdfa_dist + weights(23)=wdfa_tor + weights(24)=wdfa_nei + weights(25)=wdfa_beta write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3, - & wturn4,wturn6,wsccor + & wturn4,wturn6,wsccor,wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -169,7 +231,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ & 'WTURN6= ',f10.6,' (turns, 6th order)'/ - & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)') + & 'WSCCOR= ',f10.6,' (SC-backbone torsional correalations)'/ + & 'WDFAD= ',f10.6,' (DFA distance)'/ + & 'WDFAT= ',f10.6,' (DFA torsional)'/ + & 'WDFAN= ',f10.6,' (DFA neighbors)'/ + & 'WDFAB= ',f10.6,' (DFA beta)'/) if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', & 'between contact pairs of peptide groups' @@ -192,7 +258,11 @@ C 12/1/95 Added weight for the multi-body term WCORR bad(i,2)=scalscp*bad(i,2) enddo +#ifdef AIX + call flush_(iout) +#else call flush(iout) +#endif print *,'indpdb=',indpdb,' pdbref=',pdbref C Read sequence if not taken from the pdb file. @@ -205,6 +275,19 @@ 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 print *,nres print '(20i4)',(itype(i),i=1,nres) @@ -251,6 +334,7 @@ C Convert sequence to numeric code endif endif + nnt=1 nct=nres print *,'NNT=',NNT,' NCT=',NCT @@ -260,6 +344,47 @@ C Convert sequence to numeric code if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend nres0=nres + +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 +#ifdef DEBUG + write (iout,*) "Calling init_dfa_vars" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + call init_dfa_vars +#ifdef DEBUG + write (iout,*) 'init_dfa_vars finished!' +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + call read_dfa_info +#ifdef DEBUG + write (iout,*) 'read_dfa_info finished!' +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + endif + + if (constr_homology.gt.0) then + call read_constr_homology(print_homology_restraints) + endif + c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile @@ -327,7 +452,7 @@ c endif enddo enddo endif - call contact(.true.,ncont_ref,icont_ref) + call contact(print_contact_map,ncont_ref,icont_ref) endif c Read distance restraints if (constr_dist.gt.0) then @@ -373,7 +498,7 @@ C Read information about disulfide bridges. integer i,j C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) - print *,'ns=',ns + write(iout,*)'ns=',ns C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then @@ -431,6 +556,22 @@ c forcon(i)=fbr enddo endif endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + 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 write(iout,*) i,iss(i),dyn_ss_mask(iss(i)),"ATU" + enddo + endif + print *, "Leaving brigde read" return end c---------------------------------------------------------------------------- @@ -656,7 +797,11 @@ c write (iout,'(a)') controlcard do i=1,npair_ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) enddo +#ifdef AIX + call flush_(iout) +#else call flush(iout) +#endif if (.not.refstr .and. nfrag_.gt.0) then write (iout,*) & "ERROR: no reference structure to compute distance restraints" @@ -669,13 +814,17 @@ c write (iout,'(a)') controlcard & " but distance restraints between pairs requested" stop endif +#ifdef AIX + call flush_(iout) +#else call flush(iout) +#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) +c 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) @@ -731,8 +880,14 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) endif enddo 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 nhpb=nhpb+1 if (ibecarb(i).gt.0) then @@ -744,9 +899,473 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) endif enddo 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 enddo +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + return + end + +c====------------------------------------------------------------------- + subroutine read_constr_homology(lprn) + + include 'DIMENSIONS' +#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) + integer ilen + external ilen + logical lprn + 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 +#ifdef DEBUG + write (iout,*) "BEGIN READ HOMOLOGY INFO" +#ifdef AIX + call flush_(iout) +#else call flush(iout) +#endif +#endif + call card_concat(controlcard) + 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) + if (homol_nset.gt.1)then + call readi(controlcard,"ISET",iset,1) + call card_concat(controlcard) + read(controlcard,*) (waga_homology(i),i=1,homol_nset) + else + iset=1 + waga_homology(1)=1.0 + endif +c +#ifdef DEBUG + write(iout,*) "read_constr_homology iset",iset + write(iout,*) "waga_homology(",iset,")",waga_homology(iset) +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif +cd write (iout,*) "nnt",nnt," nct",nct +cd call flush(iout) + + + lim_odl=0 + lim_dih=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 + + write (iout,*) "CONSTR_HOMOLOGY",constr_homology + do k=1,constr_homology + + read(inp,'(a)') pdbfile +c write (iout,*) "k ",k," pdbfile ",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 + write (iout,*) "READ HOMOLOGY INFO" + write (iout,*) "read_constr_homology x: after reading pdb file" + write (iout,*) "waga_homology(",iset,")",waga_homology(iset) + write (iout,*) "waga_dist",waga_dist + write (iout,*) "waga_angle",waga_angle + write (iout,*) "waga_theta",waga_theta + write (iout,*) "waga_d",waga_d + write (iout,*) "dist_cut",dist_cut +#endif +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + +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 + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? + 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 + 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----------------------------------------------------------------------