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