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