source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-PDB/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         do jl=nnt,il-nsep
  28           nl=nl+1
  29           if (itype(il).ne.10) then
  30             ilnres=il+nres
  31           else
  32             ilnres=il
  33           endif
  34           if (itype(jl).ne.10) then
  35             jlnres=jl+nres
  36           else
  37             jlnres=jl
  38           endif
  39           qqijCM = qcontribs(il,jl,ilnres,jlnres,creff,iprot)
  40           qq = qq+qqijCM
  41           if (lprn) then
  42             write (iout,*) "qqijCM",qqijCM
  43             call flush(iout)
  44           endif
  45         enddo
  46       enddo
  47       if (lprn) then
  48         write (iout,*) "nl",nl," qq",qq
  49         call flush(iout)
  50       endif
  51       qwolyness = qq/nl
  52       return
  53       end
  54 c---------------------------------------------------------------------------
  55       double precision function qcontribs(il,jl,il1,jl1,creff,iprot)
  56       implicit none
  57       include 'DIMENSIONS'
  58       include 'DIMENSIONS.ZSCOPT'
  59       include 'COMMON.IOUNITS'
  60       include 'COMMON.COMPAR'
  61       include 'COMMON.CHAIN'
  62       include 'COMMON.INTERACT'
  63       include 'COMMON.VAR'
  64       include 'COMMON.HEADER'
  65       include 'COMMON.CLASSES'
  66       integer i,j,k,il,jl,il1,jl1,nd,iref,ib,iprot
  67       double precision dist
  68       external dist
  69       double precision creff(3,maxres2)
  70       double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
  71      &  fac1,ddave,ssij,ddqij
  72       logical lprn /.false./
  73       d0ij1=dsqrt(
  74      &           (creff(1,jl)-creff(1,il))**2+
  75      &           (creff(2,jl)-creff(2,il))**2+
  76      &           (creff(3,jl)-creff(3,il))**2)
  77       dij1=dist(il,jl)
  78       ddave=(dij1-d0ij1)**2
  79       nd=1
  80       if (jl1.ne.jl) then
  81         d0ij2=dsqrt((creff(1,jl1)-
  82      &          creff(1,il))**2+
  83      &         (creff(2,jl1)-creff(2,il))**2+
  84      &         (creff(3,jl1)-creff(3,il))**2)
  85         dij2=dist(il,jl1)
  86         ddave=ddave+(dij2-d0ij2)**2
  87         nd=nd+1
  88       endif
  89       if (il1.ne.il) then
  90         d0ij3=dsqrt((creff(1,jl)-
  91      &          creff(1,il1))**2+
  92      &         (creff(2,jl)-creff(2,il1))**2+
  93      &         (creff(3,jl)-creff(3,il1))**2)
  94         dij3=dist(il1,jl)
  95         ddave=ddave+(dij3-d0ij3)**2
  96         nd=nd+1
  97       endif
  98       if (il1.ne.il .and. jl1.ne.jl) then
  99         d0ij4=dsqrt((creff(1,jl1)-
 100      &          creff(1,il1))**2+
 101      &         (creff(2,jl1)-creff(2,il1))**2+
 102      &         (creff(3,jl1)-creff(3,il1))**2)
 103         dij4=dist(il1,jl1)
 104         ddave=ddave+(dij4-d0ij4)**2
 105         nd=nd+1
 106       endif
 107       ddave=ddave/nd
 108       if (lprn) then
 109         write (iout,*) "il",il," jl",jl,
 110      &  " itype",itype(il),itype(jl)," nd",nd
 111         write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
 112      &  " dij",dij1,dij2,dij3,dij4," ddave",ddave
 113         call flush(iout)
 114       endif
 115 c      ssij = (0.25d0*d0ij1)**2
 116       if (il.ne.il1 .and. jl.ne.jl1) then
 117 c        ssij = 16.0d0/(d0ij1*d0ij4)
 118         ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij4)
 119       else
 120 c        ssij = 16.0d0/(d0ij1*d0ij1)
 121 c        ssij = sigma2(iprot)*sigma2(iprot)/(d0ij1*d0ij1)
 122         ssij = 1.0d0/(sigma2(iprot)*sigma2(iprot))
 123       endif
 124       qcontribs = dexp(-0.5d0*ddave*ssij)
 125       return
 126       end