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