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         write(iout,*) "nss",nss
 249         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 250         totfree(i)=etot
 251 #ifdef DEBUG
 252         write (iout,*) "fdim calc", i,ii,ib,&
 253          1.0d0/(1.987d-3*beta_h(ib)),totfree(i),&
 254          entfac(ii),Fdimless(i)
 255 #endif
 256         Fdimless_(i)=Fdimless(i)
 257         totfree_(i)=totfree(i)
 258           call enerprint(energia(0)) !el
 259       enddo   ! i
 260
 261       do i=1,maxconf
 262         entfac_(i)=entfac(i)
 263       enddo
 264       do i=0,nprocs
 265         scount_(i)=scount(i)
 266       enddo
 267
 268 #ifdef MPI
 269       call MPI_Gatherv(Fdimless_(1),scount_(me),&
 270        MPI_REAL,Fdimless(1),&
 271        scount(0),idispl(0),MPI_REAL,Master,&
 272        MPI_COMM_WORLD, IERROR)
 273       call MPI_Gatherv(totfree_(1),scount(me),&
 274        MPI_DOUBLE_PRECISION,totfree(1),&
 275        scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
 276        MPI_COMM_WORLD, IERROR)
 277       call MPI_Gatherv(entfac_(indstart(me)+1),scount_(me),&
 278        MPI_DOUBLE_PRECISION,entfac(1),&
 279        scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,&
 280        MPI_COMM_WORLD, IERROR)
 281
 282       if (me.eq.Master) then
 283 #endif
 284 #ifdef DEBUG
 285         write (iout,*) "The FDIMLESS array before sorting"
 286         do i=1,ncon
 287           write (iout,*) i,fdimless(i)
 288         enddo
 289 #endif
 290         call mysort1(ncon,Fdimless,list_conf)
 291 #ifdef DEBUG
 292         write (iout,*) "The FDIMLESS array after sorting"
 293         do i=1,ncon
 294           write (iout,*) i,list_conf(i),fdimless(i)
 295         enddo
 296 #endif
 297         do i=1,ncon
 298           totfree(i)=fdimless(i)
 299         enddo
 300         qfree=0.0d0
 301         do i=1,ncon
 302           qfree=qfree+exp(-fdimless(i)+fdimless(1))
 303 !          write (iout,*) "fdimless", fdimless(i)
 304         enddo
 305         write (iout,*) "qfree",qfree !d
 306         nlist=1
 307         sumprob=0.0
 308         write (iout,*) "ncon", ncon,maxstr_proc
 309         do i=1,min0(ncon,maxstr_proc)-1
 310           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 311 #ifdef DEBUG
 312           write (iout,*) i,ib,beta_h(ib),&
 313            1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),&
 314            totfree(list_conf(i)),&
 315            -entfac(list_conf(i)),fdimless(i),sumprob
 316 #endif
 317           if (sumprob.gt.prob_limit) goto 122
 318 !          if (sumprob.gt.1.00d0) goto 122
 319           nlist=nlist+1
 320         enddo
 321   122   continue
 322 #ifdef MPI
 323       endif
 324       call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, &
 325          IERROR)
 326       call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,&
 327          IERROR)
 328 !      do iproc=0,nprocs
 329 !        write (iout,*) "iproc",iproc," indstart",indstart(iproc),
 330 !     &   " indend",indend(iproc)
 331 !      enddo
 332 #endif
 333 !write(iout,*)"koniec probabl"
 334       return
 335       end subroutine probabl
 336 !--------------------------------------------------
 337       subroutine mysort1(n, x, ipermut)
 338 !      implicit none
 339       integer :: i,j,imax,ipm,n
 340       real(kind=4) :: x(n)
 341       integer :: ipermut(n)
 342       real(kind=4) :: xtemp
 343       do i=1,n
 344         xtemp=x(i)
 345         imax=i
 346         do j=i+1,n
 347           if (x(j).lt.xtemp) then
 348             imax=j
 349             xtemp=x(j)
 350           endif
 351         enddo
 352         x(imax)=x(i)
 353         x(i)=xtemp
 354         ipm=ipermut(imax)
 355         ipermut(imax)=ipermut(i)
 356         ipermut(i)=ipm
 357       enddo
 358       return
 359       end subroutine mysort1
 360 !-----------------------------------------------------------------------------
 361 !-----------------------------------------------------------------------------
 362       end module probability