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