source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/q_measure_self.F

   1       double precision function qwolyness(creff,iprot)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'COMMON.IOUNITS'
   6       include 'COMMON.COMPAR'
   7       include 'COMMON.CHAIN'
   8       include 'COMMON.INTERACT'
   9       include 'COMMON.VAR'
  10       include 'COMMON.HEADER'
  11       integer ib,iref,iprot
  12       integer i,j,jl,ilnres,jlnres,klnres,k,l,il,kl,nl,np,ip,kp
  13       integer nsep /3/
  14       double precision dist,qm
  15       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  16       double precision qcontribs
  17       external qcontribs
  18       double precision creff(3,maxres2)
  19       logical lprn /.false./
  20       logical flag
  21       qq = 0.0d0
  22       nl=0
  23       if (lprn) then
  24       call flush(iout)
  25       endif
  26       do il=nnt+nsep,nct
  27         if (itype(il).eq.ntyp1) cycle
  28         do jl=nnt,il-nsep
  29           if (itype(jl).eq.ntyp1) cycle
  30           nl=nl+1
  31           if (itype(il).ne.10) then
  32             ilnres=il+nres
  33           else
  34             ilnres=il
  35           endif
  36           if (itype(jl).ne.10) then
  37             jlnres=jl+nres
  38           else
  39             jlnres=jl
  40           endif
  41           qqijCM = qcontribs(il,jl,ilnres,jlnres,creff,iprot)
  42           qq = qq+qqijCM
  43           if (lprn) then
  44             write (iout,*) "qqijCM",qqijCM
  45             call flush(iout)
  46           endif
  47         enddo
  48       enddo
  49       if (lprn) then
  50         write (iout,*) "nl",nl," qq",qq
  51         call flush(iout)
  52       endif
  53       qwolyness = qq/nl
  54       return
  55       end
  56 c---------------------------------------------------------------------------
  57       double precision function qcontribs(il,jl,il1,jl1,creff,iprot)
  58       implicit none
  59       include 'DIMENSIONS'
  60       include 'DIMENSIONS.ZSCOPT'
  61       include 'COMMON.IOUNITS'
  62       include 'COMMON.COMPAR'
  63       include 'COMMON.CHAIN'
  64       include 'COMMON.INTERACT'
  65       include 'COMMON.VAR'
  66       include 'COMMON.HEADER'
  67       include 'COMMON.CLASSES'
  68       integer i,j,k,il,jl,il1,jl1,nd,iref,ib,iprot
  69       double precision dist
  70       external dist
  71       double precision creff(3,maxres2)
  72       double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
  73      &  fac1,ddave,ssij,ddqij
  74       logical lprn /.false./
  75       d0ij1=dsqrt(
  76      &           (creff(1,jl)-creff(1,il))**2+
  77      &           (creff(2,jl)-creff(2,il))**2+
  78      &           (creff(3,jl)-creff(3,il))**2)
  79       dij1=dist(il,jl)
  80       ddave=(dij1-d0ij1)**2
  81       nd=1
  82       if (jl1.ne.jl) then
  83         d0ij2=dsqrt((creff(1,jl1)-
  84      &          creff(1,il))**2+
  85      &         (creff(2,jl1)-creff(2,il))**2+
  86      &         (creff(3,jl1)-creff(3,il))**2)
  87         dij2=dist(il,jl1)
  88         ddave=ddave+(dij2-d0ij2)**2
  89         nd=nd+1
  90       endif
  91       if (il1.ne.il) then
  92         d0ij3=dsqrt((creff(1,jl)-
  93      &          creff(1,il1))**2+
  94      &         (creff(2,jl)-creff(2,il1))**2+
  95      &         (creff(3,jl)-creff(3,il1))**2)
  96         dij3=dist(il1,jl)
  97         ddave=ddave+(dij3-d0ij3)**2
  98         nd=nd+1
  99       endif
 100       if (il1.ne.il .and. jl1.ne.jl) then
 101         d0ij4=dsqrt((creff(1,jl1)-
 102      &          creff(1,il1))**2+
 103      &         (creff(2,jl1)-creff(2,il1))**2+
 104      &         (creff(3,jl1)-creff(3,il1))**2)
 105         dij4=dist(il1,jl1)
 106         ddave=ddave+(dij4-d0ij4)**2
 107         nd=nd+1
 108       endif
 109       ddave=ddave/nd
 110       if (lprn) then
 111         write (iout,*) "il",il," jl",jl,
 112      &  " itype",itype(il),itype(jl)," nd",nd
 113         write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
 114      &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
 115         call flush(iout)
 116       endif
 117 c      ssij = (0.25d0*d0ij1)**2
 118       if (il.ne.il1 .and. jl.ne.jl1) then
 119 c        ssij = 16.0d0/(d0ij1*d0ij4)
 120         ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij4)
 121       else
 122 c        ssij = 16.0d0/(d0ij1*d0ij1)
 123 c        ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij1)
 124         ssij = 1.0d0/(sigma2(iprot)*sigma2(iprot))
 125       endif
 126       qcontribs = dexp(-0.5d0*ddave*ssij)
 127       return
 128       end