X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fcontact.f;fp=source%2Fwham%2Fsrc-NEWSC%2Fcontact.f;h=5b05d577e0233ac4b3ef9eea30c720cebc049f21;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/contact.f b/source/wham/src-NEWSC/contact.f new file mode 100755 index 0000000..5b05d57 --- /dev/null +++ b/source/wham/src-NEWSC/contact.f @@ -0,0 +1,171 @@ + subroutine contact(lprint,ncont,icont,ist,ien) + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + include 'COMMON.NAMES' + include 'COMMON.CALC' + include 'COMMON.CONTPAR' + include 'COMMON.LOCAL' + integer ist,ien,kkk,iti,itj,itypi,itypj,i1,i2,it1,it2 + real*8 csc,dist + real*8 cscore(maxcont),omt1(maxcont),omt2(maxcont),omt12(maxcont), + & ddsc(maxcont),ddla(maxcont),ddlb(maxcont) + integer ncont,icont(2,maxcont) + real*8 u,v,a(3),b(3),dla,dlb + logical lprint + ncont=0 + kkk=3 + if (lprint) then + do i=1,nres + write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i), + & c(3,i),dc(1,nres+i),dc(2,nres+i),dc(3,nres+i), + & dc_norm(1,nres+i),dc_norm(2,nres+i),dc_norm(3,nres+i) + enddo + endif + 110 format (a,'(',i3,')',9f8.3) + do i=ist,ien-kkk + iti=itype(i) + do j=i+kkk,ien + itj=itype(j) + itypi=iti + itypj=itj + xj = c(1,nres+j)-c(1,nres+i) + yj = c(2,nres+j)-c(2,nres+i) + zj = c(3,nres+j)-c(3,nres+i) + dxi = dc_norm(1,nres+i) + dyi = dc_norm(2,nres+i) + dzi = dc_norm(3,nres+i) + dxj = dc_norm(1,nres+j) + dyj = dc_norm(2,nres+j) + dzj = dc_norm(3,nres+j) + do k=1,3 + a(k)=dc(k,nres+i) + b(k)=dc(k,nres+j) + enddo +c write (iout,*) (a(k),k=1,3),(b(k),k=1,3) + if (icomparfunc.eq.1) then + call contfunc(csc,iti,itj) + else if (icomparfunc.eq.2) then + call scdist(csc,iti,itj) + else if (icomparfunc.eq.3 .or. icomparfunc.eq.5) then + csc = dist(nres+i,nres+j) + else if (icomparfunc.eq.4) then + call odlodc(c(1,i),c(1,j),a,b,u,v,dla,dlb,csc) + else + write (*,*) "Error - Unknown sidechain contact function" + write (iout,*) "Error - Unknown sidechain contact function" + endif + if (csc.lt.sc_cutoff(iti,itj)) then +c write(iout,*) "i",i," j",j," dla",dla,dsc(iti), +c & " dlb",dlb,dsc(itj)," csc",csc,sc_cutoff(iti,itj), +c & dxi,dyi,dzi,dxi**2+dyi**2+dzi**2, +c & dxj,dyj,dzj,dxj**2+dyj**2+dzj**2,om1,om2,om12, +c & xj,yj,zj +c write(iout,*)'egb',itypi,itypj,chi1,chi2,chip1,chip2, +c & sig0ij,rij,rrij,om1,om2,om12,chiom1,chiom2,chiom12, +c & chipom1,chipom2,chipom12,sig,eps2rt,rij_shift,e2,evdw, +c & csc + ncont=ncont+1 + cscore(ncont)=csc + icont(1,ncont)=i + icont(2,ncont)=j + omt1(ncont)=om1 + omt2(ncont)=om2 + omt12(ncont)=om12 + ddsc(ncont)=1.0d0/rij + ddla(ncont)=dla + ddlb(ncont)=dlb + 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,5f8.3,3f10.5)') + & i,restyp(it1),i1,restyp(it2),i2,cscore(i), + & sc_cutoff(it1,it2),ddsc(i),ddla(i),ddlb(i), + & omt1(i),omt2(i),omt12(i) + enddo + endif + return + end +c---------------------------------------------------------------------------- + double precision function contact_fract(ncont,ncont_ref, + & icont,icont_ref) + implicit none + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + integer i,j,nmatch + 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------------------------------------------------------------------------------ + subroutine pept_cont(lprint,ncont,icont) + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + include 'COMMON.FFIELD' + include 'COMMON.NAMES' + integer ncont,icont(2,maxcont) + integer i,j,k,kkk,i1,i2,it1,it2 + logical lprint + real*8 dist + real*8 rcomp /5.5d0/ + ncont=0 + kkk=0 + print *,'Entering pept_cont: 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 + return + end