X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=f1cc85ba77f3b38605697abcae657259ecdcde8c;hb=300ba71bc245fd18664d59de0972b36fa1bb50b7;hp=47850c2e392e350aa5985702faf08331f9723054;hpb=e11674c0a2345839beaea75df0a4c17f517e16c7;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 47850c2..f1cc85b 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -96,6 +96,7 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes @@ -372,7 +373,6 @@ C endif call reada(controlcard,"Q_NP",Q_np,0.1d0) usampl = index(controlcard,"USAMPL").gt.0 - mdpdb = index(controlcard,"MDPDB").gt.0 call reada(controlcard,"T_BATH",t_bath,300.0d0) call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) @@ -749,6 +749,10 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) + 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) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) @@ -778,11 +782,16 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(18)=scal14 weights(21)=wsccor endif + weights(25)=wdfa_dist + weights(26)=wdfa_tor + weights(27)=wdfa_nei + weights(28)=wdfa_beta if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -801,7 +810,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file)then if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', @@ -829,7 +842,8 @@ C 12/1/95 Added weight for the multi-body term WCORR if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 22 format (/'Energy-term weights (scaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -848,7 +862,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:", & temp0 @@ -860,12 +878,36 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + dyn_ss=(index(weightcard,'DYN_SS').gt.0) + do i=1,maxres + dyn_ss_mask(i)=.false. + enddo + 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 + if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss - write (iout,*) "EBR",ebr + write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth + write (iout,*)" HT",Ht print *,'indpdb=',indpdb,' pdbref=',pdbref endif if (indpdb.gt.0 .or. pdbref) then @@ -893,6 +935,9 @@ c print *,'Finished reading pdb data' call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "Before sideadd" + call intout if(me.eq.king.or..not.out1file) & write(iout,*)'Adding sidechains' maxsi=1000 @@ -909,7 +954,12 @@ c print *,'Finished reading pdb data' & i,' after ',nsi,' trials' endif enddo +C 10/03/12 Adam: Recalculate coordinates with new side chain positions + call chainbuild endif +C Following 2 lines for diagnostics; comment out if not needed +c write (iout,*) "After sideadd" +c call intout endif if (indpdb.eq.0) then C Read sequence if not taken from the pdb file. @@ -1007,6 +1057,21 @@ C 8/13/98 Set limits to generating the dihedral angles cd print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 + +C Bartek: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 + call init_dfa_vars + print*, 'init_dfa_vars finished!' + call read_dfa_info + print*, 'read_dfa_info finished!' + endif +C if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup @@ -1095,8 +1160,15 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_dist.gt.0) then call read_dist_constr - call hpb_partition endif + + + if (constr_homology.gt.0) then + call read_constr_homology + endif + + + if (nhpb.gt.0) call hpb_partition c write (iout,*) "After read_dist_constr nhpb",nhpb c call flush(iout) if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 @@ -1217,18 +1289,35 @@ C Generate distance constraints, if the PDB structure is to be regularized. write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) + if (dyn_ss) then + write(iout,*)"Running with dynamic disulfide-bond formation" + else 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) - if (me.eq.king.or..not.out1file) - & write (iout,'(2a,i3,3a,i3,a,3f10.3)') + write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') + 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. + enddo endif if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) @@ -1292,10 +1381,12 @@ C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) @@ -1306,7 +1397,8 @@ C Check whether the specified bridging residues are cystines. 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(fg_rank.eq.0) + & 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 @@ -2434,7 +2526,7 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) 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 +c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=nhpb+1 @@ -2524,6 +2616,343 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) return end c------------------------------------------------------------------------------- + + subroutine read_constr_homology + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.GEO' + include 'COMMON.INTERACT' +c +c For new homol impl +c + include '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./ +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 pdbfile,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) + 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 + + 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" + + call readpdb +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.gt.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 + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) + +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.gt.0.0d0) lim_odl=ii + if (constr_homology.gt.0) call homology_partition + if (constr_homology.gt.0) call init_int_table + write (iout,*) "homology_partition: lim_theta= ",lim_theta, + & "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 + write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d + 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 + +c ----------------------------------------------------------------- + return + end +c---------------------------------------------------------------------- + #ifdef WINIFL subroutine flush(iu) return @@ -2535,6 +2964,7 @@ c------------------------------------------------------------------------------- return end #endif + c------------------------------------------------------------------------------ subroutine copy_to_tmp(source) include "DIMENSIONS"