source/unres/src-HCD-5D/q_measure.F

   1       double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.VAR'
   8       integer i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,
   9      & secseg
  10       integer nsep /3/
  11       double precision dist,qm
  12       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  13       logical lprn /.false./
  14       logical flag
  15       double precision sigm,x
  16       sigm(x)=0.25d0*x
  17 #ifdef DEBUG
  18       write (iout,*) "qwolynes: nperm",nperm," flag",flag,
  19      &  " seg1",seg1," seg2",seg2," nsep",nsep
  20 #endif
  21       qq = 0.0d0
  22       nl=0
  23       if(flag) then
  24         do il=seg1+nsep,seg2
  25           if (itype(il).eq.ntyp1) cycle
  26           do jl=seg1,il-nsep
  27             if (itype(jl).eq.ntyp1) cycle
  28             nl=nl+1
  29             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
  30      &                 (cref(2,jl)-cref(2,il))**2+
  31      &                 (cref(3,jl)-cref(3,il))**2)
  32             dij=dist(il,jl)
  33             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  34             qq = qq+qqij
  35             if (itype(il).ne.10 .or. itype(jl).ne.10) then
  36               nl=nl+1
  37               d0ijCM=dsqrt(
  38      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
  39      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
  40      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
  41               dijCM=dist(il+nres,jl+nres)
  42               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  43               qq = qq+qqijCM
  44             endif
  45 c            write (iout,*) "il",il,itype(il)," jl",jl,itype(jl),
  46 c     &        " qqiij",qqij," qqijCM",qqijCM
  47           enddo
  48         enddo
  49 #ifdef DEBUG
  50         write (iout,*) "qwolynes: nl",nl
  51 #endif
  52         qq = qq/nl
  53       else
  54         do il=seg1,seg2
  55           if (itype(il).eq.ntyp1) cycle
  56           if((seg3-il).lt.3) then
  57              secseg=il+3
  58           else
  59              secseg=seg3
  60           endif
  61           do jl=secseg,seg4
  62             if (itype(jl).eq.ntyp1) cycle
  63             nl=nl+1
  64             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
  65      &                 (cref(2,jl)-cref(2,il))**2+
  66      &                 (cref(3,jl)-cref(3,il))**2)
  67             dij=dist(il,jl)
  68             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  69             qq = qq+qqij
  70             if (itype(il).ne.10 .or. itype(jl).ne.10) then
  71               nl=nl+1
  72               d0ijCM=dsqrt(
  73      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
  74      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
  75      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
  76               dijCM=dist(il+nres,jl+nres)
  77               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  78             endif
  79 c            write (iout,*) "il",il,itype(il)," jl",jl,itype(jl),
  80 c     &        " qqiij",qqij," qqijCM",qqijCM
  81             qq = qq+qqijCM
  82           enddo
  83         enddo
  84       qq = qq/nl
  85       endif
  86 c      write (iout,*) "qq",qq
  87       qwolynes=1.0d0-qq
  88       return
  89       end
  90 c-------------------------------------------------------------------
  91       subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
  92       implicit real*8 (a-h,o-z)
  93       include 'DIMENSIONS'
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.CHAIN'
  96       include 'COMMON.INTERACT'
  97       include 'COMMON.VAR'
  98       include 'COMMON.MD'
  99       include 'COMMON.QRESTR'
 100       integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
 101      & secseg
 102       integer nsep /3/
 103       double precision dist
 104       double precision dij,d0ij,dijCM,d0ijCM
 105       logical lprn /.false./
 106       logical flag
 107       double precision sigm,x,sim,dd0,fac,ddqij
 108       sigm(x)=0.25d0*x
 109 #ifdef DEBUG
 110       write (iout,*) "qwolynes: flag",flag," seg1 seg1",seg1,seg2,
 111      &   " nsep",nsep
 112       write (iout,*) "nperm",nperm
 113 #endif
 114       do i=0,nres
 115         do j=1,3
 116           dqwol(j,i)=0.0d0
 117           dxqwol(j,i)=0.0d0
 118         enddo
 119       enddo
 120       nl=0
 121       if(flag) then
 122         do il=seg1+nsep,seg2
 123           if (itype(il).eq.ntyp1) cycle
 124           do jl=seg1,il-nsep
 125             if (itype(jl).eq.ntyp1) cycle
 126             nl=nl+1
 127             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 128      &                 (cref(2,jl)-cref(2,il))**2+
 129      &                 (cref(3,jl)-cref(3,il))**2)
 130             dij=dist(il,jl)
 131             sim = 1.0d0/sigm(d0ij)
 132             sim = sim*sim
 133             dd0 = dij-d0ij
 134             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
 135             do k=1,3
 136               ddqij = (c(k,il)-c(k,jl))*fac
 137               dqwol(k,il)=dqwol(k,il)+ddqij
 138               dqwol(k,jl)=dqwol(k,jl)-ddqij
 139             enddo
 140             if (itype(il).ne.10 .or. itype(jl).ne.10) then
 141               nl=nl+1
 142               d0ijCM=dsqrt(
 143      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
 144      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
 145      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
 146               dijCM=dist(il+nres,jl+nres)
 147               sim = 1.0d0/sigm(d0ijCM)
 148               sim = sim*sim
 149               dd0=dijCM-d0ijCM
 150               fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
 151               do k=1,3
 152                 ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
 153                 dxqwol(k,il)=dxqwol(k,il)+ddqij
 154                 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 155               enddo
 156             endif
 157 #ifdef DEBUG
 158             write (iout,*) "prim il",il,itype(il)," jl",jl,itype(jl),
 159      &       " dqwol",(dqwol(k,il),k=1,3)," dxqwol",(dxqwol(k,il),k=1,3)
 160 #endif
 161           enddo
 162         enddo
 163       else
 164         do il=seg1,seg2
 165           if (itype(il).eq.ntyp1) cycle
 166           if((seg3-il).lt.3) then
 167              secseg=il+3
 168           else
 169              secseg=seg3
 170           endif
 171           do jl=secseg,seg4
 172             if (itype(jl).eq.ntyp1) cycle
 173             nl=nl+1
 174             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 175      &                 (cref(2,jl)-cref(2,il))**2+
 176      &                 (cref(3,jl)-cref(3,il))**2)
 177             dij=dist(il,jl)
 178             sim = 1.0d0/sigm(d0ij)
 179             sim = sim*sim
 180             dd0 = dij-d0ij
 181             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
 182             do k=1,3
 183               ddqij = (c(k,il)-c(k,jl))*fac
 184               dqwol(k,il)=dqwol(k,il)+ddqij
 185               dqwol(k,jl)=dqwol(k,jl)-ddqij
 186             enddo
 187             if (itype(il).ne.10 .or. itype(jl).ne.10) then
 188               nl=nl+1
 189               d0ijCM=dsqrt(
 190      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
 191      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
 192      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
 193               dijCM=dist(il+nres,jl+nres)
 194               sim = 1.0d0/sigm(d0ijCM)
 195               sim=sim*sim
 196               dd0 = dijCM-d0ijCM
 197               fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
 198               do k=1,3
 199                ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
 200                dxqwol(k,il)=dxqwol(k,il)+ddqij
 201                dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 202               enddo
 203             endif
 204           enddo
 205         enddo
 206       endif
 207 #ifdef DEBUG
 208       write (iout,*) "qwolynes: nl",nl
 209 #endif
 210        do i=0,nres
 211          do j=1,3
 212            dqwol(j,i)=dqwol(j,i)/nl
 213            dxqwol(j,i)=dxqwol(j,i)/nl
 214          enddo
 215        enddo
 216       return
 217       end
 218 c-------------------------------------------------------------------
 219       subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
 220       implicit real*8 (a-h,o-z)
 221       include 'DIMENSIONS'
 222       include 'COMMON.IOUNITS'
 223       include 'COMMON.CHAIN'
 224       include 'COMMON.INTERACT'
 225       include 'COMMON.VAR'
 226       integer seg1,seg2,seg3,seg4
 227       logical flag
 228       double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
 229      & qwolxan(3,0:maxres),q1,q2
 230       double precision delta /1.0d-10/
 231       do i=0,nres
 232         do j=1,3
 233           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 234           cdummy(j,i)=c(j,i)
 235           c(j,i)=c(j,i)+delta
 236           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 237           qwolan(j,i)=(q2-q1)/delta
 238           c(j,i)=cdummy(j,i)
 239         enddo
 240       enddo
 241       do i=0,nres
 242         do j=1,3
 243           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 244           cdummy(j,i+nres)=c(j,i+nres)
 245           c(j,i+nres)=c(j,i+nres)+delta
 246           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 247           qwolxan(j,i)=(q2-q1)/delta
 248           c(j,i+nres)=cdummy(j,i+nres)
 249         enddo
 250       enddo
 251 c      write(iout,*) "Numerical Q carteisan gradients backbone: "
 252 c      do i=0,nct
 253 c        write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
 254 c      enddo
 255 c      write(iout,*) "Numerical Q carteisan gradients side-chain: "
 256 c      do i=0,nct
 257 c        write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
 258 c      enddo
 259       return
 260       end
 261 c------------------------------------------------------------------------