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