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