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