source/cluster/wham/src/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(5),fTprim(5),fTbis(5),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       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
  25       double precision qfree,sumprob,eini,efree,rmsdev
  26       character*80 bxname
  27       character*2 licz1
  28       character*5 ctemper
  29       integer ilen
  30       external ilen
  31       real*4 Fdimless(maxconf)
  32       double precision energia(0:max_ene)
  33       do i=1,ncon
  34         list_conf(i)=i
  35       enddo
  36 c      do i=1,ncon
  37 c        write (iout,*) i,list_conf(i)
  38 c      enddo
  39 #ifdef MPI
  40       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  41       call daread_ccoords(indstart(me),indend(me))
  42 #endif
  43 c      write (iout,*) "ncon",ncon
  44       temper=1.0d0/(beta_h(ib)*1.987D-3)
  45 c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
  46 c      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  47 c      quotl=1.0d0
  48 c      kfacl=1.0d0
  49 c      do l=1,5
  50 c        quotl1=quotl
  51 c        quotl=quotl*quot
  52 c        kfacl=kfacl*kfac
  53 c        fT(l)=kfacl/(kfacl-1.0d0+quotl)
  54 c      enddo
  55             if (rescale_mode.eq.1) then
  56               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  57               quotl=1.0d0
  58               kfacl=1.0d0
  59               do l=1,5
  60                 quotl1=quotl
  61                 quotl=quotl*quot
  62                 kfacl=kfacl*kfac
  63                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
  64               enddo
  65             else if (rescale_mode.eq.2) then
  66               quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  67               quotl=1.0d0
  68               do l=1,5
  69                 quotl=quotl*quot
  70                 fT(l)=1.12692801104297249644d0/
  71      &             dlog(dexp(quotl)+dexp(-quotl))
  72               enddo
  73 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
  74               call flush(iout)
  75             endif
  76
  77 #ifdef MPI
  78       do i=1,scount(me)
  79         ii=i+indstart(me)-1
  80 #else
  81       do i=1,ncon
  82         ii=i
  83 #endif
  84 c        write (iout,*) "i",i," ii",ii
  85         call flush(iout)
  86         if (ib.eq.1) then
  87           do j=1,nres
  88             do k=1,3
  89               c(k,j)=allcart(k,j,i)
  90               c(k,j+nres)=allcart(k,j+nres,i)
  91             enddo
  92           enddo
  93           do k=1,3
  94             c(k,nres+1)=c(k,1)
  95             c(k,nres+nres)=c(k,nres)
  96           enddo
  97           nss=nss_all(i)
  98           do j=1,nss
  99             ihpb(j)=ihpb_all(j,i)
 100             jhpb(j)=jhpb_all(j,i)
 101           enddo
 102           call int_from_cart1(.false.)
 103           call etotal(energia(0),fT)
 104           totfree(i)=energia(0)
 105 #ifdef DEBUG
 106           write (iout,*) i," energia",(energia(j),j=0,19)
 107 #endif
 108           do k=1,max_ene
 109             enetb(k,i)=energia(k)
 110           enddo
 111         endif
 112         evdw=enetb(1,i)
 113 #ifdef SCP14
 114         evdw2_14=enetb(17,i)
 115         evdw2=enetb(2,i)+evdw2_14
 116 #else
 117         evdw2=enetb(2,i)
 118         evdw2_14=0.0d0
 119 #endif
 120 #ifdef SPLITELE
 121         ees=enetb(3,i)
 122         evdw1=enetb(16,i)
 123 #else
 124         ees=enetb(3,i)
 125         evdw1=0.0d0
 126 #endif
 127         ecorr=enetb(4,i)
 128         ecorr5=enetb(5,i)
 129         ecorr6=enetb(6,i)
 130         eel_loc=enetb(7,i)
 131         eello_turn3=enetb(8,i)
 132         eello_turn4=enetb(9,i)
 133         eturn6=enetb(10,i)
 134         ebe=enetb(11,i)
 135         escloc=enetb(12,i)
 136         etors=enetb(13,i)
 137         etors_d=enetb(14,i)
 138         ehpb=enetb(15,i)
 139         estr=enetb(18,i)
 140         esccor=enetb(19,i)
 141         edihcnstr=enetb(20,i)
 142 #ifdef SPLITELE
 143         etot=wsc*evdw+wscp*evdw2+ft(1)*welec*ees+wvdwpp*evdw1
 144      &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 145      &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 146      &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 147      &  +ft(2)*wturn3*eello_turn3
 148      &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 149      &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 150      &  +wbond*estr
 151 #else
 152         etot=wsc*evdw+wscp*evdw2+ft(1)*welec*(ees+evdw1)
 153      &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 154      &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 155      &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 156      &  +ft(2)*wturn3*eello_turn3
 157      &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 158      &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 159      &  +wbond*estr
 160 #endif
 161         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 162         totfree(i)=etot
 163 #ifdef DEBUG
 164         write (iout,*) i,ii,ib,
 165      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
 166      &   entfac(ii),Fdimless(i)
 167 #endif
 168       enddo   ! i
 169 #ifdef MPI
 170       call MPI_Gatherv(Fdimless(1),scount(me),
 171      & MPI_REAL,Fdimless(1),
 172      & scount(0),idispl(0),MPI_REAL,Master,
 173      & MPI_COMM_WORLD, IERROR)
 174       call MPI_Gatherv(totfree(1),scount(me),
 175      & MPI_DOUBLE_PRECISION,totfree(1),
 176      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 177      & MPI_COMM_WORLD, IERROR)
 178       call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
 179      & MPI_DOUBLE_PRECISION,entfac(1),
 180      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 181      & MPI_COMM_WORLD, IERROR)
 182       if (me.eq.Master) then
 183 #endif
 184 #ifdef DEBUG
 185         write (iout,*) "The FDIMLESS array before sorting"
 186         do i=1,ncon
 187           write (iout,*) i,fdimless(i)
 188         enddo
 189 #endif
 190         call mysort1(ncon,Fdimless,list_conf)
 191 #ifdef DEBUG
 192         write (iout,*) "The FDIMLESS array after sorting"
 193         do i=1,ncon
 194           write (iout,*) i,list_conf(i),fdimless(i)
 195         enddo
 196 #endif
 197         do i=1,ncon
 198           totfree(i)=fdimless(i)
 199         enddo
 200         qfree=0.0d0
 201         do i=1,ncon
 202           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 203         enddo
 204 c        write (iout,*) "qfree",qfree
 205         nlist=1
 206         sumprob=0.0
 207         do i=1,min0(ncon,maxstr_proc)-1
 208           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 209 #ifdef DEBUG
 210           write (iout,*) i,ib,beta_h(ib),
 211      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 212      &     totfree(list_conf(i)),
 213      &     -entfac(list_conf(i)),fdimless(i),sumprob
 214 #endif
 215           if (sumprob.gt.prob_limit) goto 122
 216 c          if (sumprob.gt.1.00d0) goto 122
 217           nlist=nlist+1
 218         enddo
 219   122   continue
 220 #ifdef MPI
 221       endif
 222       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 223      &   IERROR)
 224       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 225      &   IERROR)
 226 c      do iproc=0,nprocs
 227 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 228 c     &   " indend",indend(iproc)
 229 c      enddo
 230       write (iout,*) "nlist",nlist
 231 #endif
 232       return
 233       end
 234 !--------------------------------------------------
 235       subroutine mysort1(n, x, ipermut)
 236       implicit none
 237       integer i,j,imax,ipm,n
 238       real x(n)
 239       integer ipermut(n)
 240       real xtemp
 241       do i=1,n
 242         xtemp=x(i)
 243         imax=i
 244         do j=i+1,n
 245           if (x(j).lt.xtemp) then
 246             imax=j
 247             xtemp=x(j)
 248           endif
 249         enddo
 250         x(imax)=x(i)
 251         x(i)=xtemp
 252         ipm=ipermut(imax)
 253         ipermut(imax)=ipermut(i)
 254         ipermut(i)=ipm
 255       enddo
 256       return
 257       end