--- /dev/null
+ subroutine make_ensembles(islice,*)
+! construct the conformational ensembles at REMD temperatures
+ implicit none
+ include "DIMENSIONS"
+ include "DIMENSIONS.ZSCOPT"
+ include "DIMENSIONS.FREE"
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.MPI"
+ integer ierror,errcode,status(MPI_STATUS_SIZE)
+#endif
+ include "COMMON.IOUNITS"
+ include "COMMON.CONTROL"
+ include "COMMON.FREE"
+ include "COMMON.ENERGIES"
+ include "COMMON.FFIELD"
+ include "COMMON.INTERACT"
+ include "COMMON.SBRIDGE"
+ include "COMMON.CHAIN"
+ include "COMMON.PROTFILES"
+ include "COMMON.PROT"
+ real*4 csingle(3,maxres2)
+ double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
+ & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
+ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+ & escloc,
+ & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt
+ integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist
+ double precision qfree,sumprob,eini,efree,rmsdev
+ character*80 bxname
+ character*2 licz1,licz2
+ character*3 licz3,licz4
+ character*5 ctemper
+ integer ilen
+ external ilen
+ real*4 Fdimless(MaxStr)
+ double precision enepot(MaxStr)
+ integer iperm(MaxStr)
+ integer islice
+
+#ifdef MPI
+ if (me.eq.Master) then
+#endif
+ write (licz2,'(bz,i2.2)') islice
+ if (nslice.eq.1) then
+ if (.not.separate_parset) then
+ bxname = prefix(:ilen(prefix))//".bx"
+ else
+ write (licz3,'(bz,i3.3)') myparm
+ bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx"
+ endif
+ else
+ if (.not.separate_parset) then
+ bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx"
+ else
+ write (licz3,'(bz,i3.3)') myparm
+ bxname = prefix(:ilen(prefix))//"par_"//licz3//
+ & "_slice_"//licz2//".bx"
+ endif
+ endif
+ open (ientout,file=bxname,status="unknown",
+ & form="unformatted",access="direct",recl=lenrec1)
+#ifdef MPI
+ endif
+#endif
+ do iparm=1,nParmSet
+ if (iparm.ne.iparmprint) exit
+ call restore_parm(iparm)
+ do ib=1,nT_h(iparm)
+#ifdef DEBUG
+ write (iout,*) "iparm",iparm," ib",ib
+#endif
+ temper=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+c quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+c quotl=1.0d0
+c kfacl=1.0d0
+c do l=1,5
+c quotl1=quotl
+c quotl=quotl*quot
+c kfacl=kfacl*kfac
+c fT(l)=kfacl/(kfacl-1.0d0+quotl)
+c enddo
+ if (rescale_mode.eq.1) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+#if defined(FUNCTH)
+ tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=quot
+#else
+ ft(6)=1.0d0
+#endif
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ fT(l)=kfacl/(kfacl-1.0d0+quotl)
+ enddo
+ else if (rescale_mode.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+#if defined(FUNCTH)
+ tt=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0
+#elif defined(FUNCT)
+ ft(6)=quot
+#else
+ ft(6)=1.0d0
+#endif
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/
+ & dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_mode.eq.0) then
+ do l=1,5
+ fT(l)=0.0d0
+ enddo
+ else
+ write (iout,*)
+ & "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode
+ call flush(iout)
+ return1
+ endif
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+ evdw=enetb(1,i,iparm)
+ evdw_t=enetb(21,i,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,i,iparm)
+ evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,i,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,i,iparm)
+ evdw1=enetb(16,i,iparm)
+#else
+ ees=enetb(3,i,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,i,iparm)
+ ecorr5=enetb(5,i,iparm)
+ ecorr6=enetb(6,i,iparm)
+ eel_loc=enetb(7,i,iparm)
+ eello_turn3=enetb(8,i,iparm)
+ eello_turn4=enetb(9,i,iparm)
+ eturn6=enetb(10,i,iparm)
+ ebe=enetb(11,i,iparm)
+ escloc=enetb(12,i,iparm)
+ etors=enetb(13,i,iparm)
+ etors_d=enetb(14,i,iparm)
+ ehpb=enetb(15,i,iparm)
+ estr=enetb(18,i,iparm)
+ esccor=enetb(19,i,iparm)
+ edihcnstr=enetb(20,i,iparm)
+#ifdef SPLITELE
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#endif
+#ifdef MPI
+ Fdimless(i)=
+ & beta_h(ib,iparm)*etot-entfac(i)
+ potE(i,iparm)=etot
+#ifdef DEBUG
+ write (iout,*) i,indstart(me)+i-1,ib,
+ & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm),
+ & -entfac(i),Fdimless(i)
+#endif
+#else
+ Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i)
+ potE(i,iparm)=etot
+#endif
+ enddo ! i
+#ifdef MPI
+ call MPI_Gatherv(Fdimless(1),scount(me),
+ & MPI_REAL,Fdimless(1),
+ & scount(0),idispl(0),MPI_REAL,Master,
+ & WHAM_COMM, IERROR)
+#ifdef DEBUG
+ call MPI_Gatherv(potE(1,iparm),scount(me),
+ & MPI_DOUBLE_PRECISION,potE(1,iparm),
+ & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
+ & WHAM_COMM, IERROR)
+ call MPI_Gatherv(entfac(1),scount(me),
+ & MPI_DOUBLE_PRECISION,entfac(1),
+ & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
+ & WHAM_COMM, IERROR)
+#endif
+ if (me.eq.Master) then
+#ifdef DEBUG
+ write (iout,*) "The FDIMLESS array before sorting"
+ do i=1,ntot(islice)
+ write (iout,*) i,fdimless(i)
+ enddo
+#endif
+#endif
+ do i=1,ntot(islice)
+ iperm(i)=i
+ enddo
+ call mysort1(ntot(islice),Fdimless,iperm)
+#ifdef DEBUG
+ write (iout,*) "The FDIMLESS array after sorting"
+ do i=1,ntot(islice)
+ write (iout,*) i,iperm(i),fdimless(i)
+ enddo
+#endif
+ qfree=0.0d0
+ do i=1,ntot(islice)
+ qfree=qfree+exp(-fdimless(i)+fdimless(1))
+ enddo
+c write (iout,*) "qfree",qfree
+ nlist=1
+ sumprob=0.0
+ do i=1,min0(ntot(islice),ensembles)
+ sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree
+#ifdef DEBUG
+ write (iout,*) i,ib,beta_h(ib,iparm),
+ & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i),
+ & potE(iperm(i),iparm),
+ & -entfac(iperm(i)),fdimless(i),sumprob
+#endif
+ if (sumprob.gt.0.99d0) goto 122
+ nlist=nlist+1
+ enddo
+ 122 continue
+#ifdef MPI
+ endif
+ call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM,
+ & IERROR)
+ call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM,
+ & IERROR)
+ do i=1,nlist
+ ii=iperm(i)
+ iproc=0
+ do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc))
+ iproc=iproc+1
+ enddo
+ if (iproc.ge.nprocs) then
+ write (iout,*) "Fatal error: processor out of range",iproc
+ call flush(iout)
+ if (bxfile) then
+ close (ientout)
+ else
+ close (ientout,status="delete")
+ endif
+ return1
+ endif
+ ik=ii-indstart(iproc)+1
+ if (iproc.ne.Master) then
+ if (me.eq.iproc) then
+#ifdef DEBUG
+ write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik,
+ & " energy",potE(ik,iparm)
+#endif
+ call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION,
+ & Master,i,WHAM_COMM,IERROR)
+ else if (me.eq.Master) then
+ call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i,
+ & WHAM_COMM,STATUS,IERROR)
+ endif
+ else if (me.eq.Master) then
+ enepot(i)=potE(ik,iparm)
+ endif
+ enddo
+#else
+ do i=1,nlist
+ enepot(i)=potE(iperm(i),iparm)
+ enddo
+#endif
+#ifdef MPI
+ if (me.eq.Master) then
+#endif
+ write(licz3,'(bz,i3.3)') iparm
+ write(licz2,'(bz,i2.2)') islice
+ if (temper.lt.100.0d0) then
+ write(ctemper,'(f3.0)') temper
+ else if (temper.lt.1000.0) then
+ write (ctemper,'(f4.0)') temper
+ else
+ write (ctemper,'(f5.0)') temper
+ endif
+ if (nparmset.eq.1) then
+ if (separate_parset) then
+ write(licz4,'(bz,i3.3)') myparm
+ pdbname=prefix(:ilen(prefix))//"_par"//licz4
+ else
+ pdbname=prefix(:ilen(prefix))
+ endif
+ else
+ pdbname=prefix(:ilen(prefix))//"_parm_"//licz3
+ endif
+ if (nslice.eq.1) then
+ pdbname=pdbname(:ilen(pdbname))//"_T_"//
+ & ctemper(:ilen(ctemper))//"pdb"
+ else
+ pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"//
+ & ctemper(:ilen(ctemper))//"pdb"
+ endif
+ open(ipdb,file=pdbname)
+ do i=1,nlist
+ read (ientout,rec=iperm(i))
+ & ((csingle(l,k),l=1,3),k=1,nres),
+ & ((csingle(l,k+nres),l=1,3),k=nnt,nct),
+ & nss,(ihpb(k),jhpb(k),k=1,nss),
+ & eini,efree,rmsdev,iscor
+ do j=1,2*nres
+ do k=1,3
+ c(k,j)=csingle(k,j)
+ enddo
+ enddo
+ eini=fdimless(i)
+ call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ enddo ! ib
+ enddo ! iparm
+ if (bxfile) then
+ close(ientout)
+ else
+ close(ientout,status="delete")
+ endif
+ return
+ end
+!--------------------------------------------------
+ subroutine mysort1(n, x, ipermut)
+ implicit none
+ integer i,j,imax,ipm,n
+ real x(n)
+ integer ipermut(n)
+ real xtemp
+ do i=1,n
+ xtemp=x(i)
+ imax=i
+ do j=i+1,n
+ if (x(j).lt.xtemp) then
+ imax=j
+ xtemp=x(j)
+ endif
+ enddo
+ x(imax)=x(i)
+ x(i)=xtemp
+ ipm=ipermut(imax)
+ ipermut(imax)=ipermut(i)
+ ipermut(i)=ipm
+ enddo
+ return
+ end
projects / unres.git / blobdiff