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

   1       subroutine probabl(ibb,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,ehomology_constr
  26       integer i,ii,ik,iproc,iscor,j,k,l,ib,ibb,ib_start,ib_end,nlist,
  27      &  ncon
  28       double precision qfree,sumprob,eini,rmsdev
  29       real*4 probab(maxconf),totprob(maxconf),aux
  30       character*80 bxname
  31       character*2 licz1
  32       character*5 ctemper
  33       integer ilen,ijk
  34       external ilen
  35       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
  36       double precision energia(0:max_ene), totfree_buf(0:maxconf),
  37      &  entfac_buf(maxconf)
  38       if (cumul_prob) then
  39         write (iout,*)
  40      &  "Probabilities will be summed over all temperatures"
  41         ib_start=1
  42         ib_end=nT
  43       else
  44         ib_start=ibb
  45         ib_end=ibb
  46       endif
  47       probab=0.0
  48       totprob=0.0
  49       Fdimless=0.0d0
  50       Fdimless_buf=0.0d0
  51       totfree=0.0d0
  52       totfree_buf=0.0d0
  53       do i=1,ncon
  54         list_conf(i)=i
  55       enddo
  56       ib=1
  57       if (rescale_mode.eq.1) then
  58         quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  59         quotl=1.0d0
  60         kfacl=1.0d0
  61         do l=1,5
  62           quotl1=quotl
  63           quotl=quotl*quot
  64           kfacl=kfacl*kfac
  65           fT(l)=kfacl/(kfacl-1.0d0+quotl)
  66         enddo
  67 #if defined(FUNCTH)
  68         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  69      &                320.0d0
  70         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  71         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  72      &            /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  73 #elif defined(FUNCT)
  74         fT(6)=betaT/T0
  75         ftprim(6)=1.0d0/T0
  76         ftbis(6)=0.0d0
  77 #else
  78         fT(6)=1.0d0
  79         ftprim(6)=0.0d0
  80         ftbis(6)=0.0d0
  81 #endif
  82
  83       else if (rescale_mode.eq.2) then
  84         quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  85         quotl=1.0d0
  86         do l=1,5
  87           quotl=quotl*quot
  88           fT(l)=1.12692801104297249644d0/
  89      &           dlog(dexp(quotl)+dexp(-quotl))
  90         enddo
  91 c            write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
  92 c            call flush(iout)
  93 #if defined(FUNCTH)
  94             ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  95      &                320.0d0
  96             ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  97             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  98      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  99 #elif defined(FUNCT)
 100             fT(6)=betaT/T0
 101             ftprim(6)=1.0d0/T0
 102             ftbis(6)=0.0d0
 103 #else
 104             fT(6)=1.0d0
 105             ftprim(6)=0.0d0
 106             ftbis(6)=0.0d0
 107 #endif
 108           endif
 109 c      do i=1,ncon
 110 c        write (iout,*) i,list_conf(i)
 111 c      enddo
 112 #ifdef MPI
 113       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
 114       call daread_ccoords(indstart(me),indend(me))
 115 #endif
 116 C      write (iout,*) "ncon",ncon
 117 C      call flush(iout)
 118
 119 #ifdef MPI
 120       do i=1,scount(me)
 121         ii=i+indstart(me)-1
 122 #else
 123       do i=1,ncon
 124         ii=i
 125 #endif
 126         do j=1,nres
 127           do k=1,3
 128             c(k,j)=allcart(k,j,i)
 129             c(k,j+nres)=allcart(k,j+nres,i)
 130 C            write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
 131 C     &      c(k,j+nres),allcart(k,j+nres,i)
 132           enddo
 133         enddo
 134 C        write(iout,*) "out of j loop"
 135 C        call flush(iout)
 136         do k=1,3
 137           c(k,nres+1)=c(k,1)
 138           c(k,nres+nres)=c(k,nres)
 139         enddo
 140 C        write(iout,*) "after nres+nres",nss_all(i)
 141 C        call flush(iout)
 142         nss=nss_all(i)
 143         do j=1,nss
 144           ihpb(j)=ihpb_all(j,i)
 145           jhpb(j)=jhpb_all(j,i)
 146         enddo
 147         call int_from_cart1(.false.)
 148 C        write(iout,*) "before etotal"
 149 C        call flush(iout)
 150         call etotal(energia(0),fT)
 151 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 152 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 153 c          call enerprint(energia(0),fT)
 154 c          call pdbout(totfree(i),16,i)
 155 #ifdef DEBUG
 156         write (iout,*) i," energia",(energia(j),j=0,max_ene)
 157         write (iout,*) "etot", etot
 158         write (iout,*) "ft(6)", ft(6)
 159 #endif
 160         do k=1,max_ene
 161           enetb(k,i)=energia(k)
 162         enddo
 163       enddo ! i
 164 c#endif
 165 C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
 166 c        call flush(iout)
 167       do ib=ib_start,ib_end
 168           temper=1.0d0/(beta_h(ib)*1.987D-3)
 169 #ifdef MPI
 170         do i=1,scount(me)
 171           ii=i+indstart(me)-1
 172 #else
 173         do i=1,ncon
 174           ii=i
 175 #endif
 176           evdw=enetb(1,i)
 177 c          write (iout,*) evdw
 178           etot=energia(0)
 179 #ifdef SCP14
 180           evdw2_14=enetb(17,i)
 181           evdw2=enetb(2,i)+evdw2_14
 182 #else
 183           evdw2=enetb(2,i)
 184           evdw2_14=0.0d0
 185 #endif
 186 #ifdef SPLITELE
 187           ees=enetb(3,i)
 188           evdw1=enetb(16,i)
 189 #else
 190           ees=enetb(3,i)
 191           evdw1=0.0d0
 192 #endif
 193           ecorr=enetb(4,i)
 194           ecorr5=enetb(5,i)
 195           ecorr6=enetb(6,i)
 196           eel_loc=enetb(7,i)
 197           eello_turn3=enetb(8,i)
 198           eello_turn4=enetb(9,i)
 199           eturn6=enetb(10,i)
 200           ebe=enetb(11,i)
 201           escloc=enetb(12,i)
 202           etors=enetb(13,i)
 203           etors_d=enetb(14,i)
 204           ehpb=enetb(15,i)
 205           estr=enetb(18,i)
 206           esccor=enetb(19,i)
 207           edihcnstr=enetb(20,i)
 208           evdw_t=enetb(21,i)
 209           eliptran=enetb(22,i)
 210           ehomology_constr=enetb(24,i)
 211           esaxs=enetb(25,i)
 212           if (rescale_mode.eq.1) then
 213             quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
 214             quotl=1.0d0
 215             kfacl=1.0d0
 216             do l=1,5
 217               quotl1=quotl
 218               quotl=quotl*quot
 219               kfacl=kfacl*kfac
 220               fT(l)=kfacl/(kfacl-1.0d0+quotl)
 221             enddo
 222 #if defined(FUNCTH)
 223             ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 224      &                320.0d0
 225             ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 226             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 227      &            /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 228 #elif defined(FUNCT)
 229             fT(6)=betaT/T0
 230             ftprim(6)=1.0d0/T0
 231             ftbis(6)=0.0d0
 232 #else
 233             fT(6)=1.0d0
 234             ftprim(6)=0.0d0
 235             ftbis(6)=0.0d0
 236 #endif
 237
 238           else if (rescale_mode.eq.2) then
 239             quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
 240             quotl=1.0d0
 241             do l=1,5
 242               quotl=quotl*quot
 243               fT(l)=1.12692801104297249644d0/
 244      &           dlog(dexp(quotl)+dexp(-quotl))
 245             enddo
 246 c            write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
 247 c            call flush(iout)
 248 #if defined(FUNCTH)
 249             ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 250      &                320.0d0
 251             ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 252             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 253      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 254 #elif defined(FUNCT)
 255             fT(6)=betaT/T0
 256             ftprim(6)=1.0d0/T0
 257             ftbis(6)=0.0d0
 258 #else
 259             fT(6)=1.0d0
 260             ftprim(6)=0.0d0
 261             ftbis(6)=0.0d0
 262 #endif
 263           endif
 264 c          write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
 265 c          call flush(iout)
 266 c          call etotal(energia(0),fT)
 267 c#ifdef CHUJ
 268 #ifdef SPLITELE
 269           if (shield_mode.gt.0) then
 270           etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 271      &     +welec*ft(1)*ees
 272      &     +ft(1)*wvdwpp*evdw1
 273      &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 274      &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 275      &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 276      &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 277      &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 278      &     +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr
 279      &     +wliptran*eliptran+wsaxs*esaxs
 280           else
 281           etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+welec*ft(1)*ees
 282      &     +wvdwpp*evdw1
 283      &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 284      &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 285      &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 286      &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 287      &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 288      &     +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr
 289      &     +wliptran*eliptran+wsaxs*esaxs
 290           endif
 291 #else
 292           if (shield_mode.gt.0) then
 293           etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 294      &     +welec*ft(1)*(ees+evdw1)
 295      &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 296      &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 297      &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 298      &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 299      &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 300      &     +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr
 301      &     +wliptran*eliptran+wsaxs*esaxs
 302           else
 303           etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 304      &     +welec*ft(1)*(ees+evdw1)
 305      &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
 306      &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
 307      &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
 308      &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
 309      &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
 310      &     +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr
 311      &     +wliptran*eliptran+wsaxs*esaxs
 312           endif
 313 #endif
 314 #ifdef DEBUG
 315           write (iout,*) "etot2", etot
 316           write (iout,*) "evdw", wsc, evdw,evdw_t
 317           write (iout,*) "evdw2", wscp, evdw2
 318           write (iout,*) "welec", ft(1),welec,ees
 319           write (iout,*) "evdw1", wvdwpp,evdw1
 320           write (iout,*) "ebe", ebe,wang
 321 #endif
 322 c#else
 323 c          etot=energia(0)
 324 c#endif
 325 #ifdef DEBUG
 326           write (iout,'(2i5,30f10.2)') i,ii,etot,
 327      &    (enetb(ijk,i),ijk=1,max_ene)
 328 #endif
 329           Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 330 #ifdef DEBUG
 331           write (iout,*) "fdim calc", i,ii,ib,
 332      &     1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 333      &    entfac(ii),Fdimless(i)
 334 #endif
 335           Fdimless_buf(i)=Fdimless(i)
 336           totfree_buf(i)=totfree(i)
 337         enddo   ! i
 338
 339 #ifdef MPI
 340 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
 341         call MPI_Gatherv(Fdimless_buf(1),scount(me),
 342      &    MPI_REAL,Fdimless(1),
 343      &    scount(0),idispl(0),MPI_REAL,Master,
 344      &    MPI_COMM_WORLD, IERROR)
 345 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
 346         call MPI_Gatherv(totfree_buf(1),scount(me),
 347      &    MPI_DOUBLE_PRECISION,totfree(1),
 348      &    scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 349      &    MPI_COMM_WORLD, IERROR)
 350 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
 351 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 352         aux=Fdimless(1)
 353         do i=1,ncon
 354           if (Fdimless(i).lt.aux) aux=Fdimless(i)
 355         enddo
 356         sumprob=0.0d0
 357         do i=1,ncon
 358           probab(i)=exp(aux-fdimless(i))
 359           sumprob=sumprob+probab(i)
 360         enddo
 361         do i=1,ncon
 362           probab(i)=probab(i)/sumprob
 363           totprob(i)=totprob(i)+probab(i)
 364 c          write (iout,*) ib,i,probab(i),totprob(i)
 365         enddo
 366       enddo ! ib
 367       sumprob=0.0d0
 368       do i=1,ncon
 369         sumprob=sumprob+totprob(i)
 370       enddo
 371       write (iout,*) "sumprob",sumprob
 372       do i=1,ncon
 373         totprob(i)=totprob(i)/sumprob
 374       enddo
 375       do i=1,ncon
 376 c        Fdimless(i)=-dlog(Fdimless(i))
 377         Fdimless(i)=-alog(totprob(i))
 378       enddo
 379       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
 380      &    MPI_DOUBLE_PRECISION,entfac(1),
 381      &    scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 382      &    MPI_COMM_WORLD, IERROR)
 383 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 384       if (me.eq.Master) then
 385 c      WRITE (iout,*) "me.eq.Master"
 386 #endif
 387 #ifdef DEBUG
 388         write (iout,*) "The FDIMLESS array before sorting"
 389         do i=1,ncon
 390 c          write (iout,*) i,fdimless(i)
 391         enddo
 392 #endif
 393 c      WRITE (iout,*) "Wchodze do call mysort1"
 394         call mysort1(ncon,Fdimless,list_conf)
 395 c      WRITE (iout,*) "Wychodze z call mysort1"
 396 #ifdef DEBUG
 397         write (iout,*) "The FDIMLESS array after sorting"
 398         do i=1,ncon
 399           write (iout,*) i,list_conf(i),fdimless(i)
 400         enddo
 401 #endif
 402 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
 403         do i=1,ncon
 404           totfree(i)=fdimless(i)
 405         enddo
 406         qfree=0.0d0
 407         do i=1,ncon
 408           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 409 c          write (iout,*) "fdimless", fdimless(i)
 410         enddo
 411 c        write (iout,*) "qfree",qfree
 412         nlist=1
 413         sumprob=0.0
 414         write (iout,*) "ncon", ncon,maxstr_proc
 415         do i=1,min0(ncon,maxstr_proc)-1
 416           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 417 #ifdef DEBUG
 418           write (iout,*) "tu szukaj ponizej 7"
 419           write (iout,*) i,ib,beta_h(ib),
 420      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 421      &     totfree(list_conf(i)),
 422      &     -entfac(list_conf(i)),fdimless(i),sumprob
 423 #endif
 424           if (sumprob.gt.prob_limit) goto 122
 425 c          if (sumprob.gt.1.00d0) goto 122
 426           nlist=nlist+1
 427         enddo
 428   122   continue
 429 #ifdef MPI
 430       endif
 431       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 432      &   IERROR)
 433       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 434      &   IERROR)
 435 c      do iproc=0,nprocs
 436 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 437 c     &   " indend",indend(iproc)
 438 c      enddo
 439       write (iout,*) "nlist",nlist
 440 #endif
 441       if (cumul_prob) then
 442       write (iout,'(80(1h-)/i4,a,f6.3,a,10f6.1/80(1h-))')
 443      & nlist," conformations found within",prob_limit,
 444      & " probablity cut_off in the temperature range ",
 445      &  (1.0d0/(1.987d-3*beta_h(ib)),ib=1,nT)
 446       else
 447       write (iout,'(80(1h-)/i4,a,f6.3,a,f6.1/80(1h-))')
 448      & nlist," conformations found within",prob_limit,
 449      &" probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ibb))
 450       endif
 451       return
 452       end
 453 !--------------------------------------------------
 454       subroutine mysort1(n, x, ipermut)
 455       implicit none
 456       integer i,j,imax,ipm,n
 457       real x(n)
 458       integer ipermut(n)
 459       real xtemp
 460       do i=1,n
 461         xtemp=x(i)
 462         imax=i
 463         do j=i+1,n
 464           if (x(j).lt.xtemp) then
 465             imax=j
 466             xtemp=x(j)
 467           endif
 468         enddo
 469         x(imax)=x(i)
 470         x(i)=xtemp
 471         ipm=ipermut(imax)
 472         ipermut(imax)=ipermut(i)
 473         ipermut(i)=ipm
 474       enddo
 475       return
 476       end