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