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