source/cluster/wham/src-M/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
  26       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
  27       double precision qfree,sumprob,eini,efree,rmsdev
  28       character*80 bxname
  29       character*2 licz1
  30       character*5 ctemper
  31       integer ilen,ijk
  32       external ilen
  33       character*80 structure/'Structure'/
  34       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
  35       double precision energia(0:max_ene), totfree_buf(0:maxconf),
  36      &  entfac_buf(maxconf)
  37       double precision buffer(maxconf)
  38       do i=1,ncon
  39         list_conf(i)=i
  40       enddo
  41 c      do i=1,ncon
  42 c        write (iout,*) i,list_conf(i)
  43 c      enddo
  44 #ifdef MPI
  45       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  46       call daread_ccoords(indstart(me),indend(me))
  47 #endif
  48 C      write (iout,*) "ncon",ncon
  49 C      call flush(iout)
  50       temper=1.0d0/(beta_h(ib)*1.987D-3)
  51 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
  52 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  53 c      quotl=1.0d0
  54 c      kfacl=1.0d0
  55 c      do l=1,5
  56 c        quotl1=quotl
  57 c        quotl=quotl*quot
  58 c        kfacl=kfacl*kfac
  59 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
  60 c      enddo
  61 C#define DEBUG
  62             if (rescale_mode.eq.1) then
  63               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  64               quotl=1.0d0
  65               kfacl=1.0d0
  66               do l=1,5
  67                 quotl1=quotl
  68                 quotl=quotl*quot
  69                 kfacl=kfacl*kfac
  70                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
  71               enddo
  72 #if defined(FUNCTH)
  73               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  74      &                  320.0d0
  75               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  76              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  77      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  78 #elif defined(FUNCT)
  79               fT(6)=betaT/T0
  80               ftprim(6)=1.0d0/T0
  81               ftbis(6)=0.0d0
  82 #else
  83               fT(6)=1.0d0
  84               ftprim(6)=0.0d0
  85               ftbis(6)=0.0d0
  86 #endif
  87
  88             else if (rescale_mode.eq.2) then
  89               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  90               quotl=1.0d0
  91               do l=1,5
  92                 quotl=quotl*quot
  93                 fT(l)=1.12692801104297249644d0/
  94      &             dlog(dexp(quotl)+dexp(-quotl))
  95               enddo
  96 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
  97 c              call flush(iout)
  98 #if defined(FUNCTH)
  99               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 100      &                  320.0d0
 101               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 102              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 103      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 104 #elif defined(FUNCT)
 105               fT(6)=betaT/T0
 106               ftprim(6)=1.0d0/T0
 107               ftbis(6)=0.0d0
 108 #else
 109               fT(6)=1.0d0
 110               ftprim(6)=0.0d0
 111               ftbis(6)=0.0d0
 112 #endif
 113             endif
 114
 115 #ifdef MPI
 116       do i=1,scount(me)
 117         ii=i+indstart(me)-1
 118 #else
 119       do i=1,ncon
 120         ii=i
 121 #endif
 122 C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
 123 c        call flush(iout)
 124 c        if (ib.eq.1) then
 125           do j=1,nres
 126             do k=1,3
 127               c(k,j)=allcart(k,j,i)
 128               c(k,j+nres)=allcart(k,j+nres,i)
 129 C              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
 130 C     &        c(k,j+nres),allcart(k,j+nres,i)
 131             enddo
 132           enddo
 133 C          write(iout,*) "out of j loop"
 134 C          call flush(iout)
 135           do k=1,3
 136             c(k,nres+1)=c(k,1)
 137             c(k,nres+nres)=c(k,nres)
 138           enddo
 139 C          write(iout,*) "after nres+nres",nss_all(i)
 140 C          call flush(iout)
 141           nss=nss_all(i)
 142           do j=1,nss
 143             ihpb(j)=ihpb_all(j,i)
 144             jhpb(j)=jhpb_all(j,i)
 145           enddo
 146           call int_from_cart1(.false.)
 147           call etotal(energia(0),fT)
 148           if (refstr) then
 149 c            write (structure(9:),'(bz,i6.6)') i
 150             call TMscore_sub(rmsdev,gdt_ts_tb(i),
 151      &      gdt_ha_tb(i),tmscore_tb(i),Structure,.false.)
 152 #ifdef DEBUG
 153             write (iout,*) i,rmsdev,gdt_ts_tb(i),gdt_ha_tb(i),
 154      &        tmscore_tb(i)
 155 #endif
 156           endif
 157           totfree(i)=energia(0)
 158           totfree_buf(i)=totfree(i)
 159 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 160 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 161 c          call pdbout(totfree(i),16,i)
 162 c          call flush(iout)
 163 c#define DEBUG
 164 #ifdef DEBUG
 165           write (iout,*) "conformation", i
 166           call enerprint(energia(0),fT)
 167 #endif
 168 c#undef DEBUG
 169           do k=1,max_ene
 170             enetb(k,i)=energia(k)
 171           enddo
 172 c        endif
 173         evdw=enetb(1,i)
 174 c        write (iout,*) evdw
 175        etot=energia(0)
 176 #ifdef SCP14
 177         evdw2_14=enetb(17,i)
 178         evdw2=enetb(2,i)+evdw2_14
 179 #else
 180         evdw2=enetb(2,i)
 181         evdw2_14=0.0d0
 182 #endif
 183 #ifdef SPLITELE
 184         ees=enetb(3,i)
 185         evdw1=enetb(16,i)
 186 #else
 187         ees=enetb(3,i)
 188         evdw1=0.0d0
 189 #endif
 190         ecorr=enetb(4,i)
 191         ecorr5=enetb(5,i)
 192         ecorr6=enetb(6,i)
 193         eel_loc=enetb(7,i)
 194         eello_turn3=enetb(8,i)
 195         eello_turn4=enetb(9,i)
 196         eturn6=enetb(10,i)
 197         ebe=enetb(11,i)
 198         escloc=enetb(12,i)
 199         etors=enetb(13,i)
 200         etors_d=enetb(14,i)
 201         ehpb=enetb(15,i)
 202         estr=enetb(18,i)
 203         esccor=enetb(19,i)
 204         edihcnstr=enetb(20,i)
 205 c#ifdef SPLITELE
 206 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
 207 c     &ft(1)*welec*ees+wvdwpp*evdw1
 208 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 209 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 210 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 211 c     &  +ft(2)*wturn3*eello_turn3
 212 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 213 c     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 214 c     &  +wbond*estr
 215 c#else
 216 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
 217 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 218 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 219 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 220 c     &  +ft(2)*wturn3*eello_turn3
 221 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 222 c     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 223 c     &  +wbond*estr
 224 c#endif
 225 #ifdef DEBUG
 226         write (iout,*) "etot2", etot
 227         write (iout,*) "evdw", wsc, evdw,evdw_t
 228         write (iout,*) "evdw2", wscp, evdw2
 229         write (iout,*) "welec", ft(1),welec,ees
 230         write (iout,*) "evdw1", wvdwpp,evdw1
 231         write (iout,*) "ebe", ebe,wang
 232 #endif
 233         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 234         Fdimless_buf(i)=Fdimless(i)
 235         totfree(i)=etot
 236         totfree_buf(i)=totfree(i)
 237 #ifdef DEBUG
 238         write (iout,*) "fdim calc", i,ii,ib,
 239      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 240      &   entfac(ii),Fdimless(i)
 241 #endif
 242       enddo   ! i
 243
 244       do ijk=1,maxconf
 245       entfac_buf(ijk)=entfac(ijk)
 246       Fdimless_buf(ijk)=Fdimless(ijk)
 247       enddo
 248       do ijk=0,maxconf
 249       totfree_buf(ijk)=totfree(ijk)
 250       enddo
 251
 252
 253 c      scount_buf=scount(me)
 254 c      scount_buf2=scount(0)
 255
 256 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
 257
 258 #ifdef MPI
 259 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
 260       call MPI_Gatherv(Fdimless_buf(1),scount(me),
 261      & MPI_REAL,Fdimless(1),
 262      & scount(0),idispl(0),MPI_REAL,Master,
 263      & MPI_COMM_WORLD, IERROR)
 264 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
 265       call MPI_Gatherv(totfree_buf(1),scount(me),
 266      & MPI_DOUBLE_PRECISION,totfree(1),
 267      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 268      & MPI_COMM_WORLD, IERROR)
 269 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
 270       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
 271      & MPI_DOUBLE_PRECISION,entfac(1),
 272      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 273      & MPI_COMM_WORLD, IERROR)
 274 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 275       if (refstr) then
 276         do i=1,scount(me)
 277           buffer(i)=gdt_ts_tb(i)
 278         enddo
 279         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 280      &   gdt_ts_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 281      &   MPI_COMM_WORLD,IERROR)
 282         do i=1,scount(me)
 283           buffer(i)=gdt_ha_tb(i)
 284         enddo
 285         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 286      &   gdt_ha_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 287      &   MPI_COMM_WORLD,IERROR)
 288         do i=1,scount(me)
 289           buffer(i)=tmscore_tb(i)
 290         enddo
 291         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 292      &   tmscore_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 293      &   MPI_COMM_WORLD,IERROR)
 294       endif
 295       if (me.eq.Master) then
 296 c      WRITE (iout,*) "me.eq.Master"
 297 #endif
 298 #ifdef DEBUG
 299         write (iout,*) "The FDIMLESS array before sorting"
 300         do i=1,ncon
 301           write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
 302      &     gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
 303         enddo
 304 #endif
 305 c      WRITE (iout,*) "Wchodze do call mysort1"
 306         call mysort1(ncon,Fdimless,list_conf)
 307 c      WRITE (iout,*) "Wychodze z call mysort1"
 308 #ifdef DEBUG
 309         write (iout,*) "The FDIMLESS array after sorting"
 310         do i=1,ncon
 311           write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
 312      &     gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
 313         enddo
 314 #endif
 315 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
 316         do i=1,ncon
 317           totfree(i)=fdimless(i)
 318         enddo
 319         qfree=0.0d0
 320         do i=1,ncon
 321           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 322 c          write (iout,*) "fdimless", fdimless(i)
 323         enddo
 324 c        write (iout,*) "qfree",qfree
 325         nlist=1
 326         sumprob=0.0
 327         write (iout,*) "ncon", ncon,maxstr_proc
 328         do i=1,min0(ncon,maxstr_proc)-1
 329           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 330 #ifdef DEBUG
 331           write (iout,*) i,ib,beta_h(ib),
 332      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 333      &     totfree(list_conf(i)),
 334      &     -entfac(list_conf(i)),fdimless(i),sumprob
 335 #endif
 336           if (sumprob.gt.prob_limit) goto 122
 337 c          if (sumprob.gt.1.00d0) goto 122
 338           nlist=nlist+1
 339         enddo
 340   122   continue
 341 #ifdef MPI
 342       endif
 343       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 344      &   IERROR)
 345       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 346      &   IERROR)
 347 c      do iproc=0,nprocs
 348 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 349 c     &   " indend",indend(iproc)
 350 c      enddo
 351       write (iout,*) "nlist",nlist
 352 #endif
 353       return
 354       end
 355 !--------------------------------------------------
 356       subroutine mysort1(n, x, ipermut)
 357       implicit none
 358       integer i,j,imax,ipm,n
 359       real x(n)
 360       integer ipermut(n)
 361       real xtemp
 362       do i=1,n
 363         xtemp=x(i)
 364         imax=i
 365         do j=i+1,n
 366           if (x(j).lt.xtemp) then
 367             imax=j
 368             xtemp=x(j)
 369           endif
 370         enddo
 371         x(imax)=x(i)
 372         x(i)=xtemp
 373         ipm=ipermut(imax)
 374         ipermut(imax)=ipermut(i)
 375         ipermut(i)=ipm
 376       enddo
 377       return
 378       end