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