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