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