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