X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fmake_ensemble1.F;fp=source%2Fwham%2Fsrc-NEWSC%2Fmake_ensemble1.F;h=5d7b750867c2f18493eb1fdd68c0875ba4c38b30;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/wham/src-NEWSC/make_ensemble1.F b/source/wham/src-NEWSC/make_ensemble1.F
new file mode 100755
index 0000000..5d7b750
--- /dev/null
+++ b/source/wham/src-NEWSC/make_ensemble1.F
@@ -0,0 +1,375 @@
+      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