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