X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fcontact.f;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fcontact.f;h=24b11d6311fced7cdf1f2e125016dd601e9cc45c;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/contact.f b/source/unres/src_MD-M-newcorr/contact.f new file mode 100644 index 0000000..24b11d6 --- /dev/null +++ b/source/unres/src_MD-M-newcorr/contact.f @@ -0,0 +1,195 @@ + subroutine contact(lprint,ncont,icont,co) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + include 'COMMON.NAMES' + real*8 facont /1.569D0/ ! facont = (2/(1-sqrt(1-1/4)))**(1/6) + integer ncont,icont(2,maxcont) + logical lprint + ncont=0 + kkk=3 + do i=nnt+kkk,nct + iti=iabs(itype(i)) + do j=nnt,i-kkk + itj=iabs(itype(j)) + if (ipot.ne.4) then +c rcomp=sigmaii(iti,itj)+1.0D0 + rcomp=facont*sigmaii(iti,itj) + else +c rcomp=sigma(iti,itj)+1.0D0 + rcomp=facont*sigma(iti,itj) + endif +c rcomp=6.5D0 +c print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j) + if (dist(nres+i,nres+j).lt.rcomp) then + ncont=ncont+1 + icont(1,ncont)=i + icont(2,ncont)=j + endif + enddo + enddo + if (lprint) then + write (iout,'(a)') 'Contact map:' + do i=1,ncont + i1=icont(1,i) + i2=icont(2,i) + it1=itype(i1) + it2=itype(i2) + write (iout,'(i3,2x,a,i4,2x,a,i4)') + & i,restyp(it1),i1,restyp(it2),i2 + enddo + endif + co = 0.0d0 + do i=1,ncont + co = co + dfloat(iabs(icont(1,i)-icont(2,i))) + enddo + co = co / (nres*ncont) + return + end +c---------------------------------------------------------------------------- + double precision function contact_fract(ncont,ncont_ref, + & icont,icont_ref) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont) + nmatch=0 +c print *,'ncont=',ncont,' ncont_ref=',ncont_ref +c write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref) +c write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref) +c write (iout,'(20i4)') (icont(1,i),i=1,ncont) +c write (iout,'(20i4)') (icont(2,i),i=1,ncont) + do i=1,ncont + do j=1,ncont_ref + if (icont(1,i).eq.icont_ref(1,j) .and. + & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1 + enddo + enddo +c print *,' nmatch=',nmatch +c contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref)) + contact_fract=dfloat(nmatch)/dfloat(ncont_ref) + return + end +c---------------------------------------------------------------------------- + double precision function contact_fract_nn(ncont,ncont_ref, + & icont,icont_ref) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont) + nmatch=0 +c print *,'ncont=',ncont,' ncont_ref=',ncont_ref +c write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref) +c write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref) +c write (iout,'(20i4)') (icont(1,i),i=1,ncont) +c write (iout,'(20i4)') (icont(2,i),i=1,ncont) + do i=1,ncont + do j=1,ncont_ref + if (icont(1,i).eq.icont_ref(1,j) .and. + & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1 + enddo + enddo +c print *,' nmatch=',nmatch +c contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref)) + contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont) + return + end +c---------------------------------------------------------------------------- + subroutine hairpin(lprint,nharp,iharp) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + include 'COMMON.NAMES' + integer ncont,icont(2,maxcont) + integer nharp,iharp(4,maxres/3) + logical lprint,not_done + real*8 rcomp /6.0d0/ + ncont=0 + kkk=0 +c print *,'nnt=',nnt,' nct=',nct + do i=nnt,nct-3 + do k=1,3 + c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1)) + enddo + do j=i+2,nct-1 + do k=1,3 + c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1)) + enddo + if (dist(2*nres+1,2*nres+2).lt.rcomp) then + ncont=ncont+1 + icont(1,ncont)=i + icont(2,ncont)=j + endif + enddo + enddo + if (lprint) then + write (iout,'(a)') 'PP contact map:' + do i=1,ncont + i1=icont(1,i) + i2=icont(2,i) + it1=itype(i1) + it2=itype(i2) + write (iout,'(i3,2x,a,i4,2x,a,i4)') + & i,restyp(it1),i1,restyp(it2),i2 + enddo + endif +c finding hairpins + nharp=0 + do i=1,ncont + i1=icont(1,i) + j1=icont(2,i) + if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then +c write (iout,*) "found turn at ",i1,j1 + ii1=i1 + jj1=j1 + not_done=.true. + do while (not_done) + i1=i1-1 + j1=j1+1 + do j=1,ncont + if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10 + enddo + not_done=.false. + 10 continue +c write (iout,*) i1,j1,not_done + enddo + i1=i1+1 + j1=j1-1 + if (j1-i1.gt.4) then + nharp=nharp+1 + iharp(1,nharp)=i1 + iharp(2,nharp)=j1 + iharp(3,nharp)=ii1 + iharp(4,nharp)=jj1 +c write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4) + endif + endif + enddo +c do i=1,nharp +c write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4) +c enddo + if (lprint) then + write (iout,*) "Hairpins:" + do i=1,nharp + i1=iharp(1,i) + j1=iharp(2,i) + ii1=iharp(3,i) + jj1=iharp(4,i) + write (iout,*) + write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1) + write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1) +c do k=jj1,j1,-1 +c write (iout,'(a,i3,$)') restyp(itype(k)),k +c enddo + enddo + endif + return + end +c---------------------------------------------------------------------------- +