source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/qwolynes.f

   1 c---------------------------------------------------------------------------------------------------------
   2       double precision function qwolynes(ilevel,jfrag,ib,iprot)
   3       implicit none
   4       include 'DIMENSIONS'
   5       include 'DIMENSIONS.ZSCOPT'
   6       include 'DIMENSIONS.COMPAR'
   7       include 'COMMON.IOUNITS'
   8       include 'COMMON.COMPAR'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.INTERACT'
  11       include 'COMMON.VAR'
  12       include 'COMMON.PEPTCONT'
  13       include 'COMMON.CONTACTS1'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.CLASSES'
  16       integer ilevel,jfrag,ib,iprot
  17       integer i,j,jl,k,l,il,kl,nl,np,ip,kp
  18       integer nsep /3/
  19       double precision dist
  20       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  21       logical lprn /.false./
  22       double precision sigm,x
  23       sigm(x)=0.25d0*x
  24       if (lprn) write (iout,*) "QWolyes: iprot",iprot," ib",ib,
  25      &      " jfrag",jfrag," ilevel",ilevel
  26       qq = 0.0d0
  27       if (ilevel.eq.0) then
  28         if (lprn) write (iout,*) "Q computed for whole molecule"
  29         nl=0
  30         do il=nnt+nsep,nct
  31           do jl=nnt,il-nsep
  32             dij=0.0d0
  33             dijCM=0.0d0
  34             d0ij=0.0d0
  35             d0ijCM=0.0d0
  36             qqij=0.0d0
  37             qqijCM=0.0d0
  38             nl=nl+1
  39             d0ij=dsqrt((cref(1,jl,iprot)-cref(1,il,iprot))**2+
  40      &                 (cref(2,jl,iprot)-cref(2,il,iprot))**2+
  41      &                 (cref(3,jl,iprot)-cref(3,il,iprot))**2)
  42             dij=dist(il,jl)
  43             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  44             if (itype(il).ne.10 .or. itype(jl).ne.10) then
  45               nl=nl+1
  46               d0ijCM=dsqrt(
  47      &               (cref(1,jl+nres,iprot)-cref(1,il+nres,iprot))**2+
  48      &               (cref(2,jl+nres,iprot)-cref(2,il+nres,iprot))**2+
  49      &               (cref(3,jl+nres,iprot)-cref(3,il+nres,iprot))**2)
  50               dijCM=dist(il+nres,jl+nres)
  51               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  52             endif
  53             qq = qq+qqij+qqijCM
  54 c            if (lprn) then
  55 c              write (iout,*) "il",il," jl",jl,
  56 c     &         " itype",itype(il),itype(jl)
  57 c              write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
  58 c     &         " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
  59 c            endif
  60           enddo
  61         enddo
  62         qq = qq/nl
  63         if (lprn) write (iout,*) "nl",nl," qq",qq
  64       else if (ilevel.eq.1) then
  65         if (lprn) then
  66           write (iout,*) "Level",ilevel," fragment",jfrag
  67           write (iout,*) "nlist_frag",nlist_frag(jfrag,ib,iprot),
  68      &      " list_frag",(list_frag(j,jfrag,ib,iprot),
  69      &      j=1,nlist_frag(jfrag,ib,iprot))
  70         endif
  71         nl=0
  72         do i=2,nlist_frag(jfrag,ib,iprot)
  73           do j=1,i-1
  74             il=list_frag(i,jfrag,ib,iprot)
  75             jl=list_frag(j,jfrag,ib,iprot)
  76             if (iabs(il-jl).gt.nsep) then
  77               dij=0.0d0
  78               dijCM=0.0d0
  79               d0ij=0.0d0
  80               d0ijCM=0.0d0
  81               qqij=0.0d0
  82               qqijCM=0.0d0
  83               nl=nl+1
  84               d0ij=dsqrt((cref(1,jl,iprot)-cref(1,il,iprot))**2+
  85      &                   (cref(2,jl,iprot)-cref(2,il,iprot))**2+
  86      &                   (cref(3,jl,iprot)-cref(3,il,iprot))**2)
  87               dij=dist(il,jl)
  88               qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  89               if (itype(il).ne.10 .or. itype(jl).ne.10) then
  90                 nl=nl+1
  91                 d0ijCM=dsqrt(
  92      &                 (cref(1,jl+nres,iprot)-cref(1,il+nres,iprot))**2+
  93      &                 (cref(2,jl+nres,iprot)-cref(2,il+nres,iprot))**2+
  94      &                 (cref(3,jl+nres,iprot)-cref(3,il+nres,iprot))**2)
  95                 dijCM=dist(il+nres,jl+nres)
  96                qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  97               endif
  98               qq = qq+qqij+qqijCM
  99 c              if (lprn) then
 100 c                write (iout,*) "i",i," j",j," il",il," jl",jl,
 101 c     &           " itype",itype(il),itype(jl)
 102 c                write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
 103 c     &           " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
 104 c              endif
 105             endif
 106           enddo
 107         enddo
 108         qq = qq/nl
 109         if (lprn) write (iout,*) "nl",nl," qq",qq
 110       else if (ilevel.eq.2) then
 111         np=npiece(jfrag,ilevel,ib,iprot)
 112         nl=0
 113         if (lprn) then
 114           write (iout,*) "npiece",npiece(jfrag,ilevel,ib,iprot),
 115      &     " ipiece",(ipiece(i,jfrag,ilevel,ib,iprot),
 116      &       i=1,npiece(jfrag,ilevel,ib,iprot))
 117         endif
 118         do i=2,np
 119           ip=ipiece(i,jfrag,ilevel,ib,iprot)
 120           do j=1,nlist_frag(ip,ib,iprot)
 121             il=list_frag(j,ip,ib,iprot)
 122             do k=1,i-1
 123               kp=ipiece(k,jfrag,ilevel,ib,iprot)
 124               do l=1,nlist_frag(kp,ib,iprot)
 125                 kl=list_frag(l,kp,ib,iprot)
 126                 if (iabs(kl-il).gt.nsep) then
 127                   nl=nl+1
 128                   dij=0.0d0
 129                   dijCM=0.0d0
 130                   d0ij=0.0d0
 131                   d0ijCM=0.0d0
 132                   qqij=0.0d0
 133                   qqijCM=0.0d0
 134                   d0ij=dsqrt((cref(1,kl,iprot)-cref(1,il,iprot))**2+
 135      &                       (cref(2,kl,iprot)-cref(2,il,iprot))**2+
 136      &                       (cref(3,kl,iprot)-cref(3,il,iprot))**2)
 137                   dij=dist(il,kl)
 138                   qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
 139                   if (itype(il).ne.10 .or. itype(kl).ne.10) then
 140                     nl=nl+1
 141                     d0ijCM=dsqrt(
 142      &                 (cref(1,kl+nres,iprot)-cref(1,il+nres,iprot))**2+
 143      &                 (cref(2,kl+nres,iprot)-cref(2,il+nres,iprot))**2+
 144      &                 (cref(3,kl+nres,iprot)-cref(3,il+nres,iprot))**2)
 145                     dijCM=dist(il+nres,kl+nres)
 146                     qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/
 147      &                (sigm(d0ijCM)))**2)
 148                   endif
 149                   qq = qq+qqij+qqijCM
 150 c                  if (lprn) then
 151 c                    write (iout,*) "i",i," j",j," k",k," l",l," il",il,
 152 c     &                " kl",kl," itype",itype(il),itype(kl)
 153 c                    write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",
 154 c     &              d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
 155 c                  endif
 156                 endif
 157               enddo  ! l
 158             enddo    ! k
 159           enddo      ! j
 160         enddo        ! i
 161         qq = qq/nl
 162         if (lprn) write (iout,*) "nl",nl," qq",qq
 163       else
 164         write (iout,*)"Error: Q can be computed only for level 1 and 2."
 165       endif
 166       qwolynes=1.0d0-qq
 167       return
 168       end
 169