X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=3aa5fd6d1bc0cd38cc0af2f5eb5cc3533d687d51;hb=4c8b3525a22636eb28709e3f43c91127f57f8637;hp=918c761197e505f1304f369146566eb863868989;hpb=d0ddccf5d840c8ca0e2c721c3859d026d3a37aa8;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 918c761..3aa5fd6 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -2724,8 +2724,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 @@ -2743,6 +2744,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) @@ -2770,30 +2772,8 @@ 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 @@ -2801,24 +2781,7 @@ 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 -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 @@ -2829,23 +2792,53 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore rescore(k,irec)=0.0d0 goto 1301 endif - read (ientin,*,end=1401) rescore_tmp + if (read2sigma) then + read (ientin,*,end=1401) rescore2_tmp,rescore_tmp, + & idomain(k,irec) + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) + rescore2(k,irec)=0.5d0*(rescore2_tmp+0.5d0) + 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 + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) + endif 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) + + 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' + + unres_pdb=.false. + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + endif +c +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 + + + if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 ! right? without parallel. @@ -2854,6 +2847,9 @@ c do i = 1,nres ! alternative for bounds as used to set initial values i 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 + ii=ii+1 c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", @@ -2882,7 +2878,16 @@ 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 + 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 @@ -2893,25 +2898,24 @@ c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) 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) 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 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 @@ -2930,6 +2934,10 @@ c close (ientin) 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. + 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) @@ -2966,6 +2974,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), @@ -2991,12 +3003,16 @@ 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) @@ -3011,10 +3027,8 @@ 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