source/cluster/probabl.f90

   1       module probability
   2 !-----------------------------------------------------------------------------
   3       use clust_data
   4       implicit none
   5 !-----------------------------------------------------------------------------
   6 !
   7 !
   8 !-----------------------------------------------------------------------------
   9       contains
  10 !----------------------------------------------------------------------------
  11 ! probabl.f90
  12 !----------------------------------------------------------------------------
  13       subroutine probabl(ib,nlist,ncon,*)
  14 ! construct the conformational ensembles at REMD temperatures
  15 !      implicit none
  16 !      include "DIMENSIONS"
  17 !      include "sizesclu.dat"
  18       use io_units
  19       use io_clust, only: daread_ccoords
  20       use geometry_data, only: nres,c
  21       use energy_data, only: nss,ihpb,jhpb,max_ene
  22       use geometry, only: int_from_cart1
  23       use energy, only: etotal,rescale_weights
  24       use energy, only: rescale_mode,enerprint,weights
  25 #ifdef MPI
  26       use MPI_data
  27       include "mpif.h"
  28 !      include "COMMON.MPI"
  29       integer :: ierror,errcode !,status(MPI_STATUS_SIZE)
  30 #endif
  31 !      include "COMMON.IOUNITS"
  32 !      include "COMMON.FREE"
  33 !      include "COMMON.FFIELD"
  34 !      include "COMMON.INTERACT"
  35 !      include "COMMON.SBRIDGE"
  36 !      include "COMMON.CHAIN"
  37 !      include "COMMON.CLUSTER"
  38       real(kind=4) :: csingle(3,2*nres)
  39       real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
  40         eprim,ebis,temper,kfac=2.4d0,T0=300.0d0
  41       real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,&
  42             ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
  43             eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,&
  44             evdw_t
  45       integer :: i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
  46       real(kind=8) :: qfree,sumprob,eini,efree,rmsdev
  47       character(len=80) :: bxname
  48       character(len=2) :: licz1
  49       character(len=5) :: ctemper
  50 !      integer ilen
  51 !      external ilen
  52       real(kind=4) :: Fdimless(maxconf),Fdimless_(maxconf)
  53       real(kind=8) :: totfree_(0:maxconf),entfac_(maxconf)
  54       real(kind=8) :: energia(0:max_ene)
  55       integer,dimension(0:nprocs) :: scount_
  56
  57       do i=1,ncon
  58         list_conf(i)=i
  59       enddo
  60 !      do i=1,ncon
  61 !        write (iout,*) i,list_conf(i)
  62 !      enddo
  63 #ifdef MPI
  64       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
  65       call daread_ccoords(indstart(me),indend(me))
  66 #endif
  67 !      write (iout,*) "ncon",ncon
  68       temper=1.0d0/(beta_h(ib)*1.987D-3)
  69 !elwrite(iout,*)"rescale_mode", rescale_mode
  70 !      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
  71 !      quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  72 !      quotl=1.0d0
  73 !      kfacl=1.0d0
  74 !      do l=1,5
  75 !        quotl1=quotl
  76 !        quotl=quotl*quot
  77 !        kfacl=kfacl*kfac
  78 !        fT(l)=kfacl/(kfacl-1.0d0+quotl)
  79 !      enddo
  80 ! EL start old rescale-------------------------------
  81 !           if (rescale_mode.eq.1) then
  82 !             quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
  83 !             quotl=1.0d0
  84 !             kfacl=1.0d0
  85 !             do l=1,5
  86 !               quotl1=quotl
  87 !               quotl=quotl*quot
  88 !               kfacl=kfacl*kfac
  89 !               fT(l)=kfacl/(kfacl-1.0d0+quotl)
  90 !             enddo
  91 !if defined(FUNCTH)
  92 !             ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
  93 !                       320.0d0
  94 !             ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
  95 !            ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
  96 !                   /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
  97 !elif defined(FUNCT)
  98 !             fT(6)=betaT/T0
  99 !             ftprim(6)=1.0d0/T0
 100 !             ftbis(6)=0.0d0
 101 !else
 102 !             fT(6)=1.0d0
 103 !             ftprim(6)=0.0d0
 104 !             ftbis(6)=0.0d0
 105 !endif
 106 !
 107 !           else if (rescale_mode.eq.2) then
 108 !             quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
 109 !             quotl=1.0d0
 110 !             do l=1,5
 111 !               quotl=quotl*quot
 112 !               fT(l)=1.12692801104297249644d0/ &
 113 !                  dlog(dexp(quotl)+dexp(-quotl))
 114 !             enddo
 115 !el             write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3)
 116 !c             call flush(iout)
 117 !if defined(FUNCTH)
 118 !             ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
 119 !                       320.0d0
 120 !             ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 121 !            ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
 122 !                   /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 123 !elif defined(FUNCT)
 124 !             fT(6)=betaT/T0
 125 !             ftprim(6)=1.0d0/T0
 126 !             ftbis(6)=0.0d0
 127 !else
 128 !             fT(6)=1.0d0
 129 !             ftprim(6)=0.0d0
 130 !             ftbis(6)=0.0d0
 131 !endif
 132 !            endif
 133 !EL end old rescele----------------------
 134 !
 135         call rescale_weights(1.0d0/(beta_h(ib)*1.987D-3))
 136
 137 #ifdef MPI
 138       do i=1,scount(me)
 139         ii=i+indstart(me)-1
 140 #else
 141       do i=1,ncon
 142         ii=i
 143 #endif
 144 !        write (iout,*) "i",i," ii",ii
 145 !        call flush(iout)
 146         if (ib.eq.1) then
 147           do j=1,nres
 148             do k=1,3
 149               c(k,j)=allcart(k,j,i)
 150               c(k,j+nres)=allcart(k,j+nres,i)
 151             enddo
 152           enddo
 153           do k=1,3
 154             c(k,nres+1)=c(k,1)
 155             c(k,nres+nres)=c(k,nres)
 156           enddo
 157 !el          do j=1,2*nres
 158 !            do k=1,3
 159 !write(iout,*)"c, k, j",k,j,c(k,j)
 160 !            enddo
 161 !el          enddo
 162           nss=nss_all(i)
 163           do j=1,nss
 164             ihpb(j)=ihpb_all(j,i)
 165             jhpb(j)=jhpb_all(j,i)
 166           enddo
 167           call int_from_cart1(.false.)
 168 !          call etotal(energia(0),fT)
 169           call etotal(energia)
 170           totfree(i)=energia(0)
 171 !          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 172 !          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 173 !          call enerprint(energia)
 174 !          call pdbout(totfree(i),16,i)
 175 #ifdef DEBUG
 176           write (iout,*) i," energia",(energia(j),j=0,21)
 177           write (iout,*) "etot", etot
 178 !          write (iout,*) "ft(6)", ft(6)
 179 #endif
 180           do k=1,max_ene
 181             enetb(k,i)=energia(k)
 182           enddo
 183         endif
 184 !el        evdw=enetb(1,i)
 185 !        write (iout,*) evdw
 186        etot=energia(0)
 187 #ifdef SCP14
 188 !el        evdw2_14=enetb(17,i)
 189         evdw2_14=enetb(18,i)
 190         evdw2=enetb(2,i)+evdw2_14
 191 #else
 192         evdw2=enetb(2,i)
 193         evdw2_14=0.0d0
 194 #endif
 195 #ifdef SPLITELE
 196         ees=enetb(3,i)
 197         evdw1=enetb(16,i)
 198 #else
 199         ees=enetb(3,i)
 200         evdw1=0.0d0
 201 #endif
 202         ecorr=enetb(4,i)
 203         ecorr5=enetb(5,i)
 204         ecorr6=enetb(6,i)
 205         eel_loc=enetb(7,i)
 206         eello_turn3=enetb(8,i)
 207         eello_turn4=enetb(9,i)
 208         eturn6=enetb(10,i)
 209         ebe=enetb(11,i)
 210         escloc=enetb(12,i)
 211         etors=enetb(13,i)
 212         etors_d=enetb(14,i)
 213         ehpb=enetb(15,i)
 214 !        estr=enetb(18,i)
 215         estr=enetb(17,i)
 216 !        esccor=enetb(19,i)
 217         esccor=enetb(21,i)
 218 !        edihcnstr=enetb(20,i)
 219         edihcnstr=enetb(19,i)
 220 !#ifdef SPLITELE
 221 !        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
 222 !     &ft(1)*welec*ees+wvdwpp*evdw1
 223 !     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 224 !     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 225 !     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 226 !     &  +ft(2)*wturn3*eello_turn3
 227 !     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 228 !     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 229 !     &  +wbond*estr
 230 !#else
 231 !        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
 232 !     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 233 !     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 234 !     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 235 !     &  +ft(2)*wturn3*eello_turn3
 236 !     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 237 !     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 238 !     &  +wbond*estr
 239 !#endif
 240 !#ifdef DEBUG
 241 !        write (iout,*) "etot2", etot
 242 !        write (iout,*) "evdw", wsc, evdw,evdw_t
 243 !        write (iout,*) "evdw2", wscp, evdw2
 244 !        write (iout,*) "welec", ft(1),welec,ees
 245 !        write (iout,*) "evdw1", wvdwpp,evdw1
 246 !        write (iout,*) "ebe",ebe,wang
 247 !#endif
 248         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 249         totfree(i)=etot
 250 #ifdef DEBUG
 251         write (iout,*) "fdim calc", i,ii,ib,&
 252          1.0d0/(1.987d-3*beta_h(ib)),totfree(i),&
 253          entfac(ii),Fdimless(i)
 254 #endif
 255         Fdimless_(i)=Fdimless(i)
 256         totfree_(i)=totfree(i)
 257           call enerprint(energia(0)) !el
 258       enddo   ! i
 259
 260       do i=1,maxconf
 261         entfac_(i)=entfac(i)
 262       enddo
 263       do i=0,nprocs
 264         scount_(i)=scount(i)
 265       enddo
 266
 267 #ifdef MPI
 268       call MPI_Gatherv(Fdimless_(1),scount_(me),&
 269        MPI_REAL,Fdimless(1),&
 270        scount(0),idispl(0),MPI_REAL,Master,&
 271        MPI_COMM_WORLD, IERROR)
 272       call MPI_Gatherv(totfree_(1),scount(me),&
 273        MPI_DOUBLE_PRECISION,totfree(1),&
 274        scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
 275        MPI_COMM_WORLD, IERROR)
 276       call MPI_Gatherv(entfac_(indstart(me)+1),scount_(me),&
 277        MPI_DOUBLE_PRECISION,entfac(1),&
 278        scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
 279        MPI_COMM_WORLD, IERROR)
 280
 281       if (me.eq.Master) then
 282 #endif
 283 #ifdef DEBUG
 284         write (iout,*) "The FDIMLESS array before sorting"
 285         do i=1,ncon
 286           write (iout,*) i,fdimless(i)
 287         enddo
 288 #endif
 289         call mysort1(ncon,Fdimless,list_conf)
 290 #ifdef DEBUG
 291         write (iout,*) "The FDIMLESS array after sorting"
 292         do i=1,ncon
 293           write (iout,*) i,list_conf(i),fdimless(i)
 294         enddo
 295 #endif
 296         do i=1,ncon
 297           totfree(i)=fdimless(i)
 298         enddo
 299         qfree=0.0d0
 300         do i=1,ncon
 301           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 302 !          write (iout,*) "fdimless", fdimless(i)
 303         enddo
 304         write (iout,*) "qfree",qfree !d
 305         nlist=1
 306         sumprob=0.0
 307         write (iout,*) "ncon", ncon,maxstr_proc
 308         do i=1,min0(ncon,maxstr_proc)-1
 309           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 310 #ifdef DEBUG
 311           write (iout,*) i,ib,beta_h(ib),&
 312            1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),&
 313            totfree(list_conf(i)),&
 314            -entfac(list_conf(i)),fdimless(i),sumprob
 315 #endif
 316           if (sumprob.gt.prob_limit) goto 122
 317 !          if (sumprob.gt.1.00d0) goto 122
 318           nlist=nlist+1
 319         enddo
 320   122   continue
 321 #ifdef MPI
 322       endif
 323       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, &
 324          IERROR)
 325       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,&
 326          IERROR)
 327 !      do iproc=0,nprocs
 328 !        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 329 !     &   " indend",indend(iproc)
 330 !      enddo
 331 #endif
 332 !write(iout,*)"koniec probabl"
 333       return
 334       end subroutine probabl
 335 !--------------------------------------------------
 336       subroutine mysort1(n, x, ipermut)
 337 !      implicit none
 338       integer :: i,j,imax,ipm,n
 339       real(kind=4) :: x(n)
 340       integer :: ipermut(n)
 341       real(kind=4) :: xtemp
 342       do i=1,n
 343         xtemp=x(i)
 344         imax=i
 345         do j=i+1,n
 346           if (x(j).lt.xtemp) then
 347             imax=j
 348             xtemp=x(j)
 349           endif
 350         enddo
 351         x(imax)=x(i)
 352         x(i)=xtemp
 353         ipm=ipermut(imax)
 354         ipermut(imax)=ipermut(i)
 355         ipermut(i)=ipm
 356       enddo
 357       return
 358       end subroutine mysort1
 359 !-----------------------------------------------------------------------------
 360 !-----------------------------------------------------------------------------
 361       end module probability