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

   1       double precision function qwolynes(ib,iref,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 qcontrib
  17       external qcontrib
  18       logical lprn /.false./
  19       logical flag
  20       qq = 0.0d0
  21       nl=0
  22       if (lprn) then
  23       call flush(iout)
  24       endif
  25       do il=nnt+nsep,nct
  26         if (itype(il).eq.ntyp1) cycle
  27         do jl=nnt,il-nsep
  28           if (itype(jl).eq.ntyp1) cycle
  29           nl=nl+1
  30           if (itype(il).ne.10) then
  31             ilnres=il+nres
  32           else
  33             ilnres=il
  34           endif
  35           if (itype(jl).ne.10) then
  36             jlnres=jl+nres
  37           else
  38             jlnres=jl
  39           endif
  40           qqijCM = qcontrib(il,jl,ilnres,jlnres,iref,ib,iprot)
  41           qq = qq+qqijCM
  42           if (lprn) then
  43             write (iout,*) "qqijCM",qqijCM
  44             call flush(iout)
  45           endif
  46         enddo
  47       enddo
  48       if (lprn) then
  49         write (iout,*) "nl",nl," qq",qq
  50         call flush(iout)
  51       endif
  52       qwolynes = qq/nl
  53       return
  54       end
  55 c---------------------------------------------------------------------------
  56       double precision function qcontrib(il,jl,il1,jl1,iref,ib,iprot)
  57       implicit none
  58       include 'DIMENSIONS'
  59       include 'DIMENSIONS.ZSCOPT'
  60       include 'COMMON.IOUNITS'
  61       include 'COMMON.COMPAR'
  62       include 'COMMON.CHAIN'
  63       include 'COMMON.INTERACT'
  64       include 'COMMON.VAR'
  65       include 'COMMON.HEADER'
  66       include 'COMMON.CLASSES'
  67       integer i,j,k,il,jl,il1,jl1,nd,iref,ib,iprot
  68       double precision dist
  69       external dist
  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      &           (cref(1,jl,iref,ib,iprot)-cref(1,il,iref,ib,iprot))**2+
  75      &           (cref(2,jl,iref,ib,iprot)-cref(2,il,iref,ib,iprot))**2+
  76      &           (cref(3,jl,iref,ib,iprot)-cref(3,il,iref,ib,iprot))**2)
  77       dij1=dist(il,jl)
  78       ddave=(dij1-d0ij1)**2
  79       nd=1
  80       if (jl1.ne.jl) then
  81         d0ij2=dsqrt((cref(1,jl1,iref,ib,iprot)-
  82      &          cref(1,il,iref,ib,iprot))**2+
  83      &         (cref(2,jl1,iref,ib,iprot)-cref(2,il,iref,ib,iprot))**2+
  84      &         (cref(3,jl1,iref,ib,iprot)-cref(3,il,iref,ib,iprot))**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((cref(1,jl,iref,ib,iprot)-
  91      &          cref(1,il1,iref,ib,iprot))**2+
  92      &         (cref(2,jl,iref,ib,iprot)-cref(2,il1,iref,ib,iprot))**2+
  93      &         (cref(3,jl,iref,ib,iprot)-cref(3,il1,iref,ib,iprot))**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((cref(1,jl1,iref,ib,iprot)-
 100      &          cref(1,il1,iref,ib,iprot))**2+
 101      &         (cref(2,jl1,iref,ib,iprot)-cref(2,il1,iref,ib,iprot))**2+
 102      &         (cref(3,jl1,iref,ib,iprot)-cref(3,il1,iref,ib,iprot))**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       qcontrib = dexp(-0.5d0*ddave*ssij)
 125       return
 126       end