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