source/cluster/wham/src-HCD/probabl.F

   1       subroutine probabl(ib,nlist,ncon,*)
   2 ! construct the conformational ensembles at REMD temperatures
   3       implicit none
   4       include "DIMENSIONS"
   5       include "sizesclu.dat"
   6 #ifdef MPI
   7       include "mpif.h"
   8       include "COMMON.MPI"
   9       integer ierror,errcode,status(MPI_STATUS_SIZE)
  10 #endif
  11       include "COMMON.CONTROL"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.FREE"
  14       include "COMMON.FFIELD"
  15       include "COMMON.INTERACT"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.CHAIN"
  18       include "COMMON.CLUSTER"
  19       real*4 csingle(3,maxres2)
  20       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  21      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
  22       double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
  23      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  24      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  25      &      evdw_t,esaxs,eliptran,ethetacnstr,ehomology_constr,
  26      &      edfadis,edfator,edfanei,edfabet
  27       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
  28       double precision qfree,sumprob,eini,efree,rmsdev
  29       character*80 bxname
  30       character*2 licz1
  31       character*5 ctemper
  32       integer ilen,ijk
  33       external ilen
  34       character*80 structure/'Structure'/
  35       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
  36       double precision energia(0:max_ene), totfree_buf(0:maxconf),
  37      &  entfac_buf(maxconf)
  38       double precision buffer(maxconf)
  39       do i=1,ncon
  40         list_conf(i)=i
  41       enddo
  42 c      do i=1,ncon
  43 c        write (iout,*) i,list_conf(i)
  44 c      enddo
  45 #ifdef MPI
  46       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  47       call daread_ccoords(indstart(me),indend(me))
  48 #endif
  49 C      write (iout,*) "ncon",ncon
  50 C      call flush(iout)
  51       temper=1.0d0/(beta_h(ib)*1.987D-3)
  52       if (rescale_mode.eq.1) then
  53         quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  54         quotl=1.0d0
  55         kfacl=1.0d0
  56         do l=1,5
  57           quotl1=quotl
  58           quotl=quotl*quot
  59           kfacl=kfacl*kfac
  60           fT(l)=kfacl/(kfacl-1.0d0+quotl)
  61         enddo
  62 #if defined(FUNCTH)
  63         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  64      &             320.0d0
  65         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  66         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  67      &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  68 #elif defined(FUNCT)
  69         fT(6)=betaT/T0
  70         ftprim(6)=1.0d0/T0
  71         ftbis(6)=0.0d0
  72 #else
  73         fT(6)=1.0d0
  74         ftprim(6)=0.0d0
  75         ftbis(6)=0.0d0
  76 #endif
  77
  78       else if (rescale_mode.eq.2) then
  79         quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  80         quotl=1.0d0
  81         do l=1,5
  82           quotl=quotl*quot
  83           fT(l)=1.12692801104297249644d0/
  84      &       dlog(dexp(quotl)+dexp(-quotl))
  85         enddo
  86 c          write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
  87 c          call flush(iout)
  88 #if defined(FUNCTH)
  89         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  90      &             320.0d0
  91         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  92         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  93      &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  94 #elif defined(FUNCT)
  95         fT(6)=betaT/T0
  96         ftprim(6)=1.0d0/T0
  97         ftbis(6)=0.0d0
  98 #else
  99         fT(6)=1.0d0
 100         ftprim(6)=0.0d0
 101         ftbis(6)=0.0d0
 102 #endif
 103       endif
 104
 105 #ifdef MPI
 106       do i=1,scount(me)
 107         ii=i+indstart(me)-1
 108 #else
 109       do i=1,ncon
 110         ii=i
 111 #endif
 112 C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
 113 c        call flush(iout)
 114 c        if (ib.eq.1) then
 115         do j=1,nres
 116           do k=1,3
 117             c(k,j)=allcart(k,j,i)
 118             c(k,j+nres)=allcart(k,j+nres,i)
 119 C              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
 120 C     &        c(k,j+nres),allcart(k,j+nres,i)
 121           enddo
 122         enddo
 123 C          write(iout,*) "out of j loop"
 124 C          call flush(iout)
 125         do k=1,3
 126           c(k,nres+1)=c(k,1)
 127           c(k,nres+nres)=c(k,nres)
 128         enddo
 129 C          write(iout,*) "after nres+nres",nss_all(i)
 130 C          call flush(iout)
 131         nss=nss_all(i)
 132         do j=1,nss
 133           ihpb(j)=ihpb_all(j,i)
 134           jhpb(j)=jhpb_all(j,i)
 135         enddo
 136         call int_from_cart1(.false.)
 137         call etotal(energia(0),fT)
 138         if (refstr) then
 139           write (structure(9:),'(bz,i6.6)') i
 140           call TMscore_sub(rmsdev,gdt_ts_tb(i),
 141      &    gdt_ha_tb(i),tmscore_tb(i),Structure,.false.)
 142 #ifdef DEBUG
 143           write (iout,*) rmsdev,gdt_ts_tb(i),gdt_ha_tb(i),
 144      &      tmscore_tb(i)
 145 #endif
 146         endif
 147         totfree(i)=energia(0)
 148         totfree_buf(i)=totfree(i)
 149 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 150 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 151 c          call pdbout(totfree(i),16,i)
 152 c          call flush(iout)
 153 #ifdef DEBUG
 154         write (iout,*) "conformation", i
 155         call enerprint(energia(0),fT)
 156 #endif
 157         etot=energia(0)
 158         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 159         Fdimless_buf(i)=Fdimless(i)
 160         totfree(i)=etot
 161         totfree_buf(i)=totfree(i)
 162 #ifdef DEBUG
 163         write (iout,*) "fdim calc", i,ii,ib,
 164      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 165      &   entfac(ii),Fdimless(i)
 166 #endif
 167       enddo   ! i
 168
 169       do ijk=1,maxconf
 170       entfac_buf(ijk)=entfac(ijk)
 171       Fdimless_buf(ijk)=Fdimless(ijk)
 172       enddo
 173       do ijk=0,maxconf
 174       totfree_buf(ijk)=totfree(ijk)
 175       enddo
 176
 177
 178 c      scount_buf=scount(me)
 179 c      scount_buf2=scount(0)
 180
 181 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
 182
 183 #ifdef MPI
 184 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
 185       call MPI_Gatherv(Fdimless_buf(1),scount(me),
 186      & MPI_REAL,Fdimless(1),
 187      & scount(0),idispl(0),MPI_REAL,Master,
 188      & MPI_COMM_WORLD, IERROR)
 189 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
 190       call MPI_Gatherv(totfree_buf(1),scount(me),
 191      & MPI_DOUBLE_PRECISION,totfree(1),
 192      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 193      & MPI_COMM_WORLD, IERROR)
 194 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
 195       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
 196      & MPI_DOUBLE_PRECISION,entfac(1),
 197      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 198      & MPI_COMM_WORLD, IERROR)
 199 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 200       if (refstr) then
 201         do i=1,scount(me)
 202           buffer(i)=gdt_ts_tb(i)
 203         enddo
 204         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 205      &   gdt_ts_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 206      &   MPI_COMM_WORLD,IERROR)
 207         do i=1,scount(me)
 208           buffer(i)=gdt_ha_tb(i)
 209         enddo
 210         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 211      &   gdt_ha_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 212      &   MPI_COMM_WORLD,IERROR)
 213         do i=1,scount(me)
 214           buffer(i)=tmscore_tb(i)
 215         enddo
 216         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 217      &   tmscore_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 218      &   MPI_COMM_WORLD,IERROR)
 219       endif
 220       if (me.eq.Master) then
 221 c      WRITE (iout,*) "me.eq.Master"
 222 #endif
 223 #ifdef DEBUG
 224         write (iout,*) "The FDIMLESS array before sorting"
 225         do i=1,ncon
 226           write (iout,*) i,fdimless(i)
 227         enddo
 228 #endif
 229 c      WRITE (iout,*) "Wchodze do call mysort1"
 230         call mysort1(ncon,Fdimless,list_conf)
 231 c      WRITE (iout,*) "Wychodze z call mysort1"
 232 #ifdef DEBUG
 233         write (iout,*) "The FDIMLESS array after sorting"
 234         do i=1,ncon
 235           write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
 236      &     gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
 237         enddo
 238 #endif
 239 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
 240         do i=1,ncon
 241           totfree(i)=fdimless(i)
 242         enddo
 243         qfree=0.0d0
 244         do i=1,ncon
 245           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 246 c          write (iout,*) "fdimless", fdimless(i)
 247         enddo
 248 c        write (iout,*) "qfree",qfree
 249         nlist=1
 250         sumprob=0.0
 251         write (iout,*) "ncon", ncon,maxstr_proc
 252         do i=1,min0(ncon,maxstr_proc)-1
 253           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 254 #ifdef DEBUG
 255           write (iout,*) i,ib,beta_h(ib),
 256      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 257      &     totfree(list_conf(i)),
 258      &     -entfac(list_conf(i)),fdimless(i),sumprob
 259 #endif
 260           if (sumprob.gt.prob_limit) goto 122
 261 c          if (sumprob.gt.1.00d0) goto 122
 262           nlist=nlist+1
 263         enddo
 264   122   continue
 265 #ifdef MPI
 266       endif
 267       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 268      &   IERROR)
 269       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 270      &   IERROR)
 271 c      do iproc=0,nprocs
 272 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 273 c     &   " indend",indend(iproc)
 274 c      enddo
 275       write (iout,*) "nlist",nlist
 276 #endif
 277       return
 278       end
 279 !--------------------------------------------------
 280       subroutine mysort1(n, x, ipermut)
 281       implicit none
 282       integer i,j,imax,ipm,n
 283       real x(n)
 284       integer ipermut(n)
 285       real xtemp
 286       do i=1,n
 287         xtemp=x(i)
 288         imax=i
 289         do j=i+1,n
 290           if (x(j).lt.xtemp) then
 291             imax=j
 292             xtemp=x(j)
 293           endif
 294         enddo
 295         x(imax)=x(i)
 296         x(i)=xtemp
 297         ipm=ipermut(imax)
 298         ipermut(imax)=ipermut(i)
 299         ipermut(i)=ipm
 300       enddo
 301       return
 302       end