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
  31       external ilen
  32       real*4 Fdimless(maxconf),Fdimless_(maxconf)
  33       double precision energia(0:max_ene)
  34       double precision totfree_(maxconf),entfac_(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       temper=1.0d0/(beta_h(ib)*1.987D-3)
  47 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
  48 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  49 c      quotl=1.0d0
  50 c      kfacl=1.0d0
  51 c      do l=1,5
  52 c        quotl1=quotl
  53 c        quotl=quotl*quot
  54 c        kfacl=kfacl*kfac
  55 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
  56 c      enddo
  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
 110 #ifdef MPI
 111       do i=1,scount(me)
 112         ii=i+indstart(me)-1
 113 #else
 114       do i=1,ncon
 115         ii=i
 116 #endif
 117 c        write (iout,*) "i",i," ii",ii
 118 c        call flush(iout)
 119         if (ib.eq.1) then
 120           do j=1,nres
 121             do k=1,3
 122               c(k,j)=allcart(k,j,i)
 123               c(k,j+nres)=allcart(k,j+nres,i)
 124             enddo
 125           enddo
 126           do k=1,3
 127             c(k,nres+1)=c(k,1)
 128             c(k,nres+nres)=c(k,nres)
 129           enddo
 130           nss=nss_all(i)
 131           do j=1,nss
 132             ihpb(j)=ihpb_all(j,i)
 133             jhpb(j)=jhpb_all(j,i)
 134           enddo
 135           call int_from_cart1(.false.)
 136           call etotal(energia(0),fT)
 137           totfree(i)=energia(0)
 138 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 139 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 140 c          call enerprint(energia(0),fT)
 141 c          call pdbout(totfree(i),16,i)
 142 #ifdef DEBUG
 143           write (iout,*) i," energia",(energia(j),j=0,max_ene)
 144           write (iout,*) "etot", etot
 145           write (iout,*) "ft(6)", ft(6)
 146 #endif
 147           do k=1,max_ene
 148             enetb(k,i)=energia(k)
 149           enddo
 150         endif
 151         evdw=enetb(1,i)
 152 c        write (iout,*) evdw
 153        etot=energia(0)
 154 #ifdef SCP14
 155         evdw2_14=enetb(17,i)
 156         evdw2=enetb(2,i)+evdw2_14
 157 #else
 158         evdw2=enetb(2,i)
 159         evdw2_14=0.0d0
 160 #endif
 161 #ifdef SPLITELE
 162         ees=enetb(3,i)
 163         evdw1=enetb(16,i)
 164 #else
 165         ees=enetb(3,i)
 166         evdw1=0.0d0
 167 #endif
 168         ecorr=enetb(4,i)
 169         ecorr5=enetb(5,i)
 170         ecorr6=enetb(6,i)
 171         eel_loc=enetb(7,i)
 172         eello_turn3=enetb(8,i)
 173         eello_turn4=enetb(9,i)
 174         eturn6=enetb(10,i)
 175         ebe=enetb(11,i)
 176         escloc=enetb(12,i)
 177         etors=enetb(13,i)
 178         etors_d=enetb(14,i)
 179         ehpb=enetb(15,i)
 180         estr=enetb(18,i)
 181         esccor=enetb(19,i)
 182         edihcnstr=enetb(20,i)
 183         evdw_t=enetb(21,i)
 184 c#ifdef SPLITELE
 185 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
 186 c     &ft(1)*welec*ees+wvdwpp*evdw1
 187 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 188 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 189 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 190 c     &  +ft(2)*wturn3*eello_turn3
 191 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 192 c     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 193 c     &  +wbond*estr
 194 c#else
 195 c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
 196 c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 197 c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 198 c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 199 c     &  +ft(2)*wturn3*eello_turn3
 200 c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 201 c     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 202 c     &  +wbond*estr
 203 c#endif
 204 #ifdef DEBUG
 205         write (iout,*) "etot2", etot
 206         write (iout,*) "evdw", wsc, evdw,evdw_t
 207         write (iout,*) "evdw2", wscp, evdw2
 208         write (iout,*) "welec", ft(1),welec,ees
 209         write (iout,*) "evdw1", wvdwpp,evdw1
 210         write (iout,*) "ebe" ebe,wang
 211 #endif
 212 #ifdef MPI
 213         Fdimless_(i)=beta_h(ib)*etot+entfac(ii)
 214         totfree_(i)=etot
 215 #else
 216         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 217         totfree(i)=etot
 218 #endif
 219 #ifdef DEBUG
 220         write (iout,*) "fdim calc", i,ii,ib,
 221      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 222      &   entfac(ii),Fdimless(i)
 223 #endif
 224       enddo   ! i
 225 #ifdef MPI
 226       call MPI_Gatherv(Fdimless_(1),scount(me),
 227      & MPI_REAL,Fdimless(1),
 228      & scount(0),idispl(0),MPI_REAL,Master,
 229      & MPI_COMM_WORLD, IERROR)
 230       call MPI_Gatherv(totfree_(1),scount(me),
 231      & MPI_DOUBLE_PRECISION,totfree(1),
 232      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 233      & MPI_COMM_WORLD, IERROR)
 234       call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
 235      & MPI_DOUBLE_PRECISION,entfac_(1),
 236      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 237      & MPI_COMM_WORLD, IERROR)
 238       if (me.eq.Master) then
 239         do i=1,ncon
 240           entfac(i)=entfac_(i)
 241         enddo
 242 #endif
 243 #ifdef DEBUG
 244         write (iout,*) "The FDIMLESS array before sorting"
 245         do i=1,ncon
 246 c          write (iout,*) i,fdimless(i)
 247         enddo
 248 #endif
 249         call mysort1(ncon,Fdimless,list_conf)
 250 #ifdef DEBUG
 251         write (iout,*) "The FDIMLESS array after sorting"
 252         do i=1,ncon
 253           write (iout,*) i,list_conf(i),fdimless(i)
 254         enddo
 255 #endif
 256         do i=1,ncon
 257           totfree(i)=fdimless(i)
 258         enddo
 259         qfree=0.0d0
 260         do i=1,ncon
 261           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 262 c          write (iout,*) "fdimless", fdimless(i)
 263         enddo
 264 c        write (iout,*) "qfree",qfree
 265         nlist=1
 266         sumprob=0.0
 267         write (iout,*) "ncon", ncon,maxstr_proc
 268         do i=1,min0(ncon,maxstr_proc)-1
 269           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 270 #ifdef DEBUG
 271           write (iout,*) i,ib,beta_h(ib),
 272      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 273      &     totfree(list_conf(i)),
 274      &     -entfac(list_conf(i)),fdimless(i),sumprob
 275 #endif
 276           if (sumprob.gt.prob_limit) goto 122
 277 c          if (sumprob.gt.1.00d0) goto 122
 278           nlist=nlist+1
 279         enddo
 280   122   continue
 281 #ifdef MPI
 282       endif
 283       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 284      &   IERROR)
 285       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 286      &   IERROR)
 287 c      do iproc=0,nprocs
 288 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 289 c     &   " indend",indend(iproc)
 290 c      enddo
 291       write (iout,*) "nlist",nlist
 292 #endif
 293       return
 294       end
 295 !--------------------------------------------------
 296       subroutine mysort1(n, x, ipermut)
 297       implicit none
 298       integer i,j,imax,ipm,n
 299       real x(n)
 300       integer ipermut(n)
 301       real xtemp
 302       do i=1,n
 303         xtemp=x(i)
 304         imax=i
 305         do j=i+1,n
 306           if (x(j).lt.xtemp) then
 307             imax=j
 308             xtemp=x(j)
 309           endif
 310         enddo
 311         x(imax)=x(i)
 312         x(i)=xtemp
 313         ipm=ipermut(imax)
 314         ipermut(imax)=ipermut(i)
 315         ipermut(i)=ipm
 316       enddo
 317       return
 318       end