source/wham/src/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,ehomology_constr,
  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             ehomology_constr=enetb(22,i,iparm)
 166             if (homol_nset.gt.1)
 167      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
 168 #ifdef SPLITELE
 169             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 170      &      +wvdwpp*evdw1
 171      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 172      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 173      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 174      &      +ft(2)*wturn3*eello_turn3
 175      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 176      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 177      &      +wbond*estr+ehomology_constr
 178 #else
 179             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 180      &      +ft(1)*welec*(ees+evdw1)
 181      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 182      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 183      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 184      &      +ft(2)*wturn3*eello_turn3
 185      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 186      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 187      &      +wbond*estr+ehomology_constr
 188 #endif
 189 #ifdef MPI
 190             Fdimless(i)=
 191      &        beta_h(ib,iparm)*etot-entfac(i)
 192             potE(i,iparm)=etot
 193 #ifdef DEBUG
 194             write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
 195      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
 196      &       -entfac(i),Fdimless(i)
 197 #endif
 198 #else
 199             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
 200             potE(i,iparm)=etot
 201 #endif
 202           enddo   ! i
 203 #ifdef MPI
 204           call MPI_Gatherv(Fdimless(1),scount(me),
 205      &     MPI_REAL,Fdimless(1),
 206      &     scount(0),idispl(0),MPI_REAL,Master,
 207      &     WHAM_COMM, IERROR)
 208 #ifdef DEBUG
 209           call MPI_Gatherv(potE(1,iparm),scount(me),
 210      &     MPI_DOUBLE_PRECISION,potE(1,iparm),
 211      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 212      &     WHAM_COMM, IERROR)
 213           call MPI_Gatherv(entfac(1),scount(me),
 214      &     MPI_DOUBLE_PRECISION,entfac(1),
 215      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 216      &     WHAM_COMM, IERROR)
 217 #endif
 218           if (me.eq.Master) then
 219 #ifdef DEBUG
 220           write (iout,*) "The FDIMLESS array before sorting"
 221           do i=1,ntot(islice)
 222             write (iout,*) i,fdimless(i)
 223           enddo
 224 #endif
 225 #endif
 226           do i=1,ntot(islice)
 227             iperm(i)=i
 228           enddo
 229           call mysort1(ntot(islice),Fdimless,iperm)
 230 #ifdef DEBUG
 231           write (iout,*) "The FDIMLESS array after sorting"
 232           do i=1,ntot(islice)
 233             write (iout,*) i,iperm(i),fdimless(i)
 234           enddo
 235 #endif
 236           qfree=0.0d0
 237           do i=1,ntot(islice)
 238             qfree=qfree+exp(-fdimless(i)+fdimless(1))
 239           enddo
 240 c          write (iout,*) "qfree",qfree
 241           nlist=1
 242           sumprob=0.0
 243           do i=1,min0(ntot(islice),ensembles)
 244             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 245 #ifdef DEBUG
 246             write (iout,*) i,ib,beta_h(ib,iparm),
 247      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
 248      &       potE(iperm(i),iparm),
 249      &       -entfac(iperm(i)),fdimless(i),sumprob
 250 #endif
 251             if (sumprob.gt.0.99d0) goto 122
 252             nlist=nlist+1
 253           enddo
 254   122     continue
 255 #ifdef MPI
 256           endif
 257           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
 258      &       IERROR)
 259           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
 260      &       IERROR)
 261           do i=1,nlist
 262             ii=iperm(i)
 263             iproc=0
 264             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
 265               iproc=iproc+1
 266             enddo
 267             if (iproc.ge.nprocs) then
 268               write (iout,*) "Fatal error: processor out of range",iproc
 269               call flush(iout)
 270               if (bxfile) then
 271                 close (ientout)
 272               else
 273                 close (ientout,status="delete")
 274               endif
 275               return1
 276             endif
 277             ik=ii-indstart(iproc)+1
 278             if (iproc.ne.Master) then
 279               if (me.eq.iproc) then
 280 #ifdef DEBUG
 281                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
 282      &           " energy",potE(ik,iparm)
 283 #endif
 284                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
 285      &            Master,i,WHAM_COMM,IERROR)
 286               else if (me.eq.Master) then
 287                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
 288      &            WHAM_COMM,STATUS,IERROR)
 289               endif
 290             else if (me.eq.Master) then
 291               enepot(i)=potE(ik,iparm)
 292             endif
 293           enddo
 294 #else
 295           do i=1,nlist
 296             enepot(i)=potE(iperm(i),iparm)
 297           enddo
 298 #endif
 299 #ifdef MPI
 300           if (me.eq.Master) then
 301 #endif
 302           write(licz3,'(bz,i3.3)') iparm
 303           write(licz2,'(bz,i2.2)') islice
 304           if (temper.lt.100.0d0) then
 305             write(ctemper,'(f3.0)') temper
 306           else if (temper.lt.1000.0) then
 307             write (ctemper,'(f4.0)') temper
 308           else
 309             write (ctemper,'(f5.0)') temper
 310           endif
 311           if (nparmset.eq.1) then
 312             if (separate_parset) then
 313               write(licz4,'(bz,i3.3)') myparm
 314               pdbname=prefix(:ilen(prefix))//"_par"//licz4
 315             else
 316               pdbname=prefix(:ilen(prefix))
 317             endif
 318           else
 319             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
 320           endif
 321           if (nslice.eq.1) then
 322             pdbname=pdbname(:ilen(pdbname))//"_T_"//
 323      &        ctemper(:ilen(ctemper))//"pdb"
 324           else
 325             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
 326      &        ctemper(:ilen(ctemper))//"pdb"
 327           endif
 328           open(ipdb,file=pdbname)
 329           do i=1,nlist
 330             read (ientout,rec=iperm(i))
 331      &        ((csingle(l,k),l=1,3),k=1,nres),
 332      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 333      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 334      &        eini,efree,rmsdev,iscor
 335             do j=1,2*nres
 336               do k=1,3
 337                 c(k,j)=csingle(k,j)
 338               enddo
 339             enddo
 340             eini=fdimless(i)
 341             call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
 342 cc          if (temper.eq.300.0d0) then
 343 cc          write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree
 344 cc          flush(iout)
 345 cc          endif
 346           enddo
 347 cc          if (temper.eq.300.0d0) then
 348 cc          write(iout,*) 'Gral i=',i ,eini,enepot(i),efree
 349 cc          flush(iout)
 350 cc          endif
 351 #ifdef MPI
 352           endif
 353 #endif
 354         enddo     ! ib
 355       enddo       ! iparm
 356       if (bxfile) then
 357         close(ientout)
 358       else
 359         close(ientout,status="delete")
 360       endif
 361       return
 362       end
 363 !--------------------------------------------------
 364       subroutine mysort1(n, x, ipermut)
 365       implicit none
 366       integer i,j,imax,ipm,n
 367       real x(n)
 368       integer ipermut(n)
 369       real xtemp
 370       do i=1,n
 371         xtemp=x(i)
 372         imax=i
 373         do j=i+1,n
 374           if (x(j).lt.xtemp) then
 375             imax=j
 376             xtemp=x(j)
 377           endif
 378         enddo
 379         x(imax)=x(i)
 380         x(i)=xtemp
 381         ipm=ipermut(imax)
 382         ipermut(imax)=ipermut(i)
 383         ipermut(i)=ipm
 384       enddo
 385       return
 386       end