source/maxlik/src_MD_T_maxlik-NEWCORR.safe/secondary.F

   1       subroutine define_fragments(iprot)
   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.FRAG'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.CONTROL'
  11       include 'COMMON.COMPAR'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.HEADER'
  14       include 'COMMON.GEO'
  15       include 'COMMON.CONTACTS1'
  16       include 'COMMON.PEPTCONT'
  17       include 'COMMON.INTERACT'
  18       include 'COMMON.NAMES'
  19       include 'COMMON.CLASSES'
  20       integer ibatch,iprot
  21       integer i,j,ii,i1,i2,i3,i4,it1,it2,it3,it4
  22       integer nstrand,istrand(2,maxres/2)
  23       integer nhairp,ihairp(2,maxres/5)
  24       character*16 strstr(4) /'helix','hairpin','strand','strand pair'/
  25 #ifdef DEBUG
  26       write (iout,*) "Entered DEFINE_FRAGMENTS iprot",iprot
  27 #endif
  28       do ibatch=1,nclass(iprot)-1
  29 #ifdef DEBUG
  30       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel(ibatch),
  31      &               ' NC_REQ_HEL',ncreq_hel(ibatch),
  32      &               ' NC_FRAC_BET',ncfrac_bet(ibatch),
  33      &               ' NC_REQ_BET',ncreq_bet(ibatch),
  34      &               ' NC_FRAC_PAIR',ncfrac_pair(ibatch),
  35      &               ' NC_REQ_PAIR',ncreq_pair(ibatch),
  36      &  ' RMS_PAIR',irms_pair(ibatch),' SPLIT_BET',isplit_bet
  37       write (iout,*) 'NSHIFT_HEL',nshift_hel(ibatch),
  38      &               ' NSHIFT_BET',nshift_bet(ibatch),
  39      &               ' NSHIFT_STRAND',nshift_strand(ibatch),
  40      &               ' NSHIFT_PAIR',nshift_pair(ibatch)
  41       write (iout,*) 'ANGCUT_HEL',angcut_hel(ibatch)*rad2deg,
  42      &  ' MAXANG_HEL',angcut1_hel(ibatch)*rad2deg
  43       write (iout,*) 'ANGCUT_BET',angcut_bet(ibatch)*rad2deg,
  44      &               ' MAXANG_BET',angcut1_bet(ibatch)*rad2deg
  45       write (iout,*) 'ANGCUT_STRAND',angcut_strand(ibatch)*rad2deg,
  46      &               ' MAXANG_STRAND',angcut1_strand(ibatch)*rad2deg
  47       write (iout,*) 'FRAC_MIN',frac_min_set(ibatch)
  48 #endif
  49 c Find secondary structure elements (helices and beta-sheets)
  50 #ifdef DEBUG
  51       write (iout,*) "Calling SECONDARY2 iprot",iprot
  52 #endif
  53       enddo
  54       call secondary2(.true.,print_secondary,ncont_pept_ref(iprot),
  55      &  icont_pept_ref(1,1,iprot),isec_ref(1,iprot),iprot)
  56 c Define primary fragments. First include the helices.
  57       nhairp=0
  58       nstrand=0
  59 c Merge helices
  60 c AL 12/23/03 - to avoid splitting helices into very small fragments
  61       if (merge_helices) then
  62       if (print_secondary) then
  63       write (iout,*) "Before merging helices: nhfrag",nhfrag
  64       do i=1,nhfrag
  65         write (2,*) hfrag(1,i),hfrag(2,i)
  66       enddo
  67       endif
  68       i=1
  69       do while (i.lt.nhfrag)
  70         if (hfrag(1,i+1)-hfrag(2,i).le.1) then
  71           nhfrag=nhfrag-1
  72           hfrag(2,i)=hfrag(2,i+1)
  73           do j=i+1,nhfrag
  74             hfrag(1,j)=hfrag(1,j+1)
  75             hfrag(2,j)=hfrag(2,j+1)
  76           enddo
  77         endif
  78         i=i+1
  79       enddo
  80       if (print_secondary) then
  81       write (iout,*) "After merging helices: nhfrag",nhfrag
  82       do i=1,nhfrag
  83         write (2,*) hfrag(1,i),hfrag(2,i)
  84       enddo
  85       endif
  86       endif
  87       nfrag(1,iprot)=nhfrag
  88       do ibatch=1,nclass(iprot)-1
  89         do i=1,nhfrag
  90           npiece(i,1,ibatch,iprot)=1
  91           ifrag(1,1,i,ibatch,iprot)=hfrag(1,i)
  92           ifrag(2,1,i,ibatch,iprot)=hfrag(2,i)
  93           n_shift(1,i,1,ibatch,iprot)=0
  94           n_shift(2,i,1,ibatch,iprot)=nshift_hel(ibatch)
  95           ang_cut(i,ibatch,iprot)=angcut_hel(ibatch)
  96           ang_cut1(i,ibatch,iprot)=angcut1_hel(ibatch)
  97           frac_min(i,ibatch,iprot)=frac_min_set(ibatch)
  98           nc_fragm(i,1,ibatch,iprot)=ncfrac_hel(ibatch)
  99           nc_req_setf(i,1,ibatch,iprot)=ncreq_hel(ibatch)
 100           istruct(i,ibatch,iprot)=1
 101         enddo
 102       enddo
 103 #ifdef DEBUG
 104       write (iout,*) "isplit_bet",isplit_bet
 105 #endif
 106       if (isplit_bet.gt.1) then
 107 c Split beta-sheets into strands and store strands as primary fragments.
 108         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
 109         do ibatch=1,nclass(iprot)-1
 110           do i=1,nstrand
 111             ii=i+nfrag(1,iprot)
 112             npiece(ii,1,ibatch,iprot)=1
 113             ifrag(1,1,ii,ibatch,iprot)=istrand(1,i)
 114             ifrag(2,1,ii,ibatch,iprot)=istrand(2,i)
 115             n_shift(1,ii,1,ibatch,iprot)=nshift_strand(ibatch)
 116             n_shift(2,ii,1,ibatch,iprot)=nshift_strand(ibatch)
 117             ang_cut(ii,ibatch,iprot)=angcut_strand(ibatch)
 118             ang_cut1(ii,ibatch,iprot)=angcut1_strand(ibatch)
 119             frac_min(ii,ibatch,iprot)=frac_min_set(ibatch)
 120             nc_fragm(ii,1,ibatch,iprot)=0
 121             nc_req_setf(ii,1,ibatch,iprot)=0
 122             istruct(ii,ibatch,iprot)=3
 123           enddo
 124         enddo
 125         nfrag(1,iprot)=nfrag(1,iprot)+nstrand
 126       else if (isplit_bet.eq.1) then
 127 c Split only far beta-sheets; does not split hairpins.
 128         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
 129         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
 130         do ibatch=1,nclass(iprot)-1
 131           do i=1,nhairp
 132             ii=i+nfrag(1,iprot)
 133             npiece(ii,1,ibatch,iprot)=1
 134             ifrag(1,1,ii,ibatch,iprot)=ihairp(1,i)
 135             ifrag(2,1,ii,ibatch,iprot)=ihairp(2,i)
 136             n_shift(1,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 137             n_shift(2,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 138             ang_cut(ii,ibatch,iprot)=angcut_bet(ibatch)
 139             ang_cut1(ii,ibatch,iprot)=angcut1_bet(ibatch)
 140             frac_min(ii,ibatch,iprot)=frac_min_set(ibatch)
 141             nc_fragm(ii,1,ibatch,iprot)=ncfrac_bet(ibatch)
 142             nc_req_setf(ii,1,ibatch,iprot)=ncreq_bet(ibatch)
 143             istruct(ii,ibatch,iprot)=2
 144           enddo
 145         enddo
 146         nfrag(1,iprot)=nfrag(1,iprot)+nhairp
 147         do ibatch=1,nclass(iprot)-1
 148           do i=1,nstrand
 149             ii=i+nfrag(1,iprot)
 150             npiece(ii,1,ibatch,iprot)=1
 151             ifrag(1,1,ii,ibatch,iprot)=istrand(1,i)
 152             ifrag(2,1,ii,ibatch,iprot)=istrand(2,i)
 153             n_shift(1,ii,1,ibatch,iprot)=nshift_strand(ibatch)
 154             n_shift(2,ii,1,ibatch,iprot)=nshift_strand(ibatch)
 155             ang_cut(ii,ibatch,iprot)=angcut_strand(ibatch)
 156             ang_cut1(ii,ibatch,iprot)=angcut1_strand(ibatch)
 157             frac_min(ii,ibatch,iprot)=frac_min_set(ibatch)
 158             nc_fragm(ii,1,ibatch,iprot)=0
 159             nc_req_setf(ii,1,ibatch,iprot)=0
 160             istruct(ii,ibatch,iprot)=3
 161           enddo
 162         enddo
 163         nfrag(1,iprot)=nfrag(1,iprot)+nstrand
 164       else
 165 c Do not split beta-sheets; each pair of strands is a primary element.
 166         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
 167         do ibatch=1,nclass(iprot)-1
 168           do i=1,nhairp
 169             ii=i+nfrag(1,iprot)
 170             npiece(ii,1,ibatch,iprot)=1
 171             ifrag(1,1,ii,ibatch,iprot)=ihairp(1,i)
 172             ifrag(2,1,ii,ibatch,iprot)=ihairp(2,i)
 173             n_shift(1,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 174             n_shift(2,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 175             ang_cut(ii,ibatch,iprot)=angcut_bet(ibatch)
 176             ang_cut1(ii,ibatch,iprot)=angcut1_bet(ibatch)
 177             frac_min(ii,ibatch,iprot)=frac_min_set(ibatch)
 178             nc_fragm(ii,1,ibatch,iprot)=ncfrac_bet(ibatch)
 179             nc_req_setf(ii,1,ibatch,iprot)=ncreq_bet(ibatch)
 180             istruct(ii,ibatch,iprot)=2
 181           enddo
 182         enddo
 183         nfrag(1,iprot)=nfrag(1,iprot)+nhairp
 184         do ibatch=1,nclass(iprot)-1
 185           do i=1,nbfrag
 186             ii=i+nfrag(1,iprot)
 187             npiece(ii,1,ibatch,iprot)=2
 188             ifrag(1,1,ii,ibatch,iprot)=bfrag(1,i)
 189             ifrag(2,1,ii,ibatch,iprot)=bfrag(2,i)
 190             if (bfrag(3,i).lt.bfrag(4,i)) then
 191               ifrag(1,2,ii,ibatch,iprot)=bfrag(3,i)
 192               ifrag(2,2,ii,ibatch,iprot)=bfrag(4,i)
 193             else
 194               ifrag(1,2,ii,ibatch,iprot)=bfrag(4,i)
 195               ifrag(2,2,ii,ibatch,iprot)=bfrag(3,i)
 196             endif
 197             n_shift(1,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 198             n_shift(2,ii,1,ibatch,iprot)=nshift_bet(ibatch)
 199             ang_cut(ii,ibatch,iprot)=angcut_bet(ibatch)
 200             ang_cut1(ii,ibatch,iprot)=angcut1_bet(ibatch)
 201             frac_min(ii,ibatch,iprot)=frac_min_set(ibatch)
 202             nc_fragm(ii,1,ibatch,iprot)=ncfrac_bet(ibatch)
 203             nc_req_setf(ii,1,ibatch,iprot)=ncreq_bet(ibatch)
 204             istruct(ii,ibatch,iprot)=4
 205           enddo
 206         enddo
 207         nfrag(1,iprot)=nfrag(1,iprot)+nbfrag
 208       endif
 209       write (iout,*) "The following primary fragments were found:"
 210 #ifdef DEBUG
 211       write (iout,*) "Helices:",nhfrag
 212 #endif
 213       do i=1,nhfrag
 214         i1=ifrag(1,1,i,1,iprot)
 215         i2=ifrag(2,1,i,1,iprot)
 216         it1=itype(i1)
 217         it2=itype(i2)
 218 #ifdef DEBUG
 219         write (iout,'(i3,2x,a,i4,2x,a,i4)')
 220      &       i,restyp(it1),i1,restyp(it2),i2
 221 #endif
 222       enddo
 223 #ifdef DEBUG
 224       write (iout,*) "Hairpins:",nhairp
 225 #endif
 226       do i=nhfrag+1,nhfrag+nhairp
 227         i1=ifrag(1,1,i,1,iprot)
 228         i2=ifrag(2,1,i,1,iprot)
 229         it1=itype(i1)
 230         it2=itype(i2)
 231 #ifdef DEBUG
 232         write (iout,'(i3,2x,a,i4,2x,a,i4,2x)')
 233      &       i,restyp(it1),i1,restyp(it2),i2
 234 #endif
 235       enddo
 236 #ifdef DEBUG
 237       write (iout,*) "Far strand pairs:",nbfrag
 238 #endif
 239       do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
 240         i1=ifrag(1,1,i,1,iprot)
 241         i2=ifrag(2,1,i,1,iprot)
 242         it1=itype(i1)
 243         it2=itype(i2)
 244         i3=ifrag(1,2,i,1,iprot)
 245         i4=ifrag(2,2,i,1,iprot)
 246         it3=itype(i3)
 247         it4=itype(i4)
 248 #ifdef DEBUG
 249         write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)')
 250      &       i,restyp(it1),i1,restyp(it2),i2,
 251      &         restyp(it3),i3,restyp(it4),i4
 252 #endif
 253       enddo
 254 #ifdef DEBUG
 255       write (iout,*) "Strands:",nstrand
 256 #endif
 257       do i=nhfrag+nhairp+nbfrag+1,nfrag(1,iprot)
 258         i1=ifrag(1,1,i,1,iprot)
 259         i2=ifrag(2,1,i,1,iprot)
 260         it1=itype(i1)
 261         it2=itype(i2)
 262 #ifdef DEBUG
 263         write (iout,'(i3,2x,a,i4,2x,a,i4)')
 264      &       i,restyp(it1),i1,restyp(it2),i2
 265 #endif
 266       enddo
 267 #ifdef DEBUG
 268       do ibatch=1,nclass(iprot)-1
 269       write (iout,*) "Structural level",ibatch+1
 270       write (iout,*) "Fragments before sorting:"
 271       do i=1,nfrag(1)
 272         write (iout,'(20i4)') i,istruct(i,1,iprot),
 273      &    npiece(i,1,ibatch,iprot),
 274      &    (ifrag(1,j,i,1,iprot),ifrag(2,j,i,1,iprot),
 275      &    j=1,npiece(i,1,ibatch,iprot))
 276       enddo
 277       enddo
 278 #endif
 279       do ibatch=1,nclass(iprot)-1
 280       call lmysort(nfrag(1,iprot),2,maxpiece,ifrag(1,1,1,ibatch,iprot),
 281      &  npiece(1,1,ibatch,iprot),istruct(1,ibatch,iprot),
 282      &  n_shift(1,1,1,ibatch,iprot),ang_cut(1,ibatch,iprot),
 283      &  ang_cut1(1,ibatch,iprot),frac_min(1,ibatch,iprot),
 284      &  nc_fragm(1,1,ibatch,iprot),nc_req_setf(1,1,ibatch,iprot))
 285 #ifdef DEBUG
 286       write (iout,*) "Fragments after sorting:"
 287       do i=1,nfrag(1)
 288         write (iout,'(20i4)') i,istruct(i,1,iprot),
 289      &    npiece(i,1,ibatch,iprot),
 290      &    (ifrag(1,j,i,1,iprot),ifrag(2,j,i,1,iprot),
 291      &    j=1,npiece(i,1,ibatch,iprot))
 292       enddo
 293 #endif
 294       enddo
 295       call flush(iout)
 296       do i=1,nfrag(1,iprot)
 297         i1=ifrag(1,1,i,1,iprot)
 298         i2=ifrag(2,1,i,1,iprot)
 299         it1=itype(i1)
 300         it2=itype(i2)
 301 #ifdef DEBUG
 302         write (iout,*) i,ifrag(1,1,i,1,iprot),ifrag(2,1,i,1,iprot),
 303      &    it1,it2,ifrag(1,2,i,1,iprot),ifrag(2,2,i,1,iprot)
 304 #endif
 305         write (iout,'(i3,2x,a,i4,2x,a,i4,$)')
 306      &       i,restyp(it1),i1,restyp(it2),i2
 307         call flush(iout)
 308         if (npiece(i,1,1,iprot).eq.1) then
 309           write (iout,'(2x,a)') strstr(istruct(i,1,iprot))
 310         else
 311           i1=ifrag(1,2,i,1,iprot)
 312           i2=ifrag(2,2,i,1,iprot)
 313           it1=itype(i1)
 314           it2=itype(i2)
 315 #ifdef DEBUG
 316         write (iout,*) i,ifrag(1,2,i,1,iprot),ifrag(2,2,i,1,iprot),
 317      &    it1,it2
 318         call flush(iout)
 319 #endif
 320           write (iout,'(2x,a,i4,2x,a,i4,2x,a)')
 321      &       restyp(it1),i1,restyp(it2),i2,strstr(istruct(i,1,iprot))
 322         endif
 323       enddo
 324       return
 325       end
 326 c------------------------------------------------------------------------------
 327       subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
 328       implicit none
 329       include 'DIMENSIONS'
 330       include 'DIMENSIONS.ZSCOPT'
 331       include 'DIMENSIONS.COMPAR'
 332       include 'COMMON.IOUNITS'
 333       include 'COMMON.COMPAR'
 334       integer i,j,k
 335       integer nbfrag,bfrag(4,maxres/3)
 336       integer nhairp,ihairp(2,maxres/5)
 337 #ifdef DEBUG
 338       write (iout,*) "Entered find_and_remove_hairpins"
 339 #endif
 340       if (print_secondary) then
 341       write (iout,*) "Before checking for hairpins: nbfrag",nbfrag
 342       do i=1,nbfrag
 343         write (iout,*) i,(bfrag(k,i),k=1,4)
 344       enddo
 345       endif
 346       nhairp=0
 347       i=1
 348       do while (i.le.nbfrag)
 349         if (print_secondary)
 350      &  write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
 351         if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5)
 352      &  then
 353           if (print_secondary)
 354      &    write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
 355           nhairp=nhairp+1
 356           ihairp(1,nhairp)=bfrag(1,i)
 357           ihairp(2,nhairp)=bfrag(3,i)
 358           nbfrag=nbfrag-1
 359           do j=i,nbfrag
 360             do k=1,4
 361               bfrag(k,j)=bfrag(k,j+1)
 362             enddo
 363           enddo
 364         else
 365           i=i+1
 366         endif
 367       enddo
 368       if (print_secondary) then
 369       write (iout,*) "After finding hairpins:"
 370       write (iout,*) "nhairp",nhairp
 371       do i=1,nhairp
 372         write (iout,*) i,ihairp(1,i),ihairp(2,i)
 373       enddo
 374       write (iout,*) "nbfrag",nbfrag
 375       do i=1,nbfrag
 376         write (iout,*) i,(bfrag(k,i),k=1,4)
 377       enddo
 378       endif
 379       return
 380       end
 381 c------------------------------------------------------------------------------
 382       subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
 383       implicit none
 384       include 'DIMENSIONS'
 385       include 'DIMENSIONS.ZSCOPT'
 386       include 'DIMENSIONS.COMPAR'
 387       include 'COMMON.COMPAR'
 388       include 'COMMON.IOUNITS'
 389       integer i,k
 390       integer nbfrag,bfrag(4,maxres/3)
 391       integer nstrand,istrand(2,maxres/2)
 392       integer nhairp,ihairp(2,maxres/5)
 393       logical found
 394 #ifdef DEBUG
 395       write (iout,*) "Entered split_beta"
 396 #endif
 397       if (print_secondary) then
 398       write (iout,*) "Before splitting hairpins: nbfrag",nbfrag
 399       do i=1,nbfrag
 400         write (iout,*) i,(bfrag(k,i),k=1,4)
 401       enddo
 402       endif
 403       nstrand=0
 404       do i=1,nbfrag
 405         if (print_secondary)
 406      &  write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
 407         call add_strand(nstrand,istrand,nhairp,ihairp,
 408      &     bfrag(1,i),bfrag(2,i),found)
 409         if (bfrag(3,i).lt.bfrag(4,i)) then
 410           if (print_secondary)
 411      &    write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
 412           call add_strand(nstrand,istrand,nhairp,ihairp,
 413      &     bfrag(3,i),bfrag(4,i),found)
 414         else
 415           if (print_secondary)
 416      &    write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
 417           call add_strand(nstrand,istrand,nhairp,ihairp,
 418      &      bfrag(4,i),bfrag(3,i),found)
 419         endif
 420       enddo
 421       nbfrag=0
 422       if (print_secondary) then
 423       write (iout,*) "Strands found:",nstrand
 424       do i=1,nstrand
 425         write (iout,*) i,istrand(1,i),istrand(2,i)
 426       enddo
 427       endif
 428       return
 429       end
 430 c------------------------------------------------------------------------------
 431       subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
 432       implicit none
 433       include 'DIMENSIONS'
 434       include 'DIMENSIONS.ZSCOPT'
 435       include 'DIMENSIONS.COMPAR'
 436       include 'COMMON.IOUNITS'
 437       include 'COMMON.COMPAR'
 438       integer nstrand,istrand(2,maxres/2)
 439       integer nhairp,ihairp(2,maxres/5)
 440       integer is1,is2,j,idelt
 441       logical found
 442       found=.false.
 443       do j=1,nhairp
 444         idelt=(ihairp(2,j)-ihairp(1,j))/6
 445         if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
 446           if (print_secondary)
 447      &    write (iout,*) "strand",is1,is2," is part of hairpin",
 448      &      ihairp(1,j),ihairp(2,j)
 449           return
 450         endif
 451       enddo
 452       do j=1,nstrand
 453         idelt=(istrand(2,j)-istrand(1,j))/3
 454         if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt)
 455      &  then
 456 c The strand already exists in the array; update its ends if necessary.
 457           if (print_secondary)
 458      &    write (iout,*) "strand",is1,is2," found at position",j,
 459      &     ":",istrand(1,j),istrand(2,j)
 460           istrand(1,j)=min0(istrand(1,j),is1)
 461           istrand(2,j)=max0(istrand(2,j),is2)
 462           return
 463         endif
 464       enddo
 465 c The strand has not been found; add it to the array.
 466       if (print_secondary)
 467      &write (iout,*) "strand",is1,is2," added to the array."
 468       found=.true.
 469       nstrand=nstrand+1
 470       istrand(1,nstrand)=is1
 471       istrand(2,nstrand)=is2
 472       return
 473       end
 474 c------------------------------------------------------------------------------
 475       subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr,iprot)
 476       implicit none
 477       include 'DIMENSIONS'
 478       include 'DIMENSIONS.ZSCOPT'
 479       include 'COMMON.IOUNITS'
 480       include 'COMMON.FRAG'
 481       include 'COMMON.VAR'
 482       include 'COMMON.GEO'
 483       include 'COMMON.CHAIN'
 484       include 'COMMON.NAMES'
 485       include 'COMMON.INTERACT'
 486       integer i,j,ij,k,i1,j1,ii1,jj1,iii1,jjj1,ist,ien,nhelix,nbeta,
 487      &  nstrand,iprot,n310helix
 488       integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres),
 489      &  isecstr(maxres)
 490       logical lprint,lprint_sec,not_done,freeres
 491       double precision p1,p2
 492       external freeres
 493       character*1 csec(0:3) /'-','E','H','3'/
 494       if (lprint_sec) then
 495         write (iout,*) "entered SECONDARY2 iprot",iprot," ncont",ncont
 496         write (iout,*) "nstart_sup",nstart_sup(iprot),
 497      &    " nend_sup",nend_sup(iprot)
 498         write (iout,*) "The ICONT array"
 499         do i=1,ncont
 500           write (iout,*) icont(1,i),icont(2,i)
 501         enddo
 502         call flush(iout)
 503       endif
 504       do i=1,nres
 505         isecstr(i)=0
 506       enddo
 507       nbfrag=0
 508       nhfrag=0
 509       do i=1,nres
 510         isec(i,1)=0
 511         isec(i,2)=0
 512         nsec(i)=0
 513       enddo
 514
 515 c finding parallel beta
 516 cd      write (iout,*) '------- looking for parallel beta -----------'
 517       nbeta=0
 518       nstrand=0
 519       do i=1,ncont
 520         i1=icont(1,i)
 521         j1=icont(2,i)
 522         if (i1.ge.nstart_sup(iprot) .and. i1.le.nend_sup(iprot)
 523      &     .and. j1.gt.nstart_sup(iprot) .and. j1.le.nend_sup(iprot))
 524      &    then
 525 cd        write (iout,*) "parallel",i1,j1
 526         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
 527           ii1=i1
 528           jj1=j1
 529 cd          write (iout,*) i1,j1
 530           not_done=.true.
 531           do while (not_done)
 532            i1=i1+1
 533            j1=j1+1
 534             do j=1,ncont
 535               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
 536      &             freeres(i1,j1,nsec,isec)) goto 5
 537             enddo
 538             not_done=.false.
 539   5         continue
 540 cd            write (iout,*) i1,j1,not_done
 541           enddo
 542           j1=j1-1
 543           i1=i1-1
 544           if (i1-ii1.gt.1) then
 545             ii1=max0(ii1-1,1)
 546             jj1=max0(jj1-1,1)
 547             nbeta=nbeta+1
 548             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
 549      &               nbeta,ii1,i1,jj1,j1
 550
 551             nbfrag=nbfrag+1
 552             bfrag(1,nbfrag)=ii1+1
 553             bfrag(2,nbfrag)=i1+1
 554             bfrag(3,nbfrag)=jj1+1
 555             bfrag(4,nbfrag)=min0(j1+1,nres)
 556
 557             do ij=ii1,i1
 558              nsec(ij)=nsec(ij)+1
 559              isec(ij,nsec(ij))=nbeta
 560             enddo
 561             do ij=jj1,j1
 562              nsec(ij)=nsec(ij)+1
 563              isec(ij,nsec(ij))=nbeta
 564             enddo
 565
 566            if(lprint_sec) then
 567             nstrand=nstrand+1
 568             if (nbeta.le.9) then
 569               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 570      &          "DefPropRes 'strand",nstrand,
 571      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 572             else
 573               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 574      &          "DefPropRes 'strand",nstrand,
 575      &          "' 'num = ",ii1-1,"..",i1-1,"'"
 576             endif
 577             nstrand=nstrand+1
 578             if (nbeta.le.9) then
 579               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 580      &          "DefPropRes 'strand",nstrand,
 581      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 582             else
 583               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 584      &          "DefPropRes 'strand",nstrand,
 585      &          "' 'num = ",jj1-1,"..",j1-1,"'"
 586             endif
 587               write(12,'(a8,4i4)')
 588      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
 589            endif
 590           endif
 591         endif
 592         endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
 593       enddo
 594
 595 c finding antiparallel beta
 596 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
 597
 598       do i=1,ncont
 599         i1=icont(1,i)
 600         j1=icont(2,i)
 601         if (freeres(i1,j1,nsec,isec)) then
 602           ii1=i1
 603           jj1=j1
 604 cd          write (iout,*) i1,j1
 605
 606           not_done=.true.
 607           do while (not_done)
 608            i1=i1+1
 609            j1=j1-1
 610             do j=1,ncont
 611               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
 612      &             freeres(i1,j1,nsec,isec)) goto 6
 613             enddo
 614             not_done=.false.
 615   6         continue
 616 cd            write (iout,*) i1,j1,not_done
 617           enddo
 618           i1=i1-1
 619           j1=j1+1
 620           if (i1-ii1.gt.1) then
 621
 622             nbfrag=nbfrag+1
 623             bfrag(1,nbfrag)=ii1
 624             bfrag(2,nbfrag)=min0(i1+1,nres)
 625             bfrag(3,nbfrag)=min0(jj1+1,nres)
 626             bfrag(4,nbfrag)=j1
 627
 628             nbeta=nbeta+1
 629             iii1=max0(ii1-1,1)
 630             do ij=iii1,i1
 631              nsec(ij)=nsec(ij)+1
 632              if (nsec(ij).le.2) then
 633               isec(ij,nsec(ij))=nbeta
 634              endif
 635             enddo
 636             jjj1=max0(j1-1,1)
 637             do ij=jjj1,jj1
 638              nsec(ij)=nsec(ij)+1
 639              if (nsec(ij).le.2) then
 640               isec(ij,nsec(ij))=nbeta
 641              endif
 642             enddo
 643
 644
 645            if (lprint_sec) then
 646             write (iout,'(a,i3,4i4)')'antiparallel beta',
 647      &                   nbeta,ii1-1,i1,jj1,j1-1
 648             nstrand=nstrand+1
 649             if (nstrand.le.9) then
 650               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 651      &          "DefPropRes 'strand",nstrand,
 652      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 653             else
 654               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 655      &          "DefPropRes 'strand",nstrand,
 656      &          "' 'num = ",ii1-2,"..",i1-1,"'"
 657             endif
 658             nstrand=nstrand+1
 659             if (nstrand.le.9) then
 660               write(12,'(a18,i1,a9,i3,a2,i3,a1)')
 661      &          "DefPropRes 'strand",nstrand,
 662      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 663             else
 664               write(12,'(a18,i2,a9,i3,a2,i3,a1)')
 665      &          "DefPropRes 'strand",nstrand,
 666      &          "' 'num = ",j1-2,"..",jj1-1,"'"
 667             endif
 668               write(12,'(a8,4i4)')
 669      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
 670            endif
 671           endif
 672         endif
 673       enddo
 674
 675 cd      write (iout,*) "After beta:",nbfrag
 676 cd      do i=1,nbfrag
 677 cd        write (iout,*) (bfrag(j,i),j=1,4)
 678 cd      enddo
 679
 680       if (nstrand.gt.0.and.lprint_sec) then
 681         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
 682         do i=2,nstrand
 683          if (i.le.9) then
 684           write(12,'(a9,i1,$)') " | strand",i
 685          else
 686           write(12,'(a9,i2,$)') " | strand",i
 687          endif
 688         enddo
 689         write(12,'(a1)') "'"
 690       endif
 691
 692
 693 c finding alpha or 310 helix
 694
 695       nhelix=0
 696       do i=1,ncont
 697         i1=icont(1,i)
 698         j1=icont(2,i)
 699         p1=phi(i1+2)*rad2deg
 700         p2=0.0
 701         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
 702
 703
 704         if (j1.eq.i1+3 .and.
 705      &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
 706      &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
 707 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
 708 co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
 709           ii1=i1
 710           jj1=j1
 711           if (nsec(ii1).eq.0) then
 712             not_done=.true.
 713           else
 714             not_done=.false.
 715           endif
 716           do while (not_done)
 717             i1=i1+1
 718             j1=j1+1
 719             do j=1,ncont
 720               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 721             enddo
 722             not_done=.false.
 723   10        continue
 724             p1=phi(i1+2)*rad2deg
 725             p2=phi(j1+2)*rad2deg
 726             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80)
 727      &                              not_done=.false.
 728
 729 cd           write (iout,*) i1,j1,not_done,p1,p2
 730           enddo
 731           j1=j1+1
 732           if (j1-ii1.gt.4) then
 733             nhelix=nhelix+1
 734 cd            write (iout,*)'helix',nhelix,ii1,j1
 735
 736             nhfrag=nhfrag+1
 737             hfrag(1,nhfrag)=ii1
 738             hfrag(2,nhfrag)=j1
 739
 740             do ij=ii1,j1
 741              nsec(ij)=-1
 742             enddo
 743            if (lprint_sec) then
 744             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
 745             if (nhelix.le.9) then
 746               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
 747      &          "DefPropRes 'helix",nhelix,
 748      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 749             else
 750               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
 751      &          "DefPropRes 'helix",nhelix,
 752      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 753             endif
 754            endif
 755           endif
 756         endif
 757       enddo
 758
 759       if (nhelix.gt.0.and.lprint_sec) then
 760         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
 761         do i=2,nhelix
 762          if (nhelix.le.9) then
 763           write(12,'(a8,i1,$)') " | helix",i
 764          else
 765           write(12,'(a8,i2,$)') " | helix",i
 766          endif
 767         enddo
 768         write(12,'(a1)') "'"
 769       endif
 770
 771       if (lprint_sec) then
 772        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
 773        write(12,'(a20)') "XMacStand ribbon.mac"
 774       endif
 775
 776       if (lprint) then
 777
 778         write(iout,*) 'UNRES seq:'
 779         do j=1,nbfrag
 780          write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
 781         enddo
 782
 783         do j=1,nhfrag
 784          write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
 785         enddo
 786
 787       endif
 788
 789       do j=1,nbfrag
 790         do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j))
 791           isecstr(k)=1
 792         enddo
 793         do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j))
 794           isecstr(k)=1
 795         enddo
 796       enddo
 797       do j=1,nhfrag
 798         do k=hfrag(1,j),hfrag(2,j)
 799           isecstr(k)=2
 800         enddo
 801       enddo
 802 #ifdef HELIX
 803 c Find 3-10 helices
 804 c      write (iout,*) "Finding 3-10 helices"
 805       n310helix=0
 806       nh310frag=0
 807       do i=1,ncont
 808         i1=icont(1,i)
 809         j1=icont(2,i)
 810 c        write (iout,*) "i",i," i1",i1,isecstr(i1)," j1",j1,isecstr(j1)
 811         if (isecstr(i1).eq.0 .and. isecstr(i1).eq.0 .and.
 812      &     j1.eq.i1+2) then
 813           ii1=i1
 814           jj1=j1
 815           if (nsec(ii1).eq.0) then
 816             not_done=.true.
 817           else
 818             not_done=.false.
 819           endif
 820           do while (not_done)
 821             i1=i1+1
 822             j1=j1+1
 823             do j=1,ncont
 824               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 12
 825             enddo
 826             not_done=.false.
 827   12        continue
 828 c            write (iout,*) "d i1",isecstr(i1),i1," j1",j1,isecstr(i1)
 829             if (isecstr(i1).ne.0 .or. isecstr(j1).ne.0) then
 830               j1=j1-1
 831               i1=i1-1
 832               not_done=.false.
 833             endif
 834 c            write (iout,*) "not_done",not_done
 835           enddo
 836           j1=j1+1
 837 c          write (iout,*) "ii1",ii1," j1",j1
 838           if (j1-ii1.gt.4) then
 839             n310helix=n310helix+1
 840             nh310frag=nh310frag+1
 841             h310frag(1,nh310frag)=ii1
 842             h310frag(2,nh310frag)=j1
 843
 844             do ij=ii1,j1
 845              nsec(ij)=-1
 846             enddo
 847 c        write (iout,*) "n310helix",n310helix," nh310frag",nh310frag
 848            if (lprint_sec) then
 849             write (iout,'(a,i3,2i4)')
 850      &       "3-10 helix",n310helix,ii1-1,j1-1
 851             if (n310helix.le.9) then
 852               write(12,'(a17,i1,a9,i3,a2,i3,a1)')
 853      &          "DefPropRes 'helix",n310helix,
 854      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 855             else
 856               write(12,'(a17,i2,a9,i3,a2,i3,a1)')
 857      &          "DefPropRes 'helix",n310helix,
 858      &          "' 'num = ",ii1-1,"..",j1-2,"'"
 859             endif
 860            endif
 861           endif
 862         endif
 863 c        write (iout,*) "n310helix",n310helix," nh310frag",nh310frag
 864       enddo
 865 #endif
 866 c      write (iout,*) "nh310frag",nh310frag
 867       do j=1,nh310frag
 868 c        write (iout,*) "h310frag",h310frag(1,j),h310frag(2,j)
 869         do k=h310frag(1,j),h310frag(2,j)
 870           isecstr(k)=3
 871         enddo
 872       enddo
 873       if (lprint) then
 874         write (iout,*)
 875         write (iout,*) "Secondary structure"
 876         do i=1,nres,80
 877           ist=i
 878           ien=min0(i+79,nres)
 879           write (iout,*)
 880           write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
 881           write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien)
 882           write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien)
 883         enddo
 884         write (iout,*)
 885       endif
 886       return
 887       end
 888 c-------------------------------------------------
 889       logical function freeres(i,j,nsec,isec)
 890       implicit none
 891       include 'DIMENSIONS'
 892       integer i,j,k,l
 893       integer isec(maxres,4),nsec(maxres)
 894       freeres=.false.
 895
 896       if (nsec(i).gt.1.or.nsec(j).gt.1) return
 897       do k=1,nsec(i)
 898         do l=1,nsec(j)
 899           if (isec(i,k).eq.isec(j,l)) return
 900         enddo
 901       enddo
 902       freeres=.true.
 903       return
 904       end
 905