From: Cezary Czaplewski Date: Fri, 31 Jul 2015 19:33:28 +0000 (+0200) Subject: restoring read2sigma code after wrong merge X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=ee883d438c43267225445df04be634ef1bd58186;p=unres.git restoring read2sigma code after wrong merge --- diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index dd5536b..ecbb19f 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -2726,8 +2726,9 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12 + 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 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -2745,6 +2746,7 @@ c Alternative: reading from input call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call readi(controlcard,"HOMOL_NSET",homol_nset,1) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) @@ -2811,13 +2813,13 @@ 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 + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + endif c c Distance restraints c @@ -2853,24 +2855,11 @@ c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif 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) + 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 + do i = nnt,nct-2 + do j=i+2,nct if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then @@ -2891,39 +2880,34 @@ c 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 + 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) ! 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 + 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 + 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 - + endif + 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) @@ -2935,10 +2919,7 @@ c & sigma_odl_temp(i,j,k) ! not inverse because of us l_homo(k,ii)=.false. endif enddo -c read (ientin,*) sigma_odl(k,ii) ! 1st variant enddo -c lim_odl=ii -c if (constr_homology.gt.0) call homology_partition lim_odl=ii endif c @@ -2954,10 +2935,11 @@ 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. + 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) @@ -2993,6 +2975,10 @@ 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), @@ -3000,7 +2986,7 @@ 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 + & rescore(k,i-2))/3.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)* @@ -3018,12 +3004,15 @@ c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use 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? + 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) @@ -3038,11 +3027,9 @@ 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 lim_xx=nct-nnt+1 endif - close(ientin) enddo c c remove distance restraints not used in any model from the list