source/wham/src-M-SAXS/qwolynes.f

   1       double precision function qwolynes(ilevel,jfrag,kkk)
   2       implicit none
   3       include 'DIMENSIONS'
   4       include 'DIMENSIONS.ZSCOPT'
   5       include 'DIMENSIONS.COMPAR'
   6       include 'COMMON.IOUNITS'
   7       include 'COMMON.COMPAR'
   8       include 'COMMON.CHAIN'
   9       include 'COMMON.INTERACT'
  10       include 'COMMON.CONTROL'
  11       integer ilevel,jfrag,kkk
  12       integer i,j,jl,k,l,il,kl,nl,np,ip,kp
  13       integer nsep /3/
  14       integer iperm
  15       double precision dist
  16       double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
  17       logical lprn /.false./
  18       double precision sigm,x
  19       sigm(x)=0.25d0*x
  20 c      write (iout,*) "QWolyes: " jfrag",jfrag,
  21 c     &  " ilevel",ilevel
  22 c      write (iout,*) "qwolynes: permutation",kkk
  23       qq = 0.0d0
  24       if (ilevel.eq.0) then
  25         if (lprn) write (iout,*) "Q computed for whole molecule"
  26         nl=0
  27         do il=nnt+nsep,nct
  28           if (itype(il).eq.ntyp1) cycle
  29           do jl=nnt,il-nsep
  30             if (itype(jl).eq.ntyp1) cycle
  31             dij=0.0d0
  32             dijCM=0.0d0
  33             d0ij=0.0d0
  34             d0ijCM=0.0d0
  35             qqij=0.0d0
  36             qqijCM=0.0d0
  37             nl=nl+1
  38             d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
  39      &                 (cref(2,jl)-cref(2,il))**2+
  40      &                 (cref(3,jl)-cref(3,il))**2)
  41             dij=dist(iperm(il,kkk),iperm(jl,kkk))
  42             qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  43             if (itype(il).ne.10 .or. itype(jl).ne.10) then
  44               nl=nl+1
  45               d0ijCM=dsqrt(
  46      &               (cref(1,jl+nres)-cref(1,il+nres))**2+
  47      &               (cref(2,jl+nres)-cref(2,il+nres))**2+
  48      &               (cref(3,jl+nres)-cref(3,il+nres))**2)
  49               dijCM=dist(iperm(il,kkk)+nres,iperm(jl,kkk)+nres)
  50               qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  51             endif
  52             qq = qq+qqij+qqijCM
  53             if (lprn) then
  54               write (iout,*) "il",il," jl",jl,
  55      &         " itype",itype(il),itype(jl)
  56               write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
  57      &         " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
  58             endif
  59           enddo
  60         enddo
  61         qq = qq/nl
  62         if (lprn) write (iout,*) "nl",nl," qq",qq
  63       else if (ilevel.eq.1) then
  64         if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag
  65         nl=0
  66 c        write (iout,*) "nlist_frag",nlist_frag(jfrag)
  67         do i=2,nlist_frag(jfrag)
  68           do j=1,i-1
  69             il=list_frag(i,jfrag)
  70             jl=list_frag(j,jfrag)
  71             if (itype(il).eq.ntyp1.or.itype(jl).eq.ntyp1) cycle
  72             if (iabs(il-jl).gt.nsep) then
  73               dij=0.0d0
  74               dijCM=0.0d0
  75               d0ij=0.0d0
  76               d0ijCM=0.0d0
  77               qqij=0.0d0
  78               qqijCM=0.0d0
  79               nl=nl+1
  80               d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
  81      &                   (cref(2,jl)-cref(2,il))**2+
  82      &                   (cref(3,jl)-cref(3,il))**2)
  83               dij=dist(iperm(il,kkk),iperm(jl,kkk))
  84               qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
  85               if (itype(il).ne.10 .or. itype(jl).ne.10) then
  86                 nl=nl+1
  87                 d0ijCM=dsqrt(
  88      &                 (cref(1,jl+nres)-cref(1,il+nres))**2+
  89      &                 (cref(2,jl+nres)-cref(2,il+nres))**2+
  90      &                 (cref(3,jl+nres)-cref(3,il+nres))**2)
  91                 dijCM=dist(iperm(il,kkk)+nres,
  92      &                     iperm(iperm(jl,kkk),kkk)+nres)
  93                 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
  94               endif
  95               qq = qq+qqij+qqijCM
  96               if (lprn) then
  97                 write (iout,*) "i",i," j",j," il",il," jl",jl,
  98      &           " itype",itype(il),itype(jl)
  99                 write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM,
 100      &           " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
 101               endif
 102             endif
 103           enddo
 104         enddo
 105         qq = qq/nl
 106         if (lprn) write (iout,*) "nl",nl," qq",qq
 107       else if (ilevel.eq.2) then
 108         np=npiece(jfrag,ilevel)
 109         nl=0
 110         do i=2,np
 111           ip=ipiece(i,jfrag,ilevel)
 112           do j=1,nlist_frag(ip)
 113             il=list_frag(j,ip)
 114             if (itype(il).eq.ntyp1) cycle
 115             do k=1,i-1
 116               kp=ipiece(k,jfrag,ilevel)
 117               do l=1,nlist_frag(kp)
 118                 kl=list_frag(l,kp)
 119                 if (itype(kl).eq.ntyp1) cycle
 120                 if (iabs(kl-il).gt.nsep) then
 121                   nl=nl+1
 122                   dij=0.0d0
 123                   dijCM=0.0d0
 124                   d0ij=0.0d0
 125                   d0ijCM=0.0d0
 126                   qqij=0.0d0
 127                   qqijCM=0.0d0
 128                   d0ij=dsqrt((cref(1,kl)-cref(1,il))**2+
 129      &                       (cref(2,kl)-cref(2,il))**2+
 130      &                       (cref(3,kl)-cref(3,il))**2)
 131                   dij=dist(iperm(il,kkk),iperm(kl,kkk))
 132                   qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
 133                   if (itype(il).ne.10 .or. itype(kl).ne.10) then
 134                     nl=nl+1
 135                     d0ijCM=dsqrt(
 136      &                 (cref(1,kl+nres)-cref(1,il+nres))**2+
 137      &                 (cref(2,kl+nres)-cref(2,il+nres))**2+
 138      &                 (cref(3,kl+nres)-cref(3,il+nres))**2)
 139                     dijCM=dist(iperm(il,kkk)+nres,iperm(kl,kkk)+nres)
 140                     qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/
 141      &                (sigm(d0ijCM)))**2)
 142                   endif
 143                   qq = qq+qqij+qqijCM
 144                   if (lprn) then
 145                     write (iout,*) "i",i," j",j," k",k," l",l," il",il,
 146      &                " kl",kl," itype",itype(il),itype(kl)
 147                     write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM",
 148      &              d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM
 149                   endif
 150                 endif
 151               enddo  ! l
 152             enddo    ! k
 153           enddo      ! j
 154         enddo        ! i
 155         qq = qq/nl
 156         if (lprn) write (iout,*) "nl",nl," qq",qq
 157       else
 158         write (iout,*)"Error: Q can be computed only for level 1 and 2."
 159       endif
 160       qwolynes=1.0d0-qq
 161       return
 162       end
 163 c-------------------------------------------------------------------------------
 164       subroutine fragment_list
 165       implicit none
 166       include 'DIMENSIONS'
 167       include 'DIMENSIONS.ZSCOPT'
 168       include 'DIMENSIONS.COMPAR'
 169       include 'COMMON.IOUNITS'
 170       include 'COMMON.COMPAR'
 171       logical lprn /.true./
 172       integer i,ilevel,j,k,jfrag
 173       do jfrag=1,nfrag(1)
 174         nlist_frag(jfrag)=0
 175         do i=1,npiece(jfrag,1)
 176           if (lprn) write (iout,*) "jfrag=",jfrag,
 177      &      "i=",i," fragment",ifrag(1,i,jfrag),
 178      &      ifrag(2,i,jfrag)
 179           do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 180             do k=1,nlist_frag(jfrag)
 181               if (list_frag(k,jfrag).eq.j) goto 10
 182             enddo
 183             nlist_frag(jfrag)=nlist_frag(jfrag)+1
 184             list_frag(nlist_frag(jfrag),jfrag)=j
 185           enddo
 186   10      continue
 187         enddo
 188       enddo
 189       write (iout,*) "Fragment list"
 190       do j=1,nfrag(1)
 191         write (iout,*)"Fragment",j," list",(list_frag(k,j),
 192      &   k=1,nlist_frag(j))
 193       enddo
 194       return
 195       end