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         write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres),
 156      &                            ((c(l,k+nres),l=1,3),k=nnt,nct)
 157         call enerprint(energia(0),fT)
 158 #endif
 159         etot=energia(0)
 160         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 161         Fdimless_buf(i)=Fdimless(i)
 162         totfree(i)=etot
 163         totfree_buf(i)=totfree(i)
 164 #ifdef DEBUG
 165         write (iout,*) "fdim calc", i,ii,ib,
 166      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 167      &   entfac(ii),Fdimless(i)
 168 #endif
 169       enddo   ! i
 170
 171       do ijk=1,maxconf
 172       entfac_buf(ijk)=entfac(ijk)
 173       Fdimless_buf(ijk)=Fdimless(ijk)
 174       enddo
 175       do ijk=0,maxconf
 176       totfree_buf(ijk)=totfree(ijk)
 177       enddo
 178
 179
 180 c      scount_buf=scount(me)
 181 c      scount_buf2=scount(0)
 182
 183 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
 184
 185 #ifdef MPI
 186 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
 187       call MPI_Gatherv(Fdimless_buf(1),scount(me),
 188      & MPI_REAL,Fdimless(1),
 189      & scount(0),idispl(0),MPI_REAL,Master,
 190      & MPI_COMM_WORLD, IERROR)
 191 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
 192       call MPI_Gatherv(totfree_buf(1),scount(me),
 193      & MPI_DOUBLE_PRECISION,totfree(1),
 194      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 195      & MPI_COMM_WORLD, IERROR)
 196 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
 197       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
 198      & MPI_DOUBLE_PRECISION,entfac(1),
 199      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 200      & MPI_COMM_WORLD, IERROR)
 201 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 202       if (refstr) then
 203         do i=1,scount(me)
 204           buffer(i)=gdt_ts_tb(i)
 205         enddo
 206         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 207      &   gdt_ts_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 208      &   MPI_COMM_WORLD,IERROR)
 209         do i=1,scount(me)
 210           buffer(i)=gdt_ha_tb(i)
 211         enddo
 212         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 213      &   gdt_ha_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 214      &   MPI_COMM_WORLD,IERROR)
 215         do i=1,scount(me)
 216           buffer(i)=tmscore_tb(i)
 217         enddo
 218         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 219      &   tmscore_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 220      &   MPI_COMM_WORLD,IERROR)
 221       endif
 222       if (me.eq.Master) then
 223 c      WRITE (iout,*) "me.eq.Master"
 224 #endif
 225 #ifdef DEBUG
 226         write (iout,*) "The FDIMLESS array before sorting"
 227         do i=1,ncon
 228           write (iout,*) i,fdimless(i)
 229         enddo
 230 #endif
 231 c      WRITE (iout,*) "Wchodze do call mysort1"
 232         call mysort1(ncon,Fdimless,list_conf)
 233 c      WRITE (iout,*) "Wychodze z call mysort1"
 234 #ifdef DEBUG
 235         write (iout,*) "The FDIMLESS array after sorting"
 236         do i=1,ncon
 237           write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
 238      &     gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
 239         enddo
 240 #endif
 241 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
 242         do i=1,ncon
 243           totfree(i)=fdimless(i)
 244         enddo
 245         qfree=0.0d0
 246         do i=1,ncon
 247           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 248 c          write (iout,*) "fdimless", fdimless(i)
 249         enddo
 250 c        write (iout,*) "qfree",qfree
 251         nlist=1
 252         sumprob=0.0
 253         write (iout,*) "ncon", ncon,maxstr_proc
 254         do i=1,min0(ncon,maxstr_proc)-1
 255           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 256 #ifdef DEBUG
 257           write (iout,*) i,ib,beta_h(ib),
 258      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 259      &     totfree(list_conf(i)),
 260      &     -entfac(list_conf(i)),fdimless(i),sumprob
 261 #endif
 262           if (sumprob.gt.prob_limit) goto 122
 263 c          if (sumprob.gt.1.00d0) goto 122
 264           nlist=nlist+1
 265         enddo
 266   122   continue
 267 #ifdef MPI
 268       endif
 269       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 270      &   IERROR)
 271       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 272      &   IERROR)
 273 c      do iproc=0,nprocs
 274 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 275 c     &   " indend",indend(iproc)
 276 c      enddo
 277       write (iout,*) "nlist",nlist
 278 #endif
 279       return
 280       end
 281 !--------------------------------------------------
 282       subroutine mysort1(n, x, ipermut)
 283       implicit none
 284       integer i,j,imax,ipm,n
 285       real x(n)
 286       integer ipermut(n)
 287       real xtemp
 288       do i=1,n
 289         xtemp=x(i)
 290         imax=i
 291         do j=i+1,n
 292           if (x(j).lt.xtemp) then
 293             imax=j
 294             xtemp=x(j)
 295           endif
 296         enddo
 297         x(imax)=x(i)
 298         x(i)=xtemp
 299         ipm=ipermut(imax)
 300         ipermut(imax)=ipermut(i)
 301         ipermut(i)=ipm
 302       enddo
 303       return
 304       end