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       integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
 100      & secseg
 101       integer nsep /3/
 102       double precision dist
 103       double precision dij,d0ij,dijCM,d0ijCM
 104       logical lprn /.false./
 105       logical flag
 106       double precision sigm,x,sim,dd0,fac,ddqij
 107       sigm(x)=0.25d0*x
 108 #ifdef DEBUG
 109       write (iout,*) "qwolynes: flag",flag," seg1 seg1",seg1,seg2,
 110      &   " nsep",nsep
 111       write (iout,*) "nperm",nperm
 112 #endif
 113       do i=0,nres
 114         do j=1,3
 115           dqwol(j,i)=0.0d0
 116           dxqwol(j,i)=0.0d0
 117         enddo
 118       enddo
 119       nl=0
 120       if(flag) then
 121         do il=seg1+nsep,seg2
 122           if (itype(il).eq.ntyp1) cycle
 123           do jl=seg1,il-nsep
 124             if (itype(jl).eq.ntyp1) cycle
 125             nl=nl+1
 126             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 127      &                 (cref(2,jl)-cref(2,il))**2+
 128      &                 (cref(3,jl)-cref(3,il))**2)
 129             dij=dist(il,jl)
 130             sim = 1.0d0/sigm(d0ij)
 131             sim = sim*sim
 132             dd0 = dij-d0ij
 133             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
 134             do k=1,3
 135               ddqij = (c(k,il)-c(k,jl))*fac
 136               dqwol(k,il)=dqwol(k,il)+ddqij
 137               dqwol(k,jl)=dqwol(k,jl)-ddqij
 138             enddo
 139             if (itype(il).ne.10 .or. itype(jl).ne.10) then
 140               nl=nl+1
 141               d0ijCM=dsqrt(
 142      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
 143      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
 144      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
 145               dijCM=dist(il+nres,jl+nres)
 146               sim = 1.0d0/sigm(d0ijCM)
 147               sim = sim*sim
 148               dd0=dijCM-d0ijCM
 149               fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
 150               do k=1,3
 151                 ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
 152                 dxqwol(k,il)=dxqwol(k,il)+ddqij
 153                 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 154               enddo
 155             endif
 156 #ifdef DEBUG
 157             write (iout,*) "prim il",il,itype(il)," jl",jl,itype(jl),
 158      &       " dqwol",(dqwol(k,il),k=1,3)," dxqwol",(dxqwol(k,il),k=1,3)
 159 #endif
 160           enddo
 161         enddo
 162       else
 163         do il=seg1,seg2
 164           if (itype(il).eq.ntyp1) cycle
 165           if((seg3-il).lt.3) then
 166              secseg=il+3
 167           else
 168              secseg=seg3
 169           endif
 170           do jl=secseg,seg4
 171             if (itype(jl).eq.ntyp1) cycle
 172             nl=nl+1
 173             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
 174      &                 (cref(2,jl)-cref(2,il))**2+
 175      &                 (cref(3,jl)-cref(3,il))**2)
 176             dij=dist(il,jl)
 177             sim = 1.0d0/sigm(d0ij)
 178             sim = sim*sim
 179             dd0 = dij-d0ij
 180             fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
 181             do k=1,3
 182               ddqij = (c(k,il)-c(k,jl))*fac
 183               dqwol(k,il)=dqwol(k,il)+ddqij
 184               dqwol(k,jl)=dqwol(k,jl)-ddqij
 185             enddo
 186             if (itype(il).ne.10 .or. itype(jl).ne.10) then
 187               nl=nl+1
 188               d0ijCM=dsqrt(
 189      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
 190      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
 191      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
 192               dijCM=dist(il+nres,jl+nres)
 193               sim = 1.0d0/sigm(d0ijCM)
 194               sim=sim*sim
 195               dd0 = dijCM-d0ijCM
 196               fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
 197               do k=1,3
 198                ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
 199                dxqwol(k,il)=dxqwol(k,il)+ddqij
 200                dxqwol(k,jl)=dxqwol(k,jl)-ddqij
 201               enddo
 202             endif
 203           enddo
 204         enddo
 205       endif
 206 #ifdef DEBUG
 207       write (iout,*) "qwolynes: nl",nl
 208 #endif
 209        do i=0,nres
 210          do j=1,3
 211            dqwol(j,i)=dqwol(j,i)/nl
 212            dxqwol(j,i)=dxqwol(j,i)/nl
 213          enddo
 214        enddo
 215       return
 216       end
 217 c-------------------------------------------------------------------
 218       subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
 219       implicit real*8 (a-h,o-z)
 220       include 'DIMENSIONS'
 221       include 'COMMON.IOUNITS'
 222       include 'COMMON.CHAIN'
 223       include 'COMMON.INTERACT'
 224       include 'COMMON.VAR'
 225       integer seg1,seg2,seg3,seg4
 226       logical flag
 227       double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
 228      & qwolxan(3,0:maxres),q1,q2
 229       double precision delta /1.0d-10/
 230       do i=0,nres
 231         do j=1,3
 232           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 233           cdummy(j,i)=c(j,i)
 234           c(j,i)=c(j,i)+delta
 235           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 236           qwolan(j,i)=(q2-q1)/delta
 237           c(j,i)=cdummy(j,i)
 238         enddo
 239       enddo
 240       do i=0,nres
 241         do j=1,3
 242           q1=qwolynes(seg1,seg2,flag,seg3,seg4)
 243           cdummy(j,i+nres)=c(j,i+nres)
 244           c(j,i+nres)=c(j,i+nres)+delta
 245           q2=qwolynes(seg1,seg2,flag,seg3,seg4)
 246           qwolxan(j,i)=(q2-q1)/delta
 247           c(j,i+nres)=cdummy(j,i+nres)
 248         enddo
 249       enddo
 250 c      write(iout,*) "Numerical Q carteisan gradients backbone: "
 251 c      do i=0,nct
 252 c        write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
 253 c      enddo
 254 c      write(iout,*) "Numerical Q carteisan gradients side-chain: "
 255 c      do i=0,nct
 256 c        write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
 257 c      enddo
 258       return
 259       end
 260 c------------------------------------------------------------------------