source/maxlik/src_MD_T_maxlik-NEWCORR-PMF/natsimil.f

   1       subroutine natsimil(jcon,iii,iprot,lprn)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.COMPAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.COMPAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.VAR'
  11       include 'COMMON.CLASSES'
  12       include 'COMMON.VMCPAR'
  13       logical lprn
  14       integer jcon,iprot,ib
  15       integer i,iii,j,k,ik,kk,inat,lat,jfrag
  16       double precision fract
  17       double precision rmscalc,rms,rmsnat,qwolynes,gyrate,
  18      &  strand_signature
  19       character*4 liczba
  20       character*64 sbin
  21       ib = kbatch(jcon,iprot)
  22       if (lprn) then
  23         write(iout,*) "Protein",iprot," batch",ib," conformation",jcon,
  24      &    " nlevel",nlevel(iprot)
  25       endif
  26 c      if (lprn) then
  27 c        write (iout,*) "NATSIMIL: Complete reference structure"
  28 c        do i=1,nres
  29 c          write(iout,'(i4,3f10.5)') i,(cref(j,i,iprot),j=1,3)
  30 c        enddo
  31 c      endif
  32       inat=0
  33       ik = iscore(jcon,1,iprot)
  34       if (nlevel(iprot).gt.0) then
  35 c Level 1: local structure, Q value, and rmsd
  36         do j=1,nfrag(1,iprot)
  37 c Calculate the difference between the angles of the target structure and those of the native structure
  38           call angnorm(j,0,0,ang_cut1(j,ik,iprot),
  39      &    nu(inat+1,iii,iprot),fract,ib,iprot,.false.)
  40           nu(inat+2,iii,iprot)=qwolynes(1,j,ib,iprot)
  41           nu(inat+3,iii,iprot)=rmscalc(0,1,j,jcon,ib,iprot,.false.)
  42           nu(inat+3,iii,iprot)=nu(inat+3,iii,iprot)/
  43      &      len_frag(j,1,ib,iprot)
  44           nu(inat+4,iii,iprot)=gyrate(1,j,ib,iprot)
  45          nu(inat+5,iii,iprot)=1.0d0-strand_signature(j,ib,iprot,.false.)
  46           if (lprn) write (iout,'(a,i2,a,i2,a,f8.3,a,f10.5,a,f10.5,
  47      &      a,f10.5)')
  48      &      "i",1," j",j,
  49      &      " rms",nu(inat+3,iii,iprot),
  50      &      " diffang",nu(inat+1,iii,iprot)," Q",nu(inat+2,iii,iprot),
  51      &      " signature",nu(inat+5,iii,iprot)
  52 c Calculate fragment RMSD
  53           inat=inat+5
  54         enddo
  55 c Level 2: Q value and rmsd
  56         do j=1,nfrag(2,iprot)
  57           nu(inat+1,iii,iprot)=qwolynes(2,j,ib,iprot)
  58           nu(inat+2,iii,iprot)=rmscalc(0,2,j,jcon,ib,iprot,.false.)
  59           nu(inat+2,iii,iprot)=nu(inat+2,iii,iprot)/
  60      &      len_frag(j,2,ib,iprot)
  61           nu(inat+3,iii,iprot)=gyrate(2,j,ib,iprot)
  62           if (lprn) write (iout,'(a,i2,a,i2,a,f8.3,a,f10.5)')
  63      &       "i",2," j",j," rms",
  64      &       nu(inat+2,iii,iprot)," Q",nu(inat+1,iii,iprot)
  65           inat=inat+3
  66         enddo
  67 c Next levels: rmsd only
  68         do i=3,nlevel(iprot)
  69           do j=1,nfrag(i,iprot)
  70             nu(inat+1,iii,iprot)=rmscalc(0,i,j,jcon,ib,iprot,.false.)
  71             nu(inat+1,iii,iprot)=nu(inat+1,iii,iprot)
  72      &       /len_frag(j,i,ib,iprot)
  73             nu(inat+2,iii,iprot)=gyrate(i,j,ib,iprot)
  74             if (lprn) write (iout,'(a,i2,a,i2,a,f8.3)')
  75      &        "i",i," j",j," rms",nu(inat+1,iii,iprot)
  76             inat=inat+2
  77           enddo
  78         enddo
  79       endif
  80       call angnorm12(nu(inat+1,iii,iprot),iprot)
  81       nu(inat+2,iii,iprot)=qwolynes(0,1,ib,iprot)
  82       nu(inat+3,iii,iprot)=rmsnat(jcon,iprot)
  83       nu(inat+3,iii,iprot)=nu(inat+3,iii,iprot)/(nct-nnt+1)
  84       nu(inat+4,iii,iprot)=gyrate(0,1,ib,iprot)
  85       if (lprn) write (iout,'(a,f8.3,a,f10.5,a,f10.5)')
  86      & "angnorm_nat",nu(inat+1,iii,iprot),
  87      & " qnat",nu(inat+2,iii,iprot)," rmsnat",nu(inat+3,iii,iprot),
  88      & " rgy",nu(inat+4,iii,iprot)
  89       inat=inat+4
  90       if (lprn) write (iout,*) "inat",inat," natlike",natlike(iprot)
  91       RETURN
  92       END
  93 c---------------------------------------------------------------------------------------------------
  94       subroutine fragment_list(iprot)
  95       implicit none
  96       include 'DIMENSIONS'
  97       include 'DIMENSIONS.ZSCOPT'
  98       include 'DIMENSIONS.COMPAR'
  99       include 'COMMON.IOUNITS'
 100       include 'COMMON.ALLPROT'
 101       include 'COMMON.COMPAR'
 102       include 'COMMON.CLASSES'
 103       logical lprn /.false./
 104       integer ib,iprot
 105       integer i,ilevel,j,k,jfrag
 106       do ib=1,nclass(iprot)-1
 107       do jfrag=1,nfrag(1,iprot)
 108         nlist_frag(jfrag,ib,iprot)=0
 109         do i=1,npiece(jfrag,1,ib,iprot)
 110           if (lprn) write (iout,*) "Protein",iprot," batch",ib,
 111      &      " jfrag=",jfrag,
 112      &      "i=",i," fragment",ifrag(1,i,jfrag,ib,iprot),
 113      &      ifrag(2,i,jfrag,ib,iprot)
 114           do j=ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot)
 115             do k=1,nlist_frag(jfrag,ib,iprot)
 116               if (list_frag(k,jfrag,ib,iprot).eq.j) goto 10
 117             enddo
 118             nlist_frag(jfrag,ib,iprot)=nlist_frag(jfrag,ib,iprot)+1
 119             list_frag(nlist_frag(jfrag,ib,iprot),jfrag,ib,iprot)=j
 120           enddo
 121   10      continue
 122         enddo
 123       enddo
 124       if (lprn) then
 125       write (iout,*) "Fragment list for protein",iprot," batch",ib
 126       do j=1,nfrag(1,iprot)
 127         write (iout,*)"Fragment",j," list",(list_frag(k,j,ib,iprot),
 128      &   k=1,nlist_frag(j,ib,iprot))
 129       enddo
 130       endif
 131       enddo
 132       return
 133       end
 134 c-----------------------------------------------------------------------------
 135       subroutine find_near_pept_cont(iprot)
 136       implicit none
 137       include 'DIMENSIONS'
 138       include 'DIMENSIONS.ZSCOPT'
 139       include 'DIMENSIONS.COMPAR'
 140       include 'COMMON.IOUNITS'
 141       include 'COMMON.ALLPROT'
 142       include 'COMMON.COMPAR'
 143       include 'COMMON.CLASSES'
 144       include 'COMMON.PEPTCONT'
 145       include 'COMMON.INTERACT'
 146       include 'COMMON.NAMES'
 147       logical lprn /.false./
 148       integer ib,iprot
 149       integer i,ii,ilevel,j,lj,k,k1,k2,l
 150       do ib=1,nclass(iprot)-1
 151       do j=1,nfrag(1,iprot)
 152           if (lprn) write (iout,*) "Protein",iprot," batch",ib,
 153      &      " jfrag=",j," list",(list_frag(k,j,ib,iprot),
 154      &       k=1,nlist_frag(j,ib,iprot))
 155         do k=1,nlist_frag(j,ib,iprot)
 156           ii=0
 157           lj=list_frag(k,j,ib,iprot)
 158 c          write (iout,*) "j",j," k",k," lj",lj," sec",isec_ref(lj,iprot)
 159           if (isec_ref(lj,iprot).eq.1) then
 160             do l=1,ncont_pept_ref(iprot)
 161               k1 = icont_pept_ref(1,l,iprot)
 162               k2 = icont_pept_ref(2,l,iprot)
 163               if (k1.eq.lj .and. isec_ref(k2,iprot).eq.1
 164 c     &          .and. k2-lj.gt.3 ) then
 165      &          .and. iabs(k2-lj).gt.3 ) then
 166 c                write (iout,*) "k2",k2," sec",isec_ref(k2,iprot)
 167                 if (ii.lt.2) then
 168                   ii=ii+1
 169                   icont_pept_near(ii,lj,ib,iprot)=k2
 170                 else
 171                   write (iout,*)
 172      &             "Warning number of near contacts for beta fragment",
 173      &             j," residue",lj," exceeds 2, truncated."
 174                   goto 10
 175                 endif
 176               endif
 177               if (k2.eq.lj .and. isec_ref(k1,iprot).eq.1
 178 c     &          .and. k1-lj.gt.3) then
 179      &          .and. iabs(k1-lj).gt.3) then
 180 c                write (iout,*) "k1",k1," sec",isec_ref(k1,iprot)
 181                 if (ii.lt.2) then
 182                   ii=ii+1
 183                   icont_pept_near(ii,lj,ib,iprot)=k1
 184                 else
 185                   write (iout,*)
 186      &             "Warning number of near contacts for beta fragment",
 187      &             j," residue",lj," exceeds 2, truncated."
 188                   goto 10
 189                 endif
 190               endif
 191             enddo
 192           endif
 193           ncont_pept_near(lj,ib,iprot)=ii
 194         enddo
 195    10   continue
 196       enddo
 197       write (iout,*)
 198      &   "Near peptide group contacts for beta sheet fragments ",
 199      &   "protein",iprot," class",ib
 200       do j=1,nfrag(1,iprot)
 201         write (iout,*) "fragment",j
 202         do k=1,nlist_frag(j,ib,iprot)
 203           lj=list_frag(k,j,ib,iprot)
 204           write (iout,'(a4,1h(,i2,1h),2i5)') restyp(itype(lj)),lj,
 205      &      (icont_pept_near(ii,lj,ib,iprot),
 206      &       ii=1,ncont_pept_near(lj,ib,iprot))
 207         enddo
 208       enddo
 209       call flush(iout)
 210       enddo
 211       return
 212       end
 213 c---------------------------------------------------------------
 214       subroutine calc_ref_beta_signatures(iprot)
 215       implicit none
 216       include 'DIMENSIONS'
 217       include 'DIMENSIONS.ZSCOPT'
 218       include 'DIMENSIONS.COMPAR'
 219       include 'COMMON.IOUNITS'
 220       include 'COMMON.ALLPROT'
 221       include 'COMMON.COMPAR'
 222       include 'COMMON.CLASSES'
 223       include 'COMMON.PEPTCONT'
 224       include 'COMMON.INTERACT'
 225       include 'COMMON.NAMES'
 226       logical lprn /.false./
 227       integer ib,iprot
 228       integer i,ii,ilevel,j,lj,lj1,jc,jc1,k,k1,k2,l
 229       double precision signature
 230       external signature
 231       do ib=1,nclass(iprot)-1
 232       do j=1,nfrag(1,iprot)
 233           if (lprn) write (iout,*) "Protein",iprot," batch",ib,
 234      &      " jfrag=",j," list",(list_frag(k,j,ib,iprot),
 235      &       k=1,nlist_frag(j,ib,iprot))
 236         ii=0
 237         do k=1,nlist_frag(j,ib,iprot)
 238           lj=list_frag(k,j,ib,iprot)
 239           if (isec_ref(lj,iprot).eq.1) then
 240             do l=1,ncont_pept_near(lj,ib,iprot)
 241               jc=icont_pept_near(l,lj,ib,iprot)
 242               call find_dir(lj,jc,lj1,jc1,ib,iprot)
 243               sig_ref(l,lj,ib,iprot)=signature(lj,jc,lj1,jc1,lprn)
 244 c              write (iout,*) "lj",lj," jc",jc," lj1",lj1," jc1",jc1,
 245 c     &         " sig_ref",sig_ref(l,lj,ib,iprot)
 246             enddo
 247           endif
 248         enddo
 249       enddo
 250       write (iout,*)
 251      &   "Reference signatures protein",iprot," class",ib
 252       do j=1,nfrag(1,iprot)
 253         write (iout,*) "fragment",j
 254         do k=1,nlist_frag(j,ib,iprot)
 255           lj=list_frag(k,j,ib,iprot)
 256           write (iout,'(a4,1h(,i2,1h),8f10.5)') restyp(itype(lj)),lj,
 257      &      (sig_ref(ii,lj,ib,iprot),ii=1,ncont_pept_near(lj,ib,iprot))
 258         enddo
 259       enddo
 260       call flush(iout)
 261       enddo
 262
 263       return
 264       end
 265 c---------------------------------------------------------------
 266       double precision function gyrate(ilevel,jfrag,ib,iprot)
 267       implicit real*8 (a-h,o-z)
 268       include 'DIMENSIONS'
 269       include 'DIMENSIONS.ZSCOPT'
 270       include 'DIMENSIONS.COMPAR'
 271       include 'COMMON.INTERACT'
 272       include 'COMMON.CHAIN'
 273       include 'COMMON.IOUNITS'
 274       include 'COMMON.COMPAR'
 275       include 'COMMON.VAR'
 276
 277       double precision cen(3),rg
 278       logical iadded(maxres),lprn /.false./
 279       double precision creff(3,maxres2),cc(3,maxres2)
 280       integer inumber(2,maxres)
 281       common /ccc/ creff,cc,iadded,inumber
 282
 283       IF (ilevel.eq.0) then
 284
 285       do j=1,3
 286        cen(j)=0.0d0
 287       enddo
 288
 289       do i=nnt,nct
 290           do j=1,3
 291             cen(j)=cen(j)+c(j,i)
 292           enddo
 293       enddo
 294       do j=1,3
 295             cen(j)=cen(j)/dble(nct-nnt+1)
 296       enddo
 297       rg = 0.0d0
 298       do i = nnt, nct
 299         do j=1,3
 300          rg = rg + (c(j,i)-cen(j))**2
 301         enddo
 302       end do
 303
 304       gyrate = dsqrt(rg/dble(nct-nnt+1))
 305
 306       ELSE
 307
 308       ii=0
 309       do l=1,nres
 310         iadded(l)=.false.
 311       enddo
 312       do k=1,npiece(jfrag,ilevel,ib,iprot)
 313         if (ilevel.eq.1) then
 314           if (lprn)
 315      &    write (iout,*) "Level 1: jfrag=",jfrag,"k=",k,
 316      &      " adding fragment",
 317      &     ifrag(1,k,jfrag,ib,iprot),ifrag(2,k,jfrag,ib,iprot)
 318           call cprep(ifrag(1,k,jfrag,ib,iprot),
 319      &     ifrag(2,k,jfrag,ib,iprot),0,ii,iprot)
 320         else
 321           kk = ipiece(k,jfrag,ilevel,ib,iprot)
 322           do l=1,npiece(kk,1,ib,iprot)
 323             if (lprn)
 324      &      write (iout,*) "Level",i,": jfrag=",jfrag,"k=",k," kk=",kk,
 325      &       " l=",l," adding fragment",
 326      &       ifrag(1,l,kk,ib,iprot),ifrag(2,l,kk,ib,iprot)
 327             call cprep(ifrag(1,l,kk,ib,iprot),ifrag(2,l,kk,ib,iprot),
 328      &       0,ii,iprot)
 329           enddo
 330         endif
 331       enddo
 332       if (lprn) then
 333         do k=1,ii
 334           write(iout,'(5i4,2(3f10.5,5x))') ilevel,jfrag,k,inumber(1,k),
 335      &         inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
 336         enddo
 337       endif
 338
 339       do j=1,3
 340        cen(j)=0.0d0
 341       enddo
 342
 343       do i=1,ii
 344           do j=1,3
 345             cen(j)=cen(j)+cc(j,i)
 346           enddo
 347       enddo
 348       do j=1,3
 349             cen(j)=cen(j)/dble(ii)
 350       enddo
 351       rg = 0.0d0
 352       do i = 1, ii
 353         do j=1,3
 354          rg = rg + (cc(j,i)-cen(j))**2
 355         enddo
 356       end do
 357
 358       gyrate = dsqrt(rg/dble(ii))
 359
 360       ENDIF
 361
 362       return
 363
 364       end