X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=6b59d69374adb84fd4bc34735694dcc46fa943dd;hb=96d29f8250d30733e8a1f624f918388ef4f4050f;hp=dd5536b5a22d4f6aecb252832046797479d6398b;hpb=be29024edb52876c35e47eee5659b84f8c0566f6;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index dd5536b..6b59d69 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -997,6 +997,20 @@ 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 + C Assign initial virtual bond lengths do i=2,nres vbld(i)=vbl @@ -2567,6 +2581,7 @@ c call flush(iout) call readi(controlcard,"NPAIR",npair_,0) call readi(controlcard,"NDIST",ndist_,0) call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) + call reada(controlcard,'DIST2_CUT',dist2_cut,9999.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) @@ -2726,8 +2741,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 +2761,17 @@ 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) + start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) + if(.not.read2sigma.and.start_from_model) then + write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA' + start_from_model=.false. + endif + if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON' + if(start_from_model .and. rest) then + write(iout,*) 'START_FROM_MODELS is OFF' + write(iout,*) 'remove restart keyword from input' + endif if (homol_nset.gt.1)then call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) @@ -2811,13 +2838,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 @@ -2841,8 +2868,8 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore & idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp - rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0) - rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0) + rescore(k,i_tmp)=rescore_tmp + rescore2(k,i_tmp)=rescore2_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -2853,27 +2880,21 @@ 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 - - if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then - + 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) + + + 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. @@ -2883,62 +2904,35 @@ 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 + 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))* ! sigma ~ rescore ~ error + 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))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & 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 + 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 -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 +2948,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) @@ -2968,7 +2963,7 @@ 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 +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)* @@ -2993,6 +2988,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 +2999,8 @@ 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 +c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/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 +3018,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 +3041,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