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

   1       subroutine conf_compar(jcon,ib,iprot,nn,iscor,sbin,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.PEPTCONT'
  12       include 'COMMON.CONTACTS1'
  13       include 'COMMON.HEADER'
  14       include 'COMMON.CLASSES'
  15       include 'COMMON.GEO'
  16       include 'COMMON.ENERGIES'
  17       logical lprn,lprn1,lprn2,lsig_match
  18       integer ncont_frag(mmaxfrag),
  19      & icont_frag(2,maxcont,mmaxfrag),ncontsc,
  20      & icontsc(1,maxcont),nsccont_frag(mmaxfrag),
  21      & isccont_frag(2,maxcont,mmaxfrag)
  22       integer jcon,ib,iprot,iscor,nn,lat
  23       integer i,j,k,ik,kk,ind,icant,ncnat,nsec_match,nsec_nomatch,
  24      &  ishif1,ishif2,
  25      &  nc_match,ncon_match,ishif,iclass_con,ishifft_con,iclass_rms,
  26      &  iclass_q,ishifft_rms,ishiff,idig,iex,im
  27       integer isecstr(maxres)
  28       integer itemp(maxfrag)
  29       double precision rmscalc,rms,rmsnat,qwolynes
  30       double precision sig_frag(maxfrag)
  31       character*4 liczba
  32       character*64 sbin
  33       lprn1=lprn .and. print_contact
  34       lprn2=lprn .and. print_secondary
  35 c      print *,"Enter conf_compar",jcon
  36       call angnorm12(rmsang,iprot)
  37       rmsang=rmsang*rad2deg
  38 c      if (lprn) then
  39 c        write (liczba,'(bz,i4.4)') jcon
  40 c        open(ipdb,file=prefintin(:ilen(prefintin))//liczba//".pdb")
  41 c        call pdbout(energy(jcon),titel)
  42 c        write(iout,*) "Protein",iprot," conformation",jcon,
  43 c     &    " nlevel",nlevel(iprot)
  44 c      endif
  45 c      if (lprn) then
  46 c        write (iout,*) "CONF_COMPAR: Complete reference structure"
  47 c        do i=1,nres
  48 c          write(iout,'(i4,3f10.5)') i,(cref(j,i,iprot),j=1,3)
  49 c        enddo
  50 c      endif
  51 c Level 1: check secondary and supersecondary structure
  52       call elecont(lprn2,ncont,icont,nnt,nct)
  53       call secondary2(lprn1,.false.,ncont,icont,isecstr,iprot)
  54       call contact(lprn2,ncontsc,icontsc,nnt,nct)
  55       if (lprn) write(iout,*) "Assigning electrostatic contacts"
  56       call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,
  57      &   icont_frag,mask_p(1,ib,iprot),ib,iprot)
  58       if (lprn) write(iout,*) "Assigning sidechain contacts"
  59       call contacts_between_fragments(lprn,3,ncontsc,icontsc,
  60      &   nsccont_frag,isccont_frag,mask_sc(1,ib,iprot),ib,iprot)
  61       do i=1,nlevel(iprot)
  62         do j=1,isnfrag(nlevel(iprot)+1,iprot)
  63           iclass1(j,i)=0
  64         enddo
  65       enddo
  66       do j=1,nfrag(1,iprot)
  67         ind = icant(j,j)
  68         if (lprn) then
  69           write (iout,'(80(1h=))')
  70           write (iout,*) "Level",1," fragment",j
  71           write (iout,'(80(1h=))')
  72         endif
  73         rmsfrag(j,1)=rmscalc(0,1,j,jcon,ib,iprot,lprn)
  74 c Compare electrostatic contacts in the current conf with that in the native
  75 c structure.
  76         if (lprn) write (iout,*)
  77      &    "Comparing electrostatic contact map and local structure",
  78      &    " conformation",jcon," protein",iprot
  79         ncnat=ncont_frag_ref(ind,ib,iprot)
  80 c        write (iout,*) "before match_contact:",nc_fragm(j,1,ib,iprot),
  81 c     &   nc_req_setf(j,1,ib,iprot)
  82         call match_secondary(j,isecstr,nsec_match,nsec_nomatch,
  83      &    sig_frag(j),ib,iprot,lprn)
  84         if (lprn) write (iout,*) "Fragment",j," nsec_match",
  85      &    nsec_match," nsec_nomatch",nsec_nomatch,
  86      &    " length",len_frag(j,1,ib,iprot)," min_len",
  87      &    frac_sec(ib,iprot)*len_frag(j,1,ib,iprot)," sig_frag",
  88      &     sig_frag(j)
  89         if (nsec_match.ge.frac_sec(ib,iprot)*len_frag(j,1,ib,iprot))
  90      &  then
  91           iclass1(j,1)=1
  92           if (lprn) write (iout,*) "Fragment",j,
  93      &      " has correct secondary structure"
  94         else if (nsec_nomatch.ge.
  95      &      frac_sec(ib,iprot)*len_frag(j,1,ib,iprot)/2) then
  96           iclass1(j,1)=4
  97           if (lprn) write (iout,*) "Fragment",j,
  98      &      " has wrong secondary structure"
  99         else
 100           iclass1(j,1)=0
 101           if (lprn) write (iout,*) "Fragment",j,
 102      &      " has grossly incorrect secondary structure"
 103         endif
 104         if (ielecont(j,1,ib,iprot).gt.0) then
 105           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 106      &     ncont_frag_ref(ind,ib,iprot),
 107      &     icont_frag_ref(1,1,ind,ib,iprot),
 108      &     ncont_frag(ind),icont_frag(1,1,ind),
 109      &     j,n_shift(1,j,1,ib,iprot),n_shift(2,j,1,ib,iprot),
 110      &     nc_fragm(j,1,ib,iprot),nc_req_setf(j,1,ib,iprot),
 111      &     istruct(j,ib,iprot),.true.,ib,iprot,lprn)
 112         else if (isccont(j,1,ib,iprot).gt.0) then
 113           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 114      &      nsccont_frag_ref(ind,ib,iprot),
 115      &      isccont_frag_ref(1,1,ind,ib,iprot),
 116      &      nsccont_frag(ind),isccont_frag(1,1,ind),j,
 117      &      n_shift(1,j,1,ib,iprot),n_shift(2,j,1,ib,iprot),
 118      &      nc_fragm(j,1,ib,iprot),
 119      &      nc_req_setf(j,1,ib,iprot),istruct(j,ib,iprot),.true.,ib,
 120      &      iprot,lprn)
 121         else if (iloc(j,ib,iprot).gt.0) then
 122           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 123      &      0,icont_frag_ref(1,1,ind,ib,iprot),ncont_frag(ind),
 124      &      icont_frag(1,1,ind),j,n_shift(1,j,1,ib,iprot),
 125      &      n_shift(2,j,1,ib,iprot),nc_fragm(j,1,ib,iprot),
 126      &      0,istruct(j,ib,iprot),.true.,ib,iprot,lprn)
 127         else
 128           ishif=0
 129           nc_match=1
 130         endif
 131         qfrag(j,1)=qwolynes(1,j,ib,iprot)
 132         if (iqwol(j,1,ib,iprot).gt.0) then
 133           if (qfrag(j,1).le.qcutfrag(j,1,ib,iprot)) then
 134             iclass_q = 2
 135           else
 136             iclass_q = 0
 137           endif
 138         else
 139           iclass_q = 2
 140         endif
 141         lsig_match=isig_match(j,ib,iprot).eq.0 .or. sig_frag(j).gt.0.0d0
 142         if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
 143         ishif=ishif1
 144         if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 145         if (lprn) write (iout,*) "Ishift",ishif," nc_match",nc_match
 146 c        write (iout,*) iprot,j," irms",irms(j,1,ib,iprot)
 147         if (irms(j,1,ib,iprot).gt.0) then
 148           if (rmsfrag(j,1).le.rmscutfrag(1,j,1,ib,iprot)) then
 149             iclass_rms=2
 150             ishifft_rms=0
 151           else
 152             ishiff=0
 153             rms=1.0d2
 154             iclass_rms=0
 155             do while (rms.gt.rmscutfrag(1,j,1,ib,iprot) .and.
 156      &         ishiff.lt.n_shift(1,j,1,ib,iprot))
 157               ishiff=ishiff+1
 158               rms=rmscalc(-ishiff,1,j,jcon,ib,iprot,lprn)
 159 c              write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
 160 c     &          " rms",rms," rmscut",rmscutfrag(1,j,1,ib,iprot)
 161               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 162               if (rms.gt.rmscutfrag(1,j,1,ib,iprot)) then
 163                 rms=rmscalc(ishiff,1,j,jcon,ib,iprot,lprn)
 164 c                write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
 165 c     &           " rms",rms
 166               endif
 167               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 168             enddo
 169 c            write (iout,*) "After loop: rms",rms,
 170 c     &        " rmscut",rmscutfrag(1,j,1,ib,iprot)
 171 c            write (iout,*) "iclass_rms",iclass_rms
 172             if (rms.le.rmscutfrag(1,j,1,ib,iprot)) then
 173               ishifft_rms=ishiff
 174               rmsfrag(j,1)=rms
 175               iclass_rms=1
 176             endif
 177 c            write (iout,*) "iclass_rms",iclass_rms
 178           endif
 179 c          write (iout,*) "ishif",ishif
 180           if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
 181         else
 182           iclass_rms=1
 183         endif
 184 c        write (iout,*) "ishif",ishif," iclass",iclass1(j,1),
 185 c     &    " iclass_rms",iclass_rms
 186         if (iclass1(j,1).ne.4 .and. nc_match.gt.0 .and.
 187      &      iclass_rms.gt.0 .and. iclass_q.gt.0 .and. lsig_match) then
 188           if (ishif.eq.0) then
 189             iclass1(j,1)=iclass1(j,1)+6
 190           else
 191             iclass1(j,1)=iclass1(j,1)+2
 192           endif
 193         endif
 194         ncont_nat(1,j,1)=nc_match
 195         ncont_nat(2,j,1)=ncon_match
 196         ishifft(j,1)=ishif
 197       enddo
 198 c Next levels: Check arrangements of elementary fragments.
 199       do i=2,nlevel(iprot)
 200         do j=1,nfrag(i,iprot)
 201         if (i .eq. 2) ind = icant(ipiece(1,j,i,ib,iprot),
 202      &    ipiece(2,j,i,ib,iprot))
 203         if (lprn) then
 204             write (iout,'(80(1h=))')
 205             write (iout,*) "Level",i," fragment",j
 206             write (iout,'(80(1h=))')
 207         endif
 208 c If an elementary fragment doesn't exist, don't check higher hierarchy levels.
 209 c 3/4/03 AL No, we consider the existence of a composite fragment even if
 210 c the corresponding elementary fragments are incomplete.
 211 c
 212 c        do k=1,npiece(j,i,ib,iprot)
 213 c          ik=ipiece(k,j,i,ib,iprot)
 214 c          if (iclass1(ik,1).eq.0) then
 215 c            iclass1(j,i)=0
 216 c            goto 12
 217 c          endif
 218 c        enddo
 219         if (i.eq.2 .and. ielecont(j,i,ib,iprot).gt.0) then
 220           iclass_con=0
 221           ishifft_con=0
 222           if (lprn) write (iout,*)
 223      &     "Comparing electrostatic contact map: fragments",
 224      &      ipiece(1,j,i,ib,iprot),ipiece(2,j,i,ib,iprot)," ind",ind
 225           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 226      &     ncont_frag_ref(ind,ib,iprot),
 227      &     icont_frag_ref(1,1,ind,ib,iprot),
 228      &     ncont_frag(ind),icont_frag(1,1,ind),
 229      &     j,n_shift(1,j,i,ib,iprot),n_shift(2,j,i,ib,iprot),
 230      &     nc_fragm(j,i,ib,iprot),
 231      &     nc_req_setf(j,i,ib,iprot),2,.false.,ib,iprot,lprn)
 232           ishif=ishif1
 233           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 234           if (nc_match.gt.0) then
 235             if (ishif.eq.0) then
 236               iclass_con=2
 237             else
 238               iclass_con=1
 239             endif
 240           endif
 241           ncont_nat(1,j,i)=nc_match
 242           ncont_nat(2,j,i)=ncon_match
 243           ishifft_con=ishif
 244         else if (i.eq.2 .and. isccont(j,i,ib,iprot).gt.0) then
 245           iclass_con=0
 246           ishifft_con=0
 247           if (lprn) write (iout,*)
 248      &     "Comparing sidechain contact map: fragments",
 249      &    ipiece(1,j,i,ib,iprot),ipiece(2,j,i,ib,iprot)," ind",ind
 250           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 251      &     nsccont_frag_ref(ind,ib,iprot),
 252      &     isccont_frag_ref(1,1,ind,ib,iprot),
 253      &     nsccont_frag(ind),isccont_frag(1,1,ind),
 254      &     j,n_shift(1,j,i,ib,iprot),n_shift(2,j,i,ib,iprot),
 255      &     nc_fragm(j,i,ib,iprot),
 256      &     nc_req_setf(j,i,ib,iprot),2,.false.,ib,iprot,lprn)
 257           ishif=ishif1
 258           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 259           if (nc_match.gt.0) then
 260             if (ishif.eq.0) then
 261               iclass_con=2
 262             else
 263               iclass_con=1
 264             endif
 265           endif
 266           ncont_nat(1,j,i)=nc_match
 267           ncont_nat(2,j,i)=ncon_match
 268           ishifft_con=ishif
 269         else if (i.eq.2) then
 270           iclass_con=2
 271           ishifft_con=0
 272         endif
 273         if (i.eq.2) then
 274           qfrag(j,2)=qwolynes(2,j,ib,iprot)
 275           if (iqwol(j,i,ib,iprot).gt.0) then
 276             if (qfrag(j,2).le.qcutfrag(j,2,ib,iprot)) then
 277               iclass_q = 2
 278             else
 279               iclass_q = 0
 280             endif
 281           else
 282             iclass_q = 2
 283           endif
 284         else
 285           iclass_q = 2
 286         endif
 287         if (irms(j,i,ib,iprot).gt.0) then
 288           iclass_rms=0
 289           ishifft_rms=0
 290           if (lprn) write (iout,*)
 291      &     "Comparing rms: fragments",
 292      &     (ipiece(k,j,i,ib,iprot),k=1,npiece(j,i,ib,iprot))
 293           rmsfrag(j,i)=rmscalc(0,i,j,jcon,ib,iprot,lprn)
 294           if (lprn) write (iout,*) "rms",rmsfrag(j,i),
 295      &      " rmscut",rmscutfrag(1,j,i,ib,iprot)
 296           if (rmsfrag(j,i).le.rmscutfrag(1,j,i,ib,iprot)) then
 297             iclass_rms=2
 298             ishifft_rms=0
 299           else
 300             ishif=0
 301             rms=1.0d2
 302             do while (rms.gt.rmscutfrag(1,j,i,ib,iprot) .and.
 303      &         ishif.lt.n_shift(1,j,i,ib,iprot))
 304               ishif=ishif+1
 305               rms=rmscalc(-ishif,i,j,jcon,ib,iprot,lprn)
 306               if (lprn) write (iout,*) "rms",rms,
 307      &          " rmscut",rmscutfrag(1,j,i,ib,iprot)
 308 c              print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
 309               if (rms.gt.rmscutfrag(1,j,i,ib,iprot)) then
 310                 rms=rmscalc(ishif,i,j,jcon,ib,iprot,lprn)
 311 c                print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
 312               endif
 313               if (lprn) write (iout,*) "rms",rms,
 314      &          " rmscut",rmscutfrag(1,j,i,ib,iprot)
 315             enddo
 316             if (rms.le.rmscutfrag(1,j,i,ib,iprot)) then
 317               ishifft_rms=ishif
 318               rmsfrag(j,i)=rms
 319               iclass_rms=1
 320             endif
 321           endif
 322         endif
 323         if (irms(j,i,ib,iprot).eq.0 .and. ielecont(j,i,ib,iprot).eq.0
 324      &    .and. isccont(j,i,ib,iprot).eq.0
 325      &    .and. iqwol(j,i,ib,iprot).eq.0) then
 326           write (iout,*) "Error: no measure of comparison specified:",
 327      &      " level",i," part",j
 328           stop
 329         endif
 330         if (lprn)
 331      &  write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
 332         if (i.eq.2) then
 333           iclass1(j,i) = min0(iclass_con,iclass_rms,iclass_q)
 334           ishifft(j,i)= max0(ishifft_con,ishifft_rms)
 335         else if (i.gt.2) then
 336           iclass1(j,i) = iclass_rms
 337           ishifft(j,i)= ishifft_rms
 338         endif
 339    12   continue
 340         enddo
 341       enddo
 342       rms_nat=rmsnat(jcon,iprot)
 343       qnat=qwolynes(0,0,ib,iprot)
 344       if (lprn) write (iout,*) "rmsnat",rms_nat," qnat",qnat
 345 C Compute the structural class
 346       iscor=0
 347       do i=1,nlevel(iprot)
 348         IF (I.EQ.1) THEN
 349         do j=1,nfrag(i,iprot)
 350           itemp(j)=iclass1(j,i)
 351         enddo
 352         do kk=-1,1
 353           do j=1,nfrag(i,iprot)
 354             idig = 2*isnfrag(nlevel(iprot)+1,iprot)-2*isnfrag(i,iprot)
 355      &        -kk*nfrag(i,iprot)-j
 356             iex = 2**idig
 357             im=mod(itemp(j),2)
 358             itemp(j)=itemp(j)/2
 359 c            write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 360 c     &        " iclass",iclass1(j,i)," im",im
 361             if (.not.binary(iprot)) iscor=iscor+im*iex
 362             write (sbin(nn-idig:nn-idig),'(i1)') im
 363           enddo
 364         enddo
 365         ELSE
 366         do j=1,nfrag(i,iprot)
 367           idig = 2*isnfrag(nlevel(iprot)+1,iprot)-2*isnfrag(i,iprot)-j
 368           iex = 2**idig
 369           if (iclass1(j,i).gt.0) then
 370             im=1
 371           else
 372             im=0
 373           endif
 374 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 375 c     &      " iclass",iclass1(j,i)," im",im
 376           if (.not. binary(iprot)) iscor=iscor+im*iex
 377           write (sbin(nn-idig:nn-idig),'(i1)') im
 378         enddo
 379         do j=1,nfrag(i,iprot)
 380           idig = 2*isnfrag(nlevel(iprot)+1,iprot)-2*isnfrag(i,iprot)
 381      &      -nfrag(i,iprot)-j
 382           iex = 2**idig
 383           if (iclass1(j,i).gt.1) then
 384             im=1
 385           else
 386             im=0
 387           endif
 388 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 389 c     &      " iclass",iclass1(j,i)," im",im
 390           write (sbin(nn-idig:nn-idig),'(i1)') im
 391           if (.not. binary(iprot)) iscor=iscor+im*iex
 392         enddo
 393         ENDIF
 394       enddo
 395       if (lprn) then
 396       write (iout,'(i5,$)') jcon
 397       do i=1,nlevel(iprot)
 398         write (iout,'(i5,$)') i
 399         if (i.eq.1) then
 400         do j=1,nfrag(i,iprot)
 401           write (iout,'(2i4,f6.2,i3,$)') ncont_nat(1,j,i),
 402      &     ncont_nat(2,j,i),rmsfrag(j,i),ishifft(j,i)
 403         enddo
 404         else
 405         do j=1,nfrag(i,iprot)
 406           write (iout,'(f6.2,i3,$)') rmsfrag(j,i),ishifft(j,i)
 407         enddo
 408         endif
 409         write (iout,'(" ",$)')
 410         do j=1,nfrag(i,iprot)
 411           write (iout,'(i1,$)') iclass1(j,i)
 412         enddo
 413       enddo
 414       if (binary(iprot)) then
 415         write (iout,'("  ",$)')
 416         do j=1,nlevel(iprot)
 417           write (iout,'(100(i1,$))')(iclass1(k,j),k=1,
 418      &      nfrag(j,iprot),iprot)
 419           if (j.lt.nlevel(iprot)) write(iout,'(".",$)')
 420         enddo
 421         write (iout,'(f6.2)') rms_nat
 422       else
 423         write (iout,'(i10,f6.2)') iscor,rms_nat
 424       endif
 425       endif
 426
 427       if (print_class) then
 428           write(istat,'(i5,2f10.2,f8.3,2f6.3,$)')
 429      &      jcon,eini(jcon,iprot),entfac(jcon,iprot),
 430      &      rms_nat,qnat,rmsang/(nct-nnt-2)
 431           do j=1,nlevel(iprot)
 432             write(istat,'(1x,20(i3,$))')
 433      &        (ncont_nat(1,k,j),k=1,nfrag(j,iprot))
 434             if (j.eq.1) write (istat,'(1x,f4.1,$)')
 435      &        (sig_frag(k),k=1,nfrag(j,iprot))
 436             if (j.lt.3) then
 437               write (istat,'(1x,20(f5.1,f5.2,$))')
 438      &          (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j,iprot))
 439             else
 440               write(istat,'(1x,20(f5.1,$))')
 441      &          (rmsfrag(k,j),k=1,nfrag(j,iprot))
 442             endif
 443             write(istat,'(1x,20(i1,$))')
 444      &        (iclass1(k,j),k=1,nfrag(j,iprot))
 445           enddo
 446           if (binary(iprot)) then
 447             write (istat,'("  ",$)')
 448             do j=1,nlevel(iprot)
 449               write (istat,'(100(i1,$))')(iclass1(k,j),
 450      &           k=1,nfrag(j,iprot))
 451               if (j.lt.nlevel(iprot)) write(istat,'(".",$)')
 452             enddo
 453             write (istat,'(i10,i3)') iscore(jcon,0,iprot),ib
 454           else
 455             write (istat,'(2i10,i3)') iscor,iscore(jcon,0,iprot),ib
 456           endif
 457       endif
 458       RETURN
 459       END