source/wham/src-M/make_ensemble1.F

   1       subroutine make_ensembles(islice,*)
   2 ! construct the conformational ensembles at REMD temperatures
   3       implicit none
   4       include "DIMENSIONS"
   5       include "DIMENSIONS.ZSCOPT"
   6       include "DIMENSIONS.FREE"
   7 #ifdef MPI
   8       include "mpif.h"
   9       include "COMMON.MPI"
  10       integer ierror,errcode,status(MPI_STATUS_SIZE)
  11 #endif
  12       include "COMMON.IOUNITS"
  13       include "COMMON.CONTROL"
  14       include "COMMON.FREE"
  15       include "COMMON.ENERGIES"
  16       include "COMMON.FFIELD"
  17       include "COMMON.INTERACT"
  18       include "COMMON.SBRIDGE"
  19       include "COMMON.CHAIN"
  20       include "COMMON.PROTFILES"
  21       include "COMMON.PROT"
  22       real*4 csingle(3,maxres2)
  23       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  24      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
  25       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  26      &      escloc,eliptran,etube,
  27      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  28      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
  29       integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
  30       double precision qfree,sumprob,eini,efree,rmsdev
  31       character*80 bxname
  32       character*2 licz1,licz2
  33       character*3 licz3,licz4
  34       character*5 ctemper
  35       integer ilen
  36       external ilen
  37       real*4 Fdimless(MaxStr)
  38       double precision enepot(MaxStr)
  39       integer iperm(MaxStr)
  40       integer islice
  41
  42 #ifdef MPI
  43       if (me.eq.Master) then
  44 #endif
  45       write (licz2,'(bz,i2.2)') islice
  46       if (nslice.eq.1) then
  47         if (.not.separate_parset) then
  48           bxname = prefix(:ilen(prefix))//".bx"
  49         else
  50           write (licz3,'(bz,i3.3)') myparm
  51           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
  52         endif
  53       else
  54         if (.not.separate_parset) then
  55           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
  56         else
  57           write (licz3,'(bz,i3.3)') myparm
  58           bxname = prefix(:ilen(prefix))//"par_"//licz3//
  59      &      "_slice_"//licz2//".bx"
  60         endif
  61       endif
  62       open (ientout,file=bxname,status="unknown",
  63      &  form="unformatted",access="direct",recl=lenrec1)
  64 #ifdef MPI
  65       endif
  66 #endif
  67       do iparm=1,nParmSet
  68         if (iparm.ne.iparmprint) exit
  69         call restore_parm(iparm)
  70         do ib=1,nT_h(iparm)
  71 #ifdef DEBUG
  72           write (iout,*) "iparm",iparm," ib",ib
  73 #endif
  74           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
  75 c          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
  76 c          quotl=1.0d0
  77 c          kfacl=1.0d0
  78 c          do l=1,5
  79 c            quotl1=quotl
  80 c            quotl=quotl*quot
  81 c            kfacl=kfacl*kfac
  82 c            fT(l)=kfacl/(kfacl-1.0d0+quotl)
  83 c          enddo
  84             if (rescale_mode.eq.1) then
  85               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
  86 #if defined(FUNCTH)
  87               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
  88               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
  89 #elif defined(FUNCT)
  90               ft(6)=quot
  91 #else
  92               ft(6)=1.0d0
  93 #endif
  94               quotl=1.0d0
  95               kfacl=1.0d0
  96               do l=1,5
  97                 quotl1=quotl
  98                 quotl=quotl*quot
  99                 kfacl=kfacl*kfac
 100                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 101               enddo
 102             else if (rescale_mode.eq.2) then
 103               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 104 #if defined(FUNCTH)
 105               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
 106               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
 107 #elif defined(FUNCT)
 108               ft(6)=quot
 109 #else
 110               ft(6)=1.0d0
 111 #endif
 112               quotl=1.0d0
 113               do l=1,5
 114                 quotl=quotl*quot
 115                 fT(l)=1.12692801104297249644d0/
 116      &             dlog(dexp(quotl)+dexp(-quotl))
 117               enddo
 118 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 119             else if (rescale_mode.eq.0) then
 120               do l=1,5
 121                 fT(l)=0.0d0
 122               enddo
 123             else
 124               write (iout,*)
 125      &        "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
 126               call flush(iout)
 127               return1
 128             endif
 129 #ifdef MPI
 130           do i=1,scount(me1)
 131 #else
 132           do i=1,ntot(islice)
 133 #endif
 134             evdw=enetb(1,i,iparm)
 135             evdw_t=enetb(21,i,iparm)
 136 #ifdef SCP14
 137             evdw2_14=enetb(17,i,iparm)
 138             evdw2=enetb(2,i,iparm)+evdw2_14
 139 #else
 140             evdw2=enetb(2,i,iparm)
 141             evdw2_14=0.0d0
 142 #endif
 143 #ifdef SPLITELE
 144             ees=enetb(3,i,iparm)
 145             evdw1=enetb(16,i,iparm)
 146 #else
 147             ees=enetb(3,i,iparm)
 148             evdw1=0.0d0
 149 #endif
 150             ecorr=enetb(4,i,iparm)
 151             ecorr5=enetb(5,i,iparm)
 152             ecorr6=enetb(6,i,iparm)
 153             eel_loc=enetb(7,i,iparm)
 154             eello_turn3=enetb(8,i,iparm)
 155             eello_turn4=enetb(9,i,iparm)
 156             eturn6=enetb(10,i,iparm)
 157             ebe=enetb(11,i,iparm)
 158             escloc=enetb(12,i,iparm)
 159             etors=enetb(13,i,iparm)
 160             etors_d=enetb(14,i,iparm)
 161             ehpb=enetb(15,i,iparm)
 162             estr=enetb(18,i,iparm)
 163             esccor=enetb(19,i,iparm)
 164             edihcnstr=enetb(20,i,iparm)
 165             eliptran=enetb(22,i,iparm)
 166             etube=enetb(25,i,iparm)
 167
 168 #ifdef SPLITELE
 169       if (shield_mode.gt.0) then
 170             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 171      &      +ft(1)*welec*ees
 172      &      +ft(1)*wvdwpp*evdw1
 173      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 174      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 175      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 176      &      +ft(2)*wturn3*eello_turn3
 177      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 178      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 179      &      +wbond*estr+wliptran*eliptran+wtube*etube
 180       else
 181             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 182      &      +wvdwpp*evdw1
 183      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 184      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 185      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 186      &      +ft(2)*wturn3*eello_turn3
 187      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 188      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 189      &      +wbond*estr+wliptran*eliptran+wtube*etube
 190       endif
 191 #else
 192       if (shield_mode.gt.0) then
 193             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 194      &      +ft(1)*welec*(ees+evdw1)
 195      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 196      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 197      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 198      &      +ft(2)*wturn3*eello_turn3
 199      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 200      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 201      &      +wbond*estr+wliptran*eliptran+wtube*etube
 202       else
 203             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 204      &      +ft(1)*welec*(ees+evdw1)
 205      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 206      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 207      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 208      &      +ft(2)*wturn3*eello_turn3
 209      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 210      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 211      &      +wbond*estr+wliptran*eliptran+wtube*etube
 212       endif
 213
 214 #endif
 215 #ifdef MPI
 216             Fdimless(i)=
 217      &        beta_h(ib,iparm)*etot-entfac(i)
 218             potE(i,iparm)=etot
 219 #ifdef DEBUG
 220             write (iout,*) i,indstart(me)+i-1,ib,
 221      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
 222      &       -entfac(i),Fdimless(i)
 223 #endif
 224 #else
 225             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
 226             potE(i,iparm)=etot
 227 #endif
 228           enddo   ! i
 229 #ifdef MPI
 230           call MPI_Gatherv(Fdimless(1),scount(me),
 231      &     MPI_REAL,Fdimless(1),
 232      &     scount(0),idispl(0),MPI_REAL,Master,
 233      &     WHAM_COMM, IERROR)
 234 #ifdef DEBUG
 235           call MPI_Gatherv(potE(1,iparm),scount(me),
 236      &     MPI_DOUBLE_PRECISION,potE(1,iparm),
 237      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 238      &     WHAM_COMM, IERROR)
 239           call MPI_Gatherv(entfac(1),scount(me),
 240      &     MPI_DOUBLE_PRECISION,entfac(1),
 241      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 242      &     WHAM_COMM, IERROR)
 243 #endif
 244           if (me.eq.Master) then
 245 #ifdef DEBUG
 246           write (iout,*) "The FDIMLESS array before sorting"
 247           do i=1,ntot(islice)
 248             write (iout,*) i,fdimless(i)
 249           enddo
 250 #endif
 251 #endif
 252           do i=1,ntot(islice)
 253             iperm(i)=i
 254           enddo
 255           call mysort1(ntot(islice),Fdimless,iperm)
 256 #ifdef DEBUG
 257           write (iout,*) "The FDIMLESS array after sorting"
 258           do i=1,ntot(islice)
 259             write (iout,*) i,iperm(i),fdimless(i)
 260           enddo
 261 #endif
 262           qfree=0.0d0
 263           do i=1,ntot(islice)
 264             qfree=qfree+exp(-fdimless(i)+fdimless(1))
 265           enddo
 266 c          write (iout,*) "qfree",qfree
 267           nlist=1
 268           sumprob=0.0
 269           do i=1,min0(ntot(islice),ensembles)
 270             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 271 #ifdef DEBUG
 272             write (iout,*) i,ib,beta_h(ib,iparm),
 273      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
 274      &       potE(iperm(i),iparm),
 275      &       -entfac(iperm(i)),fdimless(i),sumprob
 276 #endif
 277             if (sumprob.gt.0.99d0) goto 122
 278             nlist=nlist+1
 279           enddo
 280   122     continue
 281 #ifdef MPI
 282           endif
 283           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
 284      &       IERROR)
 285           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
 286      &       IERROR)
 287           do i=1,nlist
 288             ii=iperm(i)
 289             iproc=0
 290             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
 291               iproc=iproc+1
 292             enddo
 293             if (iproc.ge.nprocs) then
 294               write (iout,*) "Fatal error: processor out of range",iproc
 295               call flush(iout)
 296               if (bxfile) then
 297                 close (ientout)
 298               else
 299                 close (ientout,status="delete")
 300               endif
 301               return1
 302             endif
 303             ik=ii-indstart(iproc)+1
 304             if (iproc.ne.Master) then
 305               if (me.eq.iproc) then
 306 #ifdef DEBUG
 307                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
 308      &           " energy",potE(ik,iparm)
 309 #endif
 310                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
 311      &            Master,i,WHAM_COMM,IERROR)
 312               else if (me.eq.Master) then
 313                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
 314      &            WHAM_COMM,STATUS,IERROR)
 315               endif
 316             else if (me.eq.Master) then
 317               enepot(i)=potE(ik,iparm)
 318             endif
 319           enddo
 320 #else
 321           do i=1,nlist
 322             enepot(i)=potE(iperm(i),iparm)
 323           enddo
 324 #endif
 325 #ifdef MPI
 326           if (me.eq.Master) then
 327 #endif
 328           write(licz3,'(bz,i3.3)') iparm
 329           write(licz2,'(bz,i2.2)') islice
 330           if (temper.lt.100.0d0) then
 331             write(ctemper,'(f3.0)') temper
 332           else if (temper.lt.1000.0) then
 333             write (ctemper,'(f4.0)') temper
 334           else
 335             write (ctemper,'(f5.0)') temper
 336           endif
 337           if (nparmset.eq.1) then
 338             if (separate_parset) then
 339               write(licz4,'(bz,i3.3)') myparm
 340               pdbname=prefix(:ilen(prefix))//"_par"//licz4
 341             else
 342               pdbname=prefix(:ilen(prefix))
 343             endif
 344           else
 345             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
 346           endif
 347           if (nslice.eq.1) then
 348             pdbname=pdbname(:ilen(pdbname))//"_T_"//
 349      &        ctemper(:ilen(ctemper))//"pdb"
 350           else
 351             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
 352      &        ctemper(:ilen(ctemper))//"pdb"
 353           endif
 354           open(ipdb,file=pdbname)
 355           do i=1,nlist
 356             read (ientout,rec=iperm(i))
 357      &        ((csingle(l,k),l=1,3),k=1,nres),
 358      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 359      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 360      &        eini,efree,rmsdev,iscor
 361             do j=1,2*nres
 362               do k=1,3
 363                 c(k,j)=csingle(k,j)
 364               enddo
 365             enddo
 366             eini=fdimless(i)
 367             call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
 368           enddo
 369 #ifdef MPI
 370           endif
 371 #endif
 372         enddo     ! ib
 373       enddo       ! iparm
 374       if (bxfile) then
 375         close(ientout)
 376       else
 377         close(ientout,status="delete")
 378       endif
 379       return
 380       end
 381 !--------------------------------------------------
 382       subroutine mysort1(n, x, ipermut)
 383       implicit none
 384       integer i,j,imax,ipm,n
 385       real x(n)
 386       integer ipermut(n)
 387       real xtemp
 388       do i=1,n
 389         xtemp=x(i)
 390         imax=i
 391         do j=i+1,n
 392           if (x(j).lt.xtemp) then
 393             imax=j
 394             xtemp=x(j)
 395           endif
 396         enddo
 397         x(imax)=x(i)
 398         x(i)=xtemp
 399         ipm=ipermut(imax)
 400         ipermut(imax)=ipermut(i)
 401         ipermut(i)=ipm
 402       enddo
 403       return
 404       end