source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/proc_cont.F

   1       subroutine proc_cont(iprot,nn,*)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.COMPAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.TIME1'
   8       include 'COMMON.SBRIDGE'
   9       include 'COMMON.COMPAR'
  10       include 'COMMON.CHAIN'
  11       include 'COMMON.HEADER'
  12       include 'COMMON.CONTACTS1'
  13       include 'COMMON.PEPTCONT'
  14       include 'COMMON.GEO'
  15       include 'COMMON.CLASSES'
  16       include 'COMMON.PROTNAME'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.INTERACT'
  19       integer i,j,k,ib,iprot,nn,ind,icant,length_frag,ndigit,len_cut,
  20      &  it1,it2,if1,if2
  21       integer ilen
  22       external ilen
  23       character*3 strstr(0:4) /'und','hel','har','std','stp'/
  24 #ifdef DEBUG
  25       write (iout,*) "proc_cont: nlevel",nlevel(iprot)
  26 #endif
  27       if (nlevel(iprot).lt.0) then
  28 #ifdef DEBUG
  29         write (iout,*) "call define_fragments"
  30 #endif
  31         call define_fragments(iprot)
  32       else
  33 #ifdef DEBUG
  34         write (iout,*) "call secondary2"
  35 #endif
  36         call secondary2(.true.,print_secondary,ncont_pept_ref(iprot),
  37      &    icont_pept_ref(1,1,iprot),isec_ref(1,iprot),iprot)
  38       endif
  39       do ib=1,nclass(iprot)-1
  40         if (print_contact) then
  41           write (iout,'(a,i2)') "Structural level",ib+1
  42           write (iout,'(80(1h=))')
  43           write (iout,*) "Electrostatic contacts"
  44         endif
  45         call contacts_between_fragments(print_contact,0,
  46      &   ncont_pept_ref(iprot),
  47      &   icont_pept_ref(1,1,iprot),ncont_frag_ref(1,ib,iprot),
  48      &   icont_frag_ref(1,1,1,ib,iprot),mask_p(1,ib,iprot),ib,iprot)
  49         if (print_contact) then
  50           write (iout,'(80(1h=))')
  51           write (iout,*) "Side chain contacts"
  52         endif
  53         call contacts_between_fragments(print_contact,0,
  54      &   ncont_ref(iprot),
  55      &   icont_ref(1,1,iprot),nsccont_frag_ref(1,ib,iprot),
  56      &   isccont_frag_ref(1,1,1,ib,iprot),mask_sc(1,ib,iprot),ib,iprot)
  57       enddo
  58       if (nlevel(iprot).lt.0) then
  59 #ifdef DEBUG
  60         write(iout,*) "rmscut_base_up",
  61      & (rmscut_base_up(ib),ib=1,nclass(iprot)-1),
  62      & " rmscut_base_low",(rmscut_base_low(ib),ib=1,nclass(iprot)-1),
  63      & " rmsup_lim",(rmsup_lim(ib),ib=1,nclass(iprot)-1)
  64 #endif
  65         do ib=1,nclass(iprot)-1
  66         do i=1,nfrag(1,iprot)
  67           ind=icant(i,i)
  68           len_cut=1000
  69           if (istruct(i,ib,iprot).le.1) then
  70             len_cut=max0(len_frag(i,1,ib,iprot)*4/5,3)
  71           else if (istruct(i,ib,iprot).eq.2.or.istruct(i,ib,iprot).eq.4)
  72      &    then
  73             len_cut=max0(len_frag(i,1,ib,iprot)*2/5,3)
  74           endif
  75 #ifdef DEBUG
  76           write (iout,*) "i",i," istruct",istruct(i,ib,iprot),
  77      &      " ncont_frag",ncont_frag_ref(ind,ib,iprot)," len_cut",
  78      &      len_cut,
  79      &      " icont_single",icont_single(ib)," iloc_single",
  80      &      iloc_single(ib)," isig_match_single",isig_match_single(ib)
  81 #endif
  82           iloc(i,ib,iprot)=iloc_single(ib)
  83           if (istruct(i,ib,iprot).eq.1)
  84      &       isig_match(i,ib,iprot)=isig_match_single(ib)
  85 #ifdef DEBUG
  86           if (iloc(i,ib,iprot).gt.0) write (iout,*)
  87      &     "Local structure used to compare structure of fragment",i,
  88      &     " of protein",iprot," to native."
  89 #endif
  90           if (istruct(i,ib,iprot).ne.3 .and. istruct(i,ib,iprot).ne.0
  91      &        .and. icont_single(ib).gt.0 .and. icontsc_single(ib).eq.0
  92      &        .and. ncont_frag_ref(ind,ib,iprot).ge.len_cut .or.
  93      &        icontp_single(ib).gt.0 .and. ncont_frag_ref(ind,ib,iprot)
  94      &        .gt.0) then
  95 #ifdef DEBUG
  96             write (iout,*) "Electrostatic contacts used to compare",
  97      &     " structure of fragment",i," of protein",iprot," to native."
  98 #endif
  99             ielecont(i,1,ib,iprot)=1
 100             isccont(i,1,ib,iprot)=0
 101           else if (icont_single(ib).gt.0 .and.
 102      &      nsccont_frag_ref(ind,ib,iprot).ge.len_cut .or.
 103      &      icontsc_single(ib).gt.0 .and.
 104      &      nsccont_frag_ref(ind,ib,iprot).gt.0) then
 105 #ifdef DEBUG
 106             write (iout,*) "Side chain contacts used to compare",
 107      &     " structure of fragment",i," of protein",iprot," to native."
 108 #endif
 109             isccont(i,1,ib,iprot)=1
 110             ielecont(i,1,ib,iprot)=0
 111           else
 112 #ifdef DEBUG
 113             write (iout,*) "Contacts not used to compare",
 114      &     " structure of fragment",i," of protein",iprot," to native."
 115 #endif
 116             ielecont(i,1,ib,iprot)=0
 117             isccont(i,1,ib,iprot)=0
 118             nc_req_setf(i,1,ib,iprot)=0
 119           endif
 120           if (iqwol_single(ib).gt.0 .or. isccont(i,1,ib,iprot).eq.0
 121      &         .and. ielecont(i,1,ib,iprot).eq.0) then
 122 #ifdef DEBUG
 123             write (iout,*) "Q used to compare",
 124      &     " structure of fragment",i," of protein",iprot," to native."
 125 #endif
 126             iqwol(i,1,ib,iprot)=1
 127             qcutfrag(i,1,ib,iprot)=qcut_single(ib)
 128             write (iout,*) "qcut_single",qcut_single(ib),
 129      &       "qcutfrag",qcutfrag(i,1,ib,iprot)
 130           else
 131 #ifdef DEBUG
 132             write (iout,*) "Q not used to compare",
 133      &     " structure of fragment",i," of protein",iprot," to native."
 134 #endif
 135             iqwol(i,1,ib,iprot)=0
 136           endif
 137           if (irms_single(ib).gt.0 .or. isccont(i,1,ib,iprot).eq.0
 138      &         .and. ielecont(i,1,ib,iprot).eq.0 .and.
 139      &         iqwol_single(ib).eq.0) then
 140 #ifdef DEBUG
 141             write (iout,*) "RMSD used to compare",
 142      &     " structure of fragment",i," of protein",iprot," to native."
 143 #endif
 144             irms(i,1,ib,iprot)=1
 145           else
 146 #ifdef DEBUG
 147             write (iout,*) "RMSD not used to compare",
 148      &     " structure of fragment",i," of protein",iprot," to native."
 149 #endif
 150             irms(i,1,ib,iprot)=0
 151           endif
 152         enddo
 153         enddo
 154       endif
 155       if (nlevel(iprot).lt.-1) then
 156         call define_pairs(iprot)
 157         nlevel(iprot) = -nlevel(iprot)
 158         if (nlevel(iprot).gt.3) nlevel(iprot)=3
 159         if (nlevel(iprot).eq.3) then
 160           nfrag(3,iprot)=1
 161           do ib=1,nclass(iprot)-1
 162           npiece(1,3,ib,iprot)=nfrag(1,iprot)
 163           do i=1,nfrag(1,iprot)
 164             ipiece(i,1,3,ib,iprot)=i
 165           enddo
 166           ielecont(1,3,ib,iprot)=0
 167           isccont(1,3,ib,iprot)=0
 168           irms(1,3,ib,iprot)=1
 169           n_shift(1,1,3,ib,iprot)=0
 170           n_shift(2,1,3,ib,iprot)=0
 171           enddo
 172         endif
 173       else if (nlevel(iprot).eq.-1) then
 174         nlevel(iprot)=1
 175       endif
 176       isnfrag(1,iprot)=0
 177       do i=1,nlevel(iprot)
 178         isnfrag(i+1,iprot)=isnfrag(i,iprot)+nfrag(i,iprot)
 179       enddo
 180 #ifdef DEBUG
 181       write (iout,*) "nfrag",(nfrag(i,iprot),i=1,nlevel(iprot))
 182       write (iout,*) "isnfrag",(isnfrag(i,iprot),i=1,nlevel(iprot)+1)
 183 #endif
 184       ndigit=3*nfrag(1,iprot)
 185       do i=2,nlevel(iprot)
 186         ndigit=ndigit+2*nfrag(i,iprot)
 187       enddo
 188 #ifdef DEBUG
 189       write (iout,*) "ndigit",ndigit
 190 #endif
 191       if (ndigit.ne.nn) then
 192         write (iout,*) "Error - length of the class mask of protein",
 193      &  iprot," is ",nn," but",ndigit,
 194      &  " fields found in the class template."
 195         return1
 196       endif
 197       if (.not.binary(iprot) .and. ndigit.gt.30) then
 198         write (iout,*) "Highest class too large for protein",iprot,
 199      &    " ; switching to binary representation."
 200         binary(iprot)=.true.
 201       endif
 202
 203       write (iout,'(/80(1h=)/a,i3,1h(,1x,a,2h)./80(1h-))')
 204      &  "Specification of fragments and cut-off criteria of protein",
 205      &  iprot,protname(iprot)(:ilen(protname(iprot)))
 206       write (iout,"(a)") "Secondary-structure codes are as follows:"
 207       write (iout,'(a)')
 208      &  "und - undefined; hel - helix; har - beta-hairpin;"
 209       write (iout,'(a)') "std - single strand; stp - pair of strands."
 210       do ib=1,nclass(iprot)-1
 211       write (iout,'(80(1h-)/a,i3)') "Structural level",ib
 212       do i=1,nlevel(iprot)
 213         do j=1,nfrag(i,iprot)
 214           length_frag = 0
 215           if (i.eq.1) then
 216             do k=1,npiece(j,i,ib,iprot)
 217               length_frag=length_frag+ifrag(2,k,j,ib,iprot)
 218      &          -ifrag(1,k,j,ib,iprot)+1
 219             enddo
 220           else
 221             do k=1,npiece(j,i,ib,iprot)
 222               length_frag=length_frag
 223      &         +len_frag(ipiece(k,j,i,ib,iprot),1,ib,iprot)
 224             enddo
 225           endif
 226           len_frag(j,i,ib,iprot)=length_frag
 227           rmscutfrag(1,j,i,ib,iprot)=rmscut_base_up(ib)*length_frag
 228           rmscutfrag(2,j,i,ib,iprot)=rmscut_base_low(ib)*length_frag
 229           if (rmscutfrag(1,j,i,ib,iprot).lt.rmsup_lim(ib))
 230      &      rmscutfrag(1,j,i,ib,iprot)=rmsup_lim(ib)
 231           if (rmscutfrag(1,j,i,ib,iprot).gt.rmsupup_lim(ib))
 232      &      rmscutfrag(1,j,i,ib,iprot)=rmsupup_lim(ib)
 233         enddo
 234       enddo
 235       write (iout,'(80(1h-)/a,i3,a,i3/)')
 236      &   "Level",1," number of fragments:",nfrag(1,iprot)
 237       write (iout,'(2a5,6a4,a5,a7,a7,a12,2a5,a5,a4,a4,2x,a)')
 238      &  'frag','len','loc','sig','eC','scC','Q','rms','qcut','rmscut',
 239      &  'shifts','angcut','frS','frC','nC','str',
 240      &  'np','specification'
 241       do j=1,nfrag(1,iprot)
 242         write (iout,'(2i5,6i4,1x,3f4.1,2i3,2f6.1,2f5.2,i5,a4,i4,$)') j,
 243      &    len_frag(j,1,ib,iprot),iloc(j,ib,iprot),
 244      &    isig_match(j,ib,iprot),ielecont(j,1,ib,iprot),
 245      &    isccont(j,1,ib,iprot),iqwol(j,1,ib,iprot),irms(j,1,ib,iprot),
 246      &    qcutfrag(j,1,ib,iprot),rmscutfrag(1,j,1,ib,iprot),
 247      &    rmscutfrag(2,j,1,ib,iprot),n_shift(1,j,1,ib,iprot),
 248      &    n_shift(2,j,1,ib,iprot),ang_cut(j,ib,iprot)*rad2deg,
 249      &    ang_cut1(j,ib,iprot)*rad2deg,frac_min(j,ib,iprot),
 250      &    nc_fragm(j,1,ib,iprot),nc_req_setf(j,1,ib,iprot),
 251      &    strstr(istruct(j,ib,iprot)),npiece(j,1,ib,iprot)
 252         do k=1,npiece(j,1,ib,iprot)
 253           if1=ifrag(1,k,j,ib,iprot)
 254           it1=itype(if1)
 255           if2=ifrag(2,k,j,ib,iprot)
 256           it2=itype(if2)
 257           write (iout,'(2x,a1,i3,1h-,a1,i3,$)') onelet(it1),if1,
 258      &      onelet(it2),if2
 259         enddo
 260         write (iout,*)
 261       enddo
 262       do i=2,nlevel(iprot)
 263         write (iout,'(80(1h-)/a,i3,a,i3/)')
 264      &   "Level",i," number of fragments:",nfrag(i,iprot)
 265         write (iout,'(2a5,4a4,a5,a7,a7,2a5,a4,2x,a)')
 266      &  'frag','len','eC','scC','Q','rms','Qcut','rmscut','shifts',
 267      &  'frC','nC','np','specification'
 268         do j=1,nfrag(i,iprot)
 269           write (iout,'(2i5,4i4,1x,3f4.1,2i3,f5.2,i5,i4,$)')
 270      &      j,len_frag(j,i,ib,iprot),ielecont(j,i,ib,iprot),
 271      &      isccont(j,i,ib,iprot),iqwol(j,i,ib,iprot),
 272      &      irms(j,i,ib,iprot),qcutfrag(j,i,ib,iprot),
 273      &      rmscutfrag(1,j,i,ib,iprot),rmscutfrag(2,j,i,ib,iprot),
 274      &      n_shift(1,j,i,ib,iprot),n_shift(2,j,i,ib,iprot),
 275      &      nc_fragm(j,i,ib,iprot),nc_req_setf(j,i,ib,iprot),
 276      &      npiece(j,i,ib,iprot)
 277           do k=1,npiece(j,i,ib,iprot)
 278             write (iout,'(i3,$)') ipiece(k,j,i,ib,iprot)
 279           enddo
 280           write (iout,*)
 281         enddo
 282       enddo
 283       enddo
 284       write (iout,'(80(1h=))')
 285       return
 286       end