source/wham/src-M/qwolynes.f

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