source/wham/src-M/conf_compar.F

   1       subroutine conf_compar(jcon,lprn,print_class)
   2       implicit real*8 (a-h,o-z)
   3 #ifdef MPI
   4       include "mpif.h"
   5 #endif
   6       include 'DIMENSIONS'
   7       include 'DIMENSIONS.ZSCOPT'
   8       include 'DIMENSIONS.COMPAR'
   9       include 'DIMENSIONS.FREE'
  10       include 'COMMON.CONTROL'
  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.FREE'
  20       include 'COMMON.ENERGIES'
  21 #ifdef MPI
  22       include 'COMMON.MPI'
  23 #endif
  24       integer ilen
  25       external ilen
  26       logical lprn,print_class
  27       integer ncont_frag(mmaxfrag),
  28      & icont_frag(2,maxcont,mmaxfrag),ncontsc,
  29      & icontsc(1,maxcont),nsccont_frag(mmaxfrag),
  30      & isccont_frag(2,maxcont,mmaxfrag)
  31       integer isecstr(maxres)
  32       integer itemp(maxfrag)
  33       character*4 liczba
  34       double precision Epot
  35 c      print *,"Enter conf_compar",jcon
  36       call angnorm12(rmsang)
  37 c Level 1: check secondary and supersecondary structure
  38       call elecont(lprn,ncont,icont,nnt,nct)
  39       if (lprn) then
  40         write (iout,*) "elecont finished"
  41         call flush(iout)
  42       endif
  43       call secondary2(lprn,.false.,ncont,icont,isecstr)
  44       if (lprn) then
  45         write (iout,*) "secondary2 finished"
  46         call flush(iout)
  47       endif
  48       call contact(lprn,ncontsc,icontsc,nnt,nct)
  49       if (lprn) then
  50          write(iout,*) "Assigning electrostatic contacts"
  51          call flush(iout)
  52       endif
  53       call contacts_between_fragments(lprn,3,ncont,icont,ncont_frag,
  54      &   icont_frag)
  55       if (lprn) then
  56         write(iout,*) "Assigning sidechain contacts"
  57         call flush(iout)
  58       endif
  59       call contacts_between_fragments(lprn,3,ncontsc,icontsc,
  60      &   nsccont_frag,isccont_frag)
  61       if (lprn) then
  62         write(iout,*) "--> After contacts_between_fragments"
  63         call flush(iout)
  64       endif
  65       do i=1,nlevel
  66         do j=1,isnfrag(nlevel+1)
  67           iclass(j,i)=0
  68         enddo
  69       enddo
  70       do j=1,nfrag(1)
  71         ind = icant(j,j)
  72         if (lprn) then
  73           write (iout,'(80(1h=))')
  74           write (iout,*) "Level",1," fragment",j
  75           write (iout,'(80(1h=))')
  76         endif
  77         call flush(iout)
  78         rmsfrag(j,1)=rmscalc(0,1,j,jcon,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         call flush(iout)
  84         ncnat=ncont_frag_ref(ind)
  85 c        write (iout,*) "before match_contact:",nc_fragm(j,1),
  86 c     &   nc_req_setf(j,1)
  87 c        call flush(iout)
  88         call match_secondary(j,isecstr,nsec_match,lprn)
  89         if (lprn) write (iout,*) "Fragment",j," nsec_match",
  90      &    nsec_match," length",len_frag(j,1)," min_len",
  91      &    frac_sec*len_frag(j,1)
  92         if (nsec_match.lt.frac_sec*len_frag(j,1)) then
  93           iclass(j,1)=0
  94           if (lprn) write (iout,*) "Fragment",j,
  95      &      " has incorrect secondary structure"
  96         else
  97           iclass(j,1)=1
  98           if (lprn) write (iout,*) "Fragment",j,
  99      &      " has correct secondary structure"
 100         endif
 101         if (ielecont(j,1).gt.0) then
 102           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 103      &      ncont_frag_ref(ind),icont_frag_ref(1,1,ind),
 104      &      ncont_frag(ind),icont_frag(1,1,ind),
 105      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
 106      &      nc_req_setf(j,1),istruct(j),.true.,lprn)
 107         else if (isccont(j,1).gt.0) then
 108           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 109      &      nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),
 110      &      nsccont_frag(ind),isccont_frag(1,1,ind),
 111      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
 112      &      nc_req_setf(j,1),istruct(j),.true.,lprn)
 113         else if (iloc(j).gt.0) then
 114 c          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
 115           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 116      &      0,icont_frag_ref(1,1,ind),
 117      &      ncont_frag(ind),icont_frag(1,1,ind),
 118      &      j,n_shift(1,j,1),n_shift(2,j,1),nc_fragm(j,1),
 119      &      0,istruct(j),.true.,lprn)
 120 c          write (iout,*) "n_shif",n_shift(1,j,1),n_shift(2,j,1)
 121         else
 122           ishif=0
 123           nc_match=1
 124         endif
 125         if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2
 126         ishif=ishif1
 127         qfrag(j,1)=qwolynes(1,j)
 128         if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 129         if (lprn) write (iout,*) "ishift",ishif," nc_match",nc_match
 130 c        write (iout,*) "j",j," ishif",ishif," rms",rmsfrag(j,1)
 131         if (irms(j,1).gt.0) then
 132           if (rmsfrag(j,1).le.rmscutfrag(1,j,1)) then
 133             iclass_rms=2
 134             ishifft_rms=0
 135           else
 136             ishiff=0
 137             rms=1.0d2
 138             iclass_rms=0
 139             do while (rms.gt.rmscutfrag(1,j,1) .and.
 140      &         ishiff.lt.n_shift(1,j,1))
 141               ishiff=ishiff+1
 142               rms=rmscalc(-ishiff,1,j,jcon,lprn)
 143 c              write(iout,*)"jcon,i,j,ishiff",jcon,i,j,-ishiff,
 144 c     &          " rms",rms," rmscut",rmscutfrag(1,j,1)
 145               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 146               if (rms.gt.rmscutfrag(1,j,1)) then
 147                 rms=rmscalc(ishiff,1,j,jcon,lprn)
 148 c                write (iout,*) "jcon,1,j,ishiff",jcon,1,j,ishiff,
 149 c     &           " rms",rms
 150               endif
 151               if (lprn) write (iout,*) "rms",rmsfrag(j,1)
 152             enddo
 153 c            write (iout,*) "After loop: rms",rms,
 154 c     &        " rmscut",rmscutfrag(1,j,1)
 155 c            write (iout,*) "iclass_rms",iclass_rms
 156             if (rms.le.rmscutfrag(1,j,1)) then
 157               ishifft_rms=ishiff
 158               rmsfrag(j,1)=rms
 159               iclass_rms=1
 160             endif
 161 c            write (iout,*) "iclass_rms",iclass_rms
 162           endif
 163 c          write (iout,*) "ishif",ishif
 164           if (iabs(ishifft_rms).gt.iabs(ishif)) ishif=ishifft_rms
 165         else
 166           iclass_rms=1
 167         endif
 168 c        write (iout,*) "ishif",ishif," iclass",iclass(j,1),
 169 c     &    " iclass_rms",iclass_rms
 170         if (nc_match.gt.0 .and. iclass_rms.gt.0) then
 171           if (ishif.eq.0) then
 172             iclass(j,1)=iclass(j,1)+6
 173           else
 174             iclass(j,1)=iclass(j,1)+2
 175           endif
 176         endif
 177         ncont_nat(1,j,1)=nc_match
 178         ncont_nat(2,j,1)=ncon_match
 179         ishifft(j,1)=ishif
 180 c        write (iout,*) "iclass",iclass(j,1)
 181       enddo
 182 c Next levels: Check arrangements of elementary fragments.
 183       do i=2,nlevel
 184         do j=1,nfrag(i)
 185         if (i .eq. 2) ind = icant(ipiece(1,j,i),ipiece(2,j,i))
 186         if (lprn) then
 187             write (iout,'(80(1h=))')
 188             write (iout,*) "Level",i," fragment",j
 189             write (iout,'(80(1h=))')
 190         endif
 191 c If an elementary fragment doesn't exist, don't check higher hierarchy levels.
 192         do k=1,npiece(j,i)
 193           ik=ipiece(k,j,i)
 194           if (iclass(ik,1).eq.0) then
 195             iclass(j,i)=0
 196             goto 12
 197           endif
 198         enddo
 199         if (i.eq.2 .and. ielecont(j,i).gt.0) then
 200           iclass_con=0
 201           ishifft_con=0
 202           if (lprn) write (iout,*)
 203      &     "Comparing electrostatic contact map: fragments",
 204      &      ipiece(1,j,i),ipiece(2,j,i)," ind",ind
 205           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 206      &     ncont_frag_ref(ind),icont_frag_ref(1,1,ind),
 207      &     ncont_frag(ind),icont_frag(1,1,ind),
 208      &     j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),
 209      &     nc_req_setf(j,i),2,.false.,lprn)
 210           ishif=ishif1
 211           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 212           if (nc_match.gt.0) then
 213             if (ishif.eq.0) then
 214               iclass_con=2
 215             else
 216               iclass_con=1
 217             endif
 218           endif
 219           ncont_nat(1,j,i)=nc_match
 220           ncont_nat(2,j,i)=ncon_match
 221           ishifft_con=ishif
 222         else if (i.eq.2 .and. isccont(j,i).gt.0) then
 223           iclass_con=0
 224           ishifft_con=0
 225           if (lprn) write (iout,*)
 226      &     "Comparing sidechain contact map: fragments",
 227      &    ipiece(1,j,i),ipiece(2,j,i)," ind",ind
 228           call match_contact(ishif1,ishif2,nc_match,ncon_match,
 229      &     nsccont_frag_ref(ind),isccont_frag_ref(1,1,ind),
 230      &     nsccont_frag(ind),isccont_frag(1,1,ind),
 231      &     j,n_shift(1,j,i),n_shift(2,j,i),nc_fragm(j,i),
 232      &     nc_req_setf(j,i),2,.false.,lprn)
 233           ishif=ishif1
 234           if (iabs(ishif2).gt.iabs(ishif1)) ishif=ishif2
 235           if (nc_match.gt.0) then
 236             if (ishif.eq.0) then
 237               iclass_con=2
 238             else
 239               iclass_con=1
 240             endif
 241           endif
 242           ncont_nat(1,j,i)=nc_match
 243           ncont_nat(2,j,i)=ncon_match
 244           ishifft_con=ishif
 245         else if (i.eq.2) then
 246           iclass_con=2
 247           ishifft_con=0
 248         endif
 249         if (i.eq.2) qfrag(j,2)=qwolynes(2,j)
 250         if (lprn) write (iout,*)
 251      &    "Comparing rms: fragments",
 252      &     (ipiece(k,j,i),k=1,npiece(j,i))
 253         rmsfrag(j,i)=rmscalc(0,i,j,jcon,lprn)
 254         if (irms(j,i).gt.0) then
 255           iclass_rms=0
 256           ishifft_rms=0
 257           if (lprn) write (iout,*) "rms",rmsfrag(j,i)
 258 c          write (iout,*) "i",i," j",j," rmsfrag",rmsfrag(j,i),
 259 c     &     " rmscutfrag",rmscutfrag(1,j,i)
 260           if (rmsfrag(j,i).le.rmscutfrag(1,j,i)) then
 261             iclass_rms=2
 262             ishifft_rms=0
 263           else
 264             ishif=0
 265             rms=1.0d2
 266             do while (rms.gt.rmscutfrag(1,j,i) .and.
 267      &         ishif.lt.n_shift(1,j,i))
 268               ishif=ishif+1
 269               rms=rmscalc(-ishif,i,j,jcon,lprn)
 270 c              print *,"jcon,i,j,ishif",jcon,i,j,-ishif," rms",rms
 271               if (lprn) write (iout,*) "rms",rmsfrag(j,i)
 272               if (rms.gt.rmscutfrag(1,j,i)) then
 273                 rms=rmscalc(ishif,i,j,jcon,lprn)
 274 c                print *,"jcon,i,j,ishif",jcon,i,j,ishif," rms",rms
 275               endif
 276               if (lprn) write (iout,*) "rms",rms
 277             enddo
 278             if (rms.le.rmscutfrag(1,j,i)) then
 279               ishifft_rms=ishif
 280               rmsfrag(j,i)=rms
 281               iclass_rms=1
 282             endif
 283           endif
 284         endif
 285         if (irms(j,i).eq.0 .and. ielecont(j,i).eq.0 .and.
 286      &    isccont(j,i).eq.0 ) then
 287           write (iout,*) "Error: no measure of comparison specified:",
 288      &      " level",i," part",j
 289           stop
 290         endif
 291         if (lprn)
 292      &  write (iout,*) "iclass_con",iclass_con," iclass_rms",iclass_rms
 293         if (i.eq.2) then
 294           iclass(j,i) = min0(iclass_con,iclass_rms)
 295           if (iabs(ishifft_rms).gt.iabs(ishifft_con)) then
 296             ishifft(j,i)=ishifft_rms
 297           else
 298             ishifft(j,i)=ishifft_con
 299           endif
 300         else if (i.gt.2) then
 301           iclass(j,i) = iclass_rms
 302           ishifft(j,i)= ishifft_rms
 303         endif
 304    12   continue
 305         enddo
 306       enddo
 307       rms_nat=rmsnat(jcon)
 308       qnat=qwolynes(0,0)
 309 C Compute the structural class
 310       iscor=0
 311       IF (.NOT. BINARY) THEN
 312       do i=1,nlevel
 313         IF (I.EQ.1) THEN
 314         do j=1,nfrag(i)
 315           itemp(j)=iclass(j,i)
 316         enddo
 317         do kk=-1,1
 318           do j=1,nfrag(i)
 319             idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-kk*nfrag(i)-j
 320             iex = 2**idig
 321             im=mod(itemp(j),2)
 322             itemp(j)=itemp(j)/2
 323 c            write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 324 c     &        " iclass",iclass(j,i)," im",im
 325             iscor=iscor+im*iex
 326           enddo
 327         enddo
 328         ELSE
 329         do j=1,nfrag(i)
 330           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-j
 331           iex = 2**idig
 332           if (iclass(j,i).gt.0) then
 333             im=1
 334           else
 335             im=0
 336           endif
 337 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 338 c     &      " iclass",iclass(j,i)," im",im
 339           iscor=iscor+im*iex
 340         enddo
 341         do j=1,nfrag(i)
 342           idig = 2*isnfrag(nlevel+1)-2*isnfrag(i)-nfrag(i)-j
 343           iex = 2**idig
 344           if (iclass(j,i).gt.1) then
 345             im=1
 346           else
 347             im=0
 348           endif
 349 c          write (iout,*) "i",i," j",j," idig",idig," iex",iex,
 350 c     &      " iclass",iclass(j,i)," im",im
 351           iscor=iscor+im*iex
 352         enddo
 353         ENDIF
 354       enddo
 355       iscore=iscor
 356       ENDIF
 357       if (print_class) then
 358 #ifdef MPI
 359           write(istat,'(i6,$)') jcon+indstart(me)-1
 360           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),
 361      &     -entfac(jcon)
 362 #else
 363           write(istat,'(i6,$)') jcon
 364           write (istat,'(f10.2,$)') (potE(jcon,k),k=1,nParmSet),
 365      &      -entfac(jcon)
 366 #endif
 367           write (istat,'(f8.3,2f6.3,$)')
 368      &      rms_nat,qnat,rmsang/(nres-3)
 369           do j=1,nlevel
 370             write(istat,'(1x,$,20(i3,$))')
 371      &        (ncont_nat(1,k,j),k=1,nfrag(j))
 372             if (j.lt.3) then
 373               write(istat,'(1x,$,20(f5.1,f5.2$))')
 374      &          (rmsfrag(k,j),qfrag(k,j),k=1,nfrag(j))
 375             else
 376               write(istat,'(1x,$,20(f5.1$))')
 377      &          (rmsfrag(k,j),k=1,nfrag(j))
 378             endif
 379             write(istat,'(1x,$,20(i1,$))')
 380      &        (iclass(k,j),k=1,nfrag(j))
 381           enddo
 382           if (binary) then
 383             write (istat,'("  ",$)')
 384             do j=1,nlevel
 385               write (istat,'(100(i1,$))')(iclass(k,j),
 386      &           k=1,nfrag(j))
 387               if (j.lt.nlevel) write(iout,'(".",$)')
 388             enddo
 389             write (istat,*)
 390           else
 391             write (istat,'(i10)') iscore
 392           endif
 393       endif
 394       RETURN
 395       END