source/cluster/wham/src-M-SAXS-homology/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.CONTROL"
  12       include "COMMON.IOUNITS"
  13       include "COMMON.FREE"
  14       include "COMMON.FFIELD"
  15       include "COMMON.INTERACT"
  16       include "COMMON.SBRIDGE"
  17       include "COMMON.CHAIN"
  18       include "COMMON.CLUSTER"
  19       real*4 csingle(3,maxres2)
  20       double precision fT(6),fTprim(6),fTbis(6),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      &      evdw_t,esaxs,eliptran,ethetacnstr,ehomology_constr,
  26      &      edfadis,edfator,edfanei,edfabet
  27       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
  28       double precision qfree,sumprob,eini,efree,rmsdev
  29       character*80 bxname
  30       character*2 licz1
  31       character*5 ctemper
  32       integer ilen,ijk
  33       external ilen
  34       character*80 structure/'Structure'/
  35       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
  36       double precision energia(0:max_ene), totfree_buf(0:maxconf),
  37      &  entfac_buf(maxconf)
  38       double precision buffer(maxconf)
  39       do i=1,ncon
  40         list_conf(i)=i
  41       enddo
  42 c      do i=1,ncon
  43 c        write (iout,*) i,list_conf(i)
  44 c      enddo
  45 #ifdef MPI
  46       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  47       call daread_ccoords(indstart(me),indend(me))
  48 #endif
  49 C      write (iout,*) "ncon",ncon
  50 C      call flush(iout)
  51       temper=1.0d0/(beta_h(ib)*1.987D-3)
  52       if (rescale_mode.eq.1) then
  53         quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  54         quotl=1.0d0
  55         kfacl=1.0d0
  56         do l=1,5
  57           quotl1=quotl
  58           quotl=quotl*quot
  59           kfacl=kfacl*kfac
  60           fT(l)=kfacl/(kfacl-1.0d0+quotl)
  61         enddo
  62 #if defined(FUNCTH)
  63         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  64      &             320.0d0
  65         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  66         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  67      &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  68 #elif defined(FUNCT)
  69         fT(6)=betaT/T0
  70         ftprim(6)=1.0d0/T0
  71         ftbis(6)=0.0d0
  72 #else
  73         fT(6)=1.0d0
  74         ftprim(6)=0.0d0
  75         ftbis(6)=0.0d0
  76 #endif
  77
  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          call flush(iout)
  88 #if defined(FUNCTH)
  89         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
  90      &             320.0d0
  91         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  92         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
  93      &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  94 #elif defined(FUNCT)
  95         fT(6)=betaT/T0
  96         ftprim(6)=1.0d0/T0
  97         ftbis(6)=0.0d0
  98 #else
  99         fT(6)=1.0d0
 100         ftprim(6)=0.0d0
 101         ftbis(6)=0.0d0
 102 #endif
 103       endif
 104
 105 #ifdef MPI
 106       do i=1,scount(me)
 107         ii=i+indstart(me)-1
 108 #else
 109       do i=1,ncon
 110         ii=i
 111 #endif
 112 C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
 113 c        call flush(iout)
 114 c        if (ib.eq.1) then
 115         do j=1,nres
 116           do k=1,3
 117             c(k,j)=allcart(k,j,i)
 118             c(k,j+nres)=allcart(k,j+nres,i)
 119 C              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
 120 C     &        c(k,j+nres),allcart(k,j+nres,i)
 121           enddo
 122         enddo
 123 C          write(iout,*) "out of j loop"
 124 C          call flush(iout)
 125         do k=1,3
 126           c(k,nres+1)=c(k,1)
 127           c(k,nres+nres)=c(k,nres)
 128         enddo
 129 C          write(iout,*) "after nres+nres",nss_all(i)
 130 C          call flush(iout)
 131         nss=nss_all(i)
 132         do j=1,nss
 133           ihpb(j)=ihpb_all(j,i)
 134           jhpb(j)=jhpb_all(j,i)
 135         enddo
 136         call int_from_cart1(.false.)
 137         call etotal(energia(0),fT)
 138         if (refstr) then
 139           write (structure(9:),'(bz,i6.6)') i
 140           call TMscore_sub(rmsdev,gdt_ts_tb(i),
 141      &    gdt_ha_tb(i),tmscore_tb(i),Structure,.false.)
 142 #define DEBUG
 143 #ifdef DEBUG
 144           write (iout,*) rmsdev,gdt_ts_tb(i),gdt_ha_tb(i),
 145      &      tmscore_tb(i)
 146 #endif
 147         endif
 148         totfree(i)=energia(0)
 149         totfree_buf(i)=totfree(i)
 150 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 151 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 152 c          call pdbout(totfree(i),16,i)
 153 c          call flush(iout)
 154 #ifdef DEBUG
 155         write (iout,*) "conformation", i
 156         call enerprint(energia(0),fT)
 157 #endif
 158         etot=energia(0)
 159         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 160         Fdimless_buf(i)=Fdimless(i)
 161         totfree(i)=etot
 162         totfree_buf(i)=totfree(i)
 163 #ifdef DEBUG
 164         write (iout,*) "fdim calc", 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
 170       do ijk=1,maxconf
 171       entfac_buf(ijk)=entfac(ijk)
 172       Fdimless_buf(ijk)=Fdimless(ijk)
 173       enddo
 174       do ijk=0,maxconf
 175       totfree_buf(ijk)=totfree(ijk)
 176       enddo
 177
 178
 179 c      scount_buf=scount(me)
 180 c      scount_buf2=scount(0)
 181
 182 c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
 183
 184 #ifdef MPI
 185 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
 186       call MPI_Gatherv(Fdimless_buf(1),scount(me),
 187      & MPI_REAL,Fdimless(1),
 188      & scount(0),idispl(0),MPI_REAL,Master,
 189      & MPI_COMM_WORLD, IERROR)
 190 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
 191       call MPI_Gatherv(totfree_buf(1),scount(me),
 192      & MPI_DOUBLE_PRECISION,totfree(1),
 193      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 194      & MPI_COMM_WORLD, IERROR)
 195 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
 196       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
 197      & MPI_DOUBLE_PRECISION,entfac(1),
 198      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 199      & MPI_COMM_WORLD, IERROR)
 200 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
 201       if (refstr) then
 202         do i=1,scount(me)
 203           buffer(i)=gdt_ts_tb(i)
 204         enddo
 205         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 206      &   gdt_ts_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 207      &   MPI_COMM_WORLD,IERROR)
 208         do i=1,scount(me)
 209           buffer(i)=gdt_ha_tb(i)
 210         enddo
 211         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 212      &   gdt_ha_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 213      &   MPI_COMM_WORLD,IERROR)
 214         do i=1,scount(me)
 215           buffer(i)=tmscore_tb(i)
 216         enddo
 217         call MPI_Gatherv(buffer(1),scount(me),MPI_DOUBLE_PRECISION,
 218      &   tmscore_tb(1),scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 219      &   MPI_COMM_WORLD,IERROR)
 220       endif
 221       if (me.eq.Master) then
 222 c      WRITE (iout,*) "me.eq.Master"
 223 #endif
 224 #ifdef DEBUG
 225         write (iout,*) "The FDIMLESS array before sorting"
 226         do i=1,ncon
 227           write (iout,*) i,fdimless(i)
 228         enddo
 229 #endif
 230 c      WRITE (iout,*) "Wchodze do call mysort1"
 231         call mysort1(ncon,Fdimless,list_conf)
 232 c      WRITE (iout,*) "Wychodze z call mysort1"
 233 #ifdef DEBUG
 234         write (iout,*) "The FDIMLESS array after sorting"
 235         do i=1,ncon
 236           write (iout,'(2i5,4f10.5)') i,list_conf(i),fdimless(i),
 237      &     gdt_ts_tb(i),gdt_ha_tb(i),tmscore_tb(i)
 238         enddo
 239 #endif
 240 c      WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)"
 241         do i=1,ncon
 242           totfree(i)=fdimless(i)
 243         enddo
 244         qfree=0.0d0
 245         do i=1,ncon
 246           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 247 c          write (iout,*) "fdimless", fdimless(i)
 248         enddo
 249 c        write (iout,*) "qfree",qfree
 250         nlist=1
 251         sumprob=0.0
 252         write (iout,*) "ncon", ncon,maxstr_proc
 253         do i=1,min0(ncon,maxstr_proc)-1
 254           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 255 #ifdef DEBUG
 256           write (iout,*) i,ib,beta_h(ib),
 257      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
 258      &     totfree(list_conf(i)),
 259      &     -entfac(list_conf(i)),fdimless(i),sumprob
 260 #endif
 261           if (sumprob.gt.prob_limit) goto 122
 262 c          if (sumprob.gt.1.00d0) goto 122
 263           nlist=nlist+1
 264         enddo
 265   122   continue
 266 #ifdef MPI
 267       endif
 268       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD,
 269      &   IERROR)
 270       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,
 271      &   IERROR)
 272 c      do iproc=0,nprocs
 273 c        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 274 c     &   " indend",indend(iproc)
 275 c      enddo
 276       write (iout,*) "nlist",nlist
 277 #endif
 278       return
 279       end
 280 !--------------------------------------------------
 281       subroutine mysort1(n, x, ipermut)
 282       implicit none
 283       integer i,j,imax,ipm,n
 284       real x(n)
 285       integer ipermut(n)
 286       real xtemp
 287       do i=1,n
 288         xtemp=x(i)
 289         imax=i
 290         do j=i+1,n
 291           if (x(j).lt.xtemp) then
 292             imax=j
 293             xtemp=x(j)
 294           endif
 295         enddo
 296         x(imax)=x(i)
 297         x(i)=xtemp
 298         ipm=ipermut(imax)
 299         ipermut(imax)=ipermut(i)
 300         ipermut(i)=ipm
 301       enddo
 302       return
 303       end