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