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