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      &      edfadis,edfator,edfanei,edfabet
  30       integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
  31       double precision qfree,sumprob,eini,efree,rmsdev
  32       character*80 bxname
  33       character*2 licz1,licz2
  34       character*3 licz3,licz4
  35       character*5 ctemper
  36       integer ilen
  37       external ilen
  38       real*4 Fdimless(MaxStr),Fdimless_(MaxStr)
  39       double precision enepot(MaxStr)
  40       integer iperm(MaxStr)
  41       integer islice
  42
  43 #ifdef MPI
  44       if (me.eq.Master) then
  45 #endif
  46       write (licz2,'(bz,i2.2)') islice
  47       if (nslice.eq.1) then
  48         if (.not.separate_parset) then
  49           bxname = prefix(:ilen(prefix))//".bx"
  50         else
  51           write (licz3,'(bz,i3.3)') myparm
  52           bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
  53         endif
  54       else
  55         if (.not.separate_parset) then
  56           bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
  57         else
  58           write (licz3,'(bz,i3.3)') myparm
  59           bxname = prefix(:ilen(prefix))//"par_"//licz3//
  60      &      "_slice_"//licz2//".bx"
  61         endif
  62       endif
  63       open (ientout,file=bxname,status="unknown",
  64      &  form="unformatted",access="direct",recl=lenrec1)
  65 #ifdef MPI
  66       endif
  67 #endif
  68       do iparm=1,nParmSet
  69         if (iparm.ne.iparmprint) exit
  70         call restore_parm(iparm)
  71         do ib=1,nT_h(iparm)
  72 #ifdef DEBUG
  73           write (iout,*) "iparm",iparm," ib",ib
  74 #endif
  75           temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
  76 c          quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
  77 c          quotl=1.0d0
  78 c          kfacl=1.0d0
  79 c          do l=1,5
  80 c            quotl1=quotl
  81 c            quotl=quotl*quot
  82 c            kfacl=kfacl*kfac
  83 c            fT(l)=kfacl/(kfacl-1.0d0+quotl)
  84 c          enddo
  85             if (rescale_mode.eq.1) then
  86               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
  87 #if defined(FUNCTH)
  88               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
  89               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
  90 #elif defined(FUNCT)
  91               ft(6)=quot
  92 #else
  93               ft(6)=1.0d0
  94 #endif
  95               quotl=1.0d0
  96               kfacl=1.0d0
  97               do l=1,5
  98                 quotl1=quotl
  99                 quotl=quotl*quot
 100                 kfacl=kfacl*kfac
 101                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 102               enddo
 103             else if (rescale_mode.eq.2) then
 104               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 105 #if defined(FUNCTH)
 106               tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
 107               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
 108 #elif defined(FUNCT)
 109               ft(6)=quot
 110 #else
 111               ft(6)=1.0d0
 112 #endif
 113               quotl=1.0d0
 114               do l=1,5
 115                 quotl=quotl*quot
 116                 fT(l)=1.12692801104297249644d0/
 117      &             dlog(dexp(quotl)+dexp(-quotl))
 118               enddo
 119 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 120             else if (rescale_mode.eq.0) then
 121               do l=1,5
 122                 fT(l)=0.0d0
 123               enddo
 124             else
 125               write (iout,*)
 126      &        "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
 127               call flush(iout)
 128               return1
 129             endif
 130 #ifdef MPI
 131           do i=1,scount(me1)
 132 #else
 133           do i=1,ntot(islice)
 134 #endif
 135             evdw=enetb(1,i,iparm)
 136             evdw_t=enetb(21,i,iparm)
 137 #ifdef SCP14
 138             evdw2_14=enetb(17,i,iparm)
 139             evdw2=enetb(2,i,iparm)+evdw2_14
 140 #else
 141             evdw2=enetb(2,i,iparm)
 142             evdw2_14=0.0d0
 143 #endif
 144 #ifdef SPLITELE
 145             ees=enetb(3,i,iparm)
 146             evdw1=enetb(16,i,iparm)
 147 #else
 148             ees=enetb(3,i,iparm)
 149             evdw1=0.0d0
 150 #endif
 151             ecorr=enetb(4,i,iparm)
 152             ecorr5=enetb(5,i,iparm)
 153             ecorr6=enetb(6,i,iparm)
 154             eel_loc=enetb(7,i,iparm)
 155             eello_turn3=enetb(8,i,iparm)
 156             eello_turn4=enetb(9,i,iparm)
 157             eturn6=enetb(10,i,iparm)
 158             ebe=enetb(11,i,iparm)
 159             escloc=enetb(12,i,iparm)
 160             etors=enetb(13,i,iparm)
 161             etors_d=enetb(14,i,iparm)
 162             ehpb=enetb(15,i,iparm)
 163             estr=enetb(18,i,iparm)
 164             esccor=enetb(19,i,iparm)
 165             edihcnstr=enetb(20,i,iparm)
 166             ehomology_constr=enetb(22,i,iparm)
 167             edfadis=enetb(23,i,iparm)
 168             edfator=enetb(24,i,iparm)
 169             edfanei=enetb(25,i,iparm)
 170             edfabet=enetb(26,i,iparm)
 171             if (homol_nset.gt.1)
 172      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
 173 #ifdef SPLITELE
 174             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 175      &      +wvdwpp*evdw1
 176      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 177      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 178      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 179      &      +ft(2)*wturn3*eello_turn3
 180      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 181      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 182      &      +wbond*estr+ehomology_constr
 183      &      +wdfa_dist*edfadis
 184      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 185 #else
 186             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 187      &      +ft(1)*welec*(ees+evdw1)
 188      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 189      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 190      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 191      &      +ft(2)*wturn3*eello_turn3
 192      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 193      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 194      &      +wbond*estr+ehomology_constr
 195      &      +wdfa_dist*edfadis
 196      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 197 #endif
 198 #ifdef MPI
 199             Fdimless_(i)=
 200      &        beta_h(ib,iparm)*etot-entfac(i)
 201             potE(i,iparm)=etot
 202 #ifdef DEBUG
 203             write (iout,*) 'EEE',i,indstart(me)+i-1,ib,
 204      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
 205      &       -entfac(i),Fdimless(i)
 206 #endif
 207 #else
 208             Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
 209             potE(i,iparm)=etot
 210 #endif
 211           enddo   ! i
 212 #ifdef MPI
 213           call MPI_Gatherv(Fdimless_(1),scount(me),
 214      &     MPI_REAL,Fdimless(1),
 215      &     scount(0),idispl(0),MPI_REAL,Master,
 216      &     WHAM_COMM, IERROR)
 217 #ifdef DEBUG
 218           call MPI_Gatherv(potE(1,iparm),scount(me),
 219      &     MPI_DOUBLE_PRECISION,potE(1,iparm),
 220      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 221      &     WHAM_COMM, IERROR)
 222           call MPI_Gatherv(entfac(1),scount(me),
 223      &     MPI_DOUBLE_PRECISION,entfac(1),
 224      &     scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
 225      &     WHAM_COMM, IERROR)
 226 #endif
 227           if (me.eq.Master) then
 228 #ifdef DEBUG
 229           write (iout,*) "The FDIMLESS array before sorting"
 230           do i=1,ntot(islice)
 231             write (iout,*) i,fdimless(i)
 232           enddo
 233 #endif
 234 #endif
 235           do i=1,ntot(islice)
 236             iperm(i)=i
 237           enddo
 238           call mysort1(ntot(islice),Fdimless,iperm)
 239 #ifdef DEBUG
 240           write (iout,*) "The FDIMLESS array after sorting"
 241           do i=1,ntot(islice)
 242             write (iout,*) i,iperm(i),fdimless(i)
 243           enddo
 244 #endif
 245           qfree=0.0d0
 246           do i=1,ntot(islice)
 247             qfree=qfree+exp(-fdimless(i)+fdimless(1))
 248           enddo
 249 c          write (iout,*) "qfree",qfree
 250           nlist=1
 251           sumprob=0.0
 252           do i=1,min0(ntot(islice),ensembles)
 253             sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
 254 #ifdef DEBUG
 255             write (iout,*) i,ib,beta_h(ib,iparm),
 256      &       1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
 257      &       potE(iperm(i),iparm),
 258      &       -entfac(iperm(i)),fdimless(i),sumprob
 259 #endif
 260             if (sumprob.gt.0.99d0) goto 122
 261             nlist=nlist+1
 262           enddo
 263   122     continue
 264 #ifdef MPI
 265           endif
 266           call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
 267      &       IERROR)
 268           call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
 269      &       IERROR)
 270           do i=1,nlist
 271             ii=iperm(i)
 272             iproc=0
 273             do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
 274               iproc=iproc+1
 275             enddo
 276             if (iproc.ge.nprocs) then
 277               write (iout,*) "Fatal error: processor out of range",iproc
 278               call flush(iout)
 279               if (bxfile) then
 280                 close (ientout)
 281               else
 282                 close (ientout,status="delete")
 283               endif
 284               return1
 285             endif
 286             ik=ii-indstart(iproc)+1
 287             if (iproc.ne.Master) then
 288               if (me.eq.iproc) then
 289 #ifdef DEBUG
 290                 write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
 291      &           " energy",potE(ik,iparm)
 292 #endif
 293                 call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
 294      &            Master,i,WHAM_COMM,IERROR)
 295               else if (me.eq.Master) then
 296                 call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
 297      &            WHAM_COMM,STATUS,IERROR)
 298               endif
 299             else if (me.eq.Master) then
 300               enepot(i)=potE(ik,iparm)
 301             endif
 302           enddo
 303 #else
 304           do i=1,nlist
 305             enepot(i)=potE(iperm(i),iparm)
 306           enddo
 307 #endif
 308 #ifdef MPI
 309           if (me.eq.Master) then
 310 #endif
 311           write(licz3,'(bz,i3.3)') iparm
 312           write(licz2,'(bz,i2.2)') islice
 313           if (temper.lt.100.0d0) then
 314             write(ctemper,'(f3.0)') temper
 315           else if (temper.lt.1000.0) then
 316             write (ctemper,'(f4.0)') temper
 317           else
 318             write (ctemper,'(f5.0)') temper
 319           endif
 320           if (nparmset.eq.1) then
 321             if (separate_parset) then
 322               write(licz4,'(bz,i3.3)') myparm
 323               pdbname=prefix(:ilen(prefix))//"_par"//licz4
 324             else
 325               pdbname=prefix(:ilen(prefix))
 326             endif
 327           else
 328             pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
 329           endif
 330           if (nslice.eq.1) then
 331             pdbname=pdbname(:ilen(pdbname))//"_T_"//
 332      &        ctemper(:ilen(ctemper))//"pdb"
 333           else
 334             pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
 335      &        ctemper(:ilen(ctemper))//"pdb"
 336           endif
 337           open(ipdb,file=pdbname)
 338           do i=1,nlist
 339             read (ientout,rec=iperm(i))
 340      &        ((csingle(l,k),l=1,3),k=1,nres),
 341      &        ((csingle(l,k+nres),l=1,3),k=nnt,nct),
 342      &        nss,(ihpb(k),jhpb(k),k=1,nss),
 343      &        eini,efree,rmsdev,iscor
 344             do j=1,2*nres
 345               do k=1,3
 346                 c(k,j)=csingle(k,j)
 347               enddo
 348             enddo
 349             eini=fdimless(i)
 350             call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
 351 cc          if (temper.eq.300.0d0) then
 352 cc          write(iout,*) 'Gral i=',iperm(i) ,eini,enepot(i),efree
 353 cc          flush(iout)
 354 cc          endif
 355           enddo
 356 cc          if (temper.eq.300.0d0) then
 357 cc          write(iout,*) 'Gral i=',i ,eini,enepot(i),efree
 358 cc          flush(iout)
 359 cc          endif
 360 #ifdef MPI
 361           endif
 362 #endif
 363         enddo     ! ib
 364       enddo       ! iparm
 365       if (bxfile) then
 366         close(ientout)
 367       else
 368         close(ientout,status="delete")
 369       endif
 370       return
 371       end
 372 !--------------------------------------------------
 373       subroutine mysort1(n, x, ipermut)
 374       implicit none
 375       integer i,j,imax,ipm,n
 376       real x(n)
 377       integer ipermut(n)
 378       real xtemp
 379       do i=1,n
 380         xtemp=x(i)
 381         imax=i
 382         do j=i+1,n
 383           if (x(j).lt.xtemp) then
 384             imax=j
 385             xtemp=x(j)
 386           endif
 387         enddo
 388         x(imax)=x(i)
 389         x(i)=xtemp
 390         ipm=ipermut(imax)
 391         ipermut(imax)=ipermut(i)
 392         ipermut(i)=ipm
 393       enddo
 394       return
 395       end