cluster average over multiple temperatures

diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL
index bb11a5e..6894685 100644 (file)
--- a/source/cluster/wham/src-M/COMMON.CONTROL
+++ b/source/cluster/wham/src-M/COMMON.CONTROL
@@ -13,14 +13,15 @@
      & with_dihed_constr,with_theta_constr,out1file,
      & print_homology_restraints,
      & print_contact_map,print_homology_models,
-     & read2sigma,l_homo,read_homol_frag,dist1cut
+     & read2sigma,l_homo,read_homol_frag,dist1cut,cumul_prob
       common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2,
      & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
      & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr,
      & efree,iopt,nstart,nend,symetr,
      & tor_mode,shield_mode,
      & constr_dist,out1file,
-     & constr_homology,homol_nset,read2sigma,read_homol_frag,dist1cut
+     & constr_homology,homol_nset,read2sigma,read_homol_frag,dist1cut,
+     & cumul_prob
       common /homol/  waga_homology(10),
      & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut,
      & iset,ihset,l_homo(max_template,maxdim),

diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F
index 3c7e430..5b3cfe8 100644 (file)
--- a/source/cluster/wham/src-M/main_clust.F
+++ b/source/cluster/wham/src-M/main_clust.F
@@ -41,6 +41,7 @@ C
       double precision hrtime,mintime,sectime
       logical eof
       external ilen
+      integer nTend
 #ifdef MPI
       call MPI_Init( IERROR )
       call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR )
@@ -123,7 +124,12 @@ c      call flush(iout)
       write (iout,*) 'from read_coords: ncon',ncon
       
       write (iout,*) "nT",nT
-      do iT=1,nT
+      if (cumul_prob) then
+        nTend=1
+      else
+        nTend=nT
+      endif
+      do iT=1,nTend
       write (iout,*) "iT",iT
 #ifdef MPI
       call work_partition(.true.,ncon)

diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F
index 22db88f..198c5f5 100644 (file)
--- a/source/cluster/wham/src-M/probabl.F
+++ b/source/cluster/wham/src-M/probabl.F
@@ -1,4 +1,4 @@
-      subroutine probabl(ib,nlist,ncon,*)
+      subroutine probabl(ibb,nlist,ncon,*)
 ! construct the conformational ensembles at REMD temperatures
       implicit none
       include "DIMENSIONS"
@@ -8,6 +8,7 @@
       include "COMMON.MPI"
       integer ierror,errcode,status(MPI_STATUS_SIZE) 
 #endif
+      include "COMMON.CONTROL"
       include "COMMON.IOUNITS"
       include "COMMON.FREE"
       include "COMMON.FFIELD"
@@ -21,9 +22,11 @@
       double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc,
      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
-     &      evdw_t,esaxs
-      integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
-      double precision qfree,sumprob,eini,efree,rmsdev
+     &      evdw_t,esaxs,eliptran,ehomology_constr
+      integer i,ii,ik,iproc,iscor,j,k,l,ib,ibb,ib_start,ib_end,nlist,
+     &  ncon
+      double precision qfree,sumprob,eini,rmsdev
+      real*4 probab(maxconf),totprob(maxconf),aux
       character*80 bxname
       character*2 licz1
       character*5 ctemper
@@ -32,82 +35,86 @@
       real*4 Fdimless(maxconf), Fdimless_buf(maxconf)
       double precision energia(0:max_ene), totfree_buf(0:maxconf),
      &  entfac_buf(maxconf)
+      if (cumul_prob) then
+        write (iout,*) 
+     &  "Probabilities will be summed over all temperatures"
+        ib_start=1
+        ib_end=nT
+      else
+        ib_start=ibb
+        ib_end=ibb
+      endif
+      probab=0.0
+      totprob=0.0
+      Fdimless=0.0d0
+      Fdimless_buf=0.0d0
+      totfree=0.0d0
+      totfree_buf=0.0d0
       do i=1,ncon
         list_conf(i)=i
       enddo
-c      do i=1,ncon
-c        write (iout,*) i,list_conf(i)
-c      enddo
-#ifdef MPI
-      write (iout,*) me," indstart",indstart(me)," indend",indend(me)
-      call daread_ccoords(indstart(me),indend(me))
-#endif
-C      write (iout,*) "ncon",ncon
-C      call flush(iout)
-      temper=1.0d0/(beta_h(ib)*1.987D-3)
-c      write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper
-c      quot=1.0d0/(T0*beta_h(ib)*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
-C#define DEBUG
-            if (rescale_mode.eq.1) then
-              quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
-              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
+      ib=1
+      if (rescale_mode.eq.1) then
+        quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
+        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
 #if defined(FUNCTH)
-              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
-     &                  320.0d0
-              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
-             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
-     &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+        ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &                320.0d0
+        ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+        ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+     &            /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 #elif defined(FUNCT)
-              fT(6)=betaT/T0
-              ftprim(6)=1.0d0/T0
-              ftbis(6)=0.0d0
+        fT(6)=betaT/T0
+        ftprim(6)=1.0d0/T0
+        ftbis(6)=0.0d0
 #else
-              fT(6)=1.0d0
-              ftprim(6)=0.0d0
-              ftbis(6)=0.0d0
+        fT(6)=1.0d0
+        ftprim(6)=0.0d0
+        ftbis(6)=0.0d0
 #endif
 
-            else if (rescale_mode.eq.2) then
-              quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
-              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)*1.987D-3),ft
-c              call flush(iout)
+      else if (rescale_mode.eq.2) then
+        quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
+        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)*1.987D-3),ft
+c            call flush(iout)
 #if defined(FUNCTH)
-              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
-     &                  320.0d0
-              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
-             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+            ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &                320.0d0
+            ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+            ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 #elif defined(FUNCT)
-              fT(6)=betaT/T0
-              ftprim(6)=1.0d0/T0
-              ftbis(6)=0.0d0
+            fT(6)=betaT/T0
+            ftprim(6)=1.0d0/T0
+            ftbis(6)=0.0d0
 #else
-              fT(6)=1.0d0
-              ftprim(6)=0.0d0
-              ftbis(6)=0.0d0
+            fT(6)=1.0d0
+            ftprim(6)=0.0d0
+            ftbis(6)=0.0d0
+#endif
+          endif
+c      do i=1,ncon
+c        write (iout,*) i,list_conf(i)
+c      enddo
+#ifdef MPI
+      write (iout,*) me," indstart",indstart(me)," indend",indend(me)
+      call daread_ccoords(indstart(me),indend(me))
 #endif
-            endif
+C      write (iout,*) "ncon",ncon
+C      call flush(iout)
 
 #ifdef MPI
       do i=1,scount(me)
@@ -116,151 +123,263 @@ c              call flush(iout)
       do i=1,ncon
         ii=i
 #endif
-C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
-c        call flush(iout)
-        if (ib.eq.1) then
-          do j=1,nres
-            do k=1,3
-              c(k,j)=allcart(k,j,i)
-              c(k,j+nres)=allcart(k,j+nres,i)
-C              write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
-C     &        c(k,j+nres),allcart(k,j+nres,i)
-            enddo
-          enddo
-C          write(iout,*) "out of j loop"
-C          call flush(iout)
+        do j=1,nres
           do k=1,3
-            c(k,nres+1)=c(k,1)
-            c(k,nres+nres)=c(k,nres)
+            c(k,j)=allcart(k,j,i)
+            c(k,j+nres)=allcart(k,j+nres,i)
+C            write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j),
+C     &      c(k,j+nres),allcart(k,j+nres,i)
           enddo
-C          write(iout,*) "after nres+nres",nss_all(i)
-C          call flush(iout)
-          nss=nss_all(i)
-          do j=1,nss
-            ihpb(j)=ihpb_all(j,i)
-            jhpb(j)=jhpb_all(j,i)
-          enddo 
-          call int_from_cart1(.false.)
-C          write(iout,*) "before etotal"
-C          call flush(iout)
-          call etotal(energia(0),fT)
-          totfree(i)=energia(0)         
-          totfree_buf(i)=totfree(i)
+        enddo
+C        write(iout,*) "out of j loop"
+C        call flush(iout)
+        do k=1,3
+          c(k,nres+1)=c(k,1)
+          c(k,nres+nres)=c(k,nres)
+        enddo
+C        write(iout,*) "after nres+nres",nss_all(i)
+C        call flush(iout)
+        nss=nss_all(i)
+        do j=1,nss
+          ihpb(j)=ihpb_all(j,i)
+          jhpb(j)=jhpb_all(j,i)
+        enddo 
+        call int_from_cart1(.false.)
+C        write(iout,*) "before etotal"
+C        call flush(iout)
+        call etotal(energia(0),fT)
 c          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
 c          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
 c          call enerprint(energia(0),fT)
 c          call pdbout(totfree(i),16,i)
 #ifdef DEBUG
-          write (iout,*) i," energia",(energia(j),j=0,max_ene)
-          write (iout,*) "etot", etot
-          write (iout,*) "ft(6)", ft(6)
+        write (iout,*) i," energia",(energia(j),j=0,max_ene)
+        write (iout,*) "etot", etot
+        write (iout,*) "ft(6)", ft(6)
 #endif
-          do k=1,max_ene
-            enetb(k,i)=energia(k)
-          enddo
-        endif
-        evdw=enetb(1,i)
-c        write (iout,*) evdw
-       etot=energia(0)
+        do k=1,max_ene
+          enetb(k,i)=energia(k)
+        enddo
+      enddo ! i
+c#endif
+C        write (iout,*) "i",i," ii",ii,"ib",ib,scount(me)
+c        call flush(iout)
+      do ib=ib_start,ib_end
+          temper=1.0d0/(beta_h(ib)*1.987D-3)
+#ifdef MPI
+        do i=1,scount(me)
+          ii=i+indstart(me)-1
+#else
+        do i=1,ncon
+          ii=i
+#endif
+          evdw=enetb(1,i)
+c          write (iout,*) evdw
+          etot=energia(0)
 #ifdef SCP14
-        evdw2_14=enetb(17,i)
-        evdw2=enetb(2,i)+evdw2_14
+          evdw2_14=enetb(17,i)
+          evdw2=enetb(2,i)+evdw2_14
+#else
+          evdw2=enetb(2,i)
+          evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+          ees=enetb(3,i)
+          evdw1=enetb(16,i)
+#else
+          ees=enetb(3,i)
+          evdw1=0.0d0
+#endif
+          ecorr=enetb(4,i)
+          ecorr5=enetb(5,i)
+          ecorr6=enetb(6,i)
+          eel_loc=enetb(7,i)
+          eello_turn3=enetb(8,i)
+          eello_turn4=enetb(9,i)
+          eturn6=enetb(10,i)
+          ebe=enetb(11,i)
+          escloc=enetb(12,i)
+          etors=enetb(13,i)
+          etors_d=enetb(14,i)
+          ehpb=enetb(15,i)
+          estr=enetb(18,i)
+          esccor=enetb(19,i)
+          edihcnstr=enetb(20,i)
+          evdw_t=enetb(21,i)
+          eliptran=enetb(22,i)
+          ehomology_constr=enetb(24,i)
+          esaxs=enetb(25,i)
+          if (rescale_mode.eq.1) then
+            quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
+            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
+#if defined(FUNCTH)
+            ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &                320.0d0
+            ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+            ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+     &            /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+            fT(6)=betaT/T0
+            ftprim(6)=1.0d0/T0
+            ftbis(6)=0.0d0
 #else
-        evdw2=enetb(2,i)
-        evdw2_14=0.0d0
+            fT(6)=1.0d0
+            ftprim(6)=0.0d0
+            ftbis(6)=0.0d0
 #endif
+
+          else if (rescale_mode.eq.2) then
+            quot=1.0d0/(T0*beta_h(ib)*1.987D-3)
+            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)*1.987D-3),ft
+c            call flush(iout)
+#if defined(FUNCTH)
+            ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &                320.0d0
+            ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+            ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+     &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+            fT(6)=betaT/T0
+            ftprim(6)=1.0d0/T0
+            ftbis(6)=0.0d0
+#else
+            fT(6)=1.0d0
+            ftprim(6)=0.0d0
+            ftbis(6)=0.0d0
+#endif
+          endif
+c          write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
+c          call flush(iout)
+c          call etotal(energia(0),fT)
+c#ifdef CHUJ
 #ifdef SPLITELE
-        ees=enetb(3,i)
-        evdw1=enetb(16,i)
+          if (shield_mode.gt.0) then
+          etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &     +welec*ft(1)*ees
+     &     +ft(1)*wvdwpp*evdw1
+     &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
+     &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
+     &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
+     &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
+     &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
+     &     +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr
+     &     +wliptran*eliptran+wsaxs*esaxs
+          else
+          etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+welec*ft(1)*ees
+     &     +wvdwpp*evdw1
+     &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
+     &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
+     &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
+     &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
+     &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
+     &     +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr
+     &     +wliptran*eliptran+wsaxs*esaxs
+          endif
 #else
-        ees=enetb(3,i)
-        evdw1=0.0d0
+          if (shield_mode.gt.0) then
+          etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &     +welec*ft(1)*(ees+evdw1)
+     &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
+     &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
+     &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
+     &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
+     &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
+     &     +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr
+     &     +wliptran*eliptran+wsaxs*esaxs
+          else
+          etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+     &     +welec*ft(1)*(ees+evdw1)
+     &     +wang*ebe+wtor*ft(1)*etors+wscloc*escloc
+     &     +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5
+     &     +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4
+     &     +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6
+     &     +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d
+     &     +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr
+     &     +wliptran*eliptran+wsaxs*esaxs
+          endif
 #endif
-        ecorr=enetb(4,i)
-        ecorr5=enetb(5,i)
-        ecorr6=enetb(6,i)
-        eel_loc=enetb(7,i)
-        eello_turn3=enetb(8,i)
-        eello_turn4=enetb(9,i)
-        eturn6=enetb(10,i)
-        ebe=enetb(11,i)
-        escloc=enetb(12,i)
-        etors=enetb(13,i)
-        etors_d=enetb(14,i)
-        ehpb=enetb(15,i)
-        estr=enetb(18,i)
-        esccor=enetb(19,i)
-        edihcnstr=enetb(20,i)
-        esaxs=enetb(25,i)
-c#ifdef SPLITELE
-c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+
-c     &ft(1)*welec*ees+wvdwpp*evdw1
-c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
-c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
-c     &  +ft(2)*wturn3*eello_turn3
-c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
-c     &  +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-c     &  +wbond*estr
+#ifdef DEBUG
+          write (iout,*) "etot2", etot
+          write (iout,*) "evdw", wsc, evdw,evdw_t
+          write (iout,*) "evdw2", wscp, evdw2
+          write (iout,*) "welec", ft(1),welec,ees
+          write (iout,*) "evdw1", wvdwpp,evdw1
+          write (iout,*) "ebe", ebe,wang
+#endif        
 c#else
-c        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1)
-c     &  +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-c     &  +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
-c     &  +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
-c     &  +ft(2)*wturn3*eello_turn3
-c     &  +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
-c     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-c     &  +wbond*estr
+c          etot=energia(0)
 c#endif
 #ifdef DEBUG
-        write (iout,*) "etot2", etot
-        write (iout,*) "evdw", wsc, evdw,evdw_t
-        write (iout,*) "evdw2", wscp, evdw2
-        write (iout,*) "welec", ft(1),welec,ees
-        write (iout,*) "evdw1", wvdwpp,evdw1
-        write (iout,*) "ebe", ebe,wang
-#endif        
-        Fdimless(i)=beta_h(ib)*etot+entfac(ii)
-        Fdimless_buf(i)=Fdimless(i)
-        totfree(i)=etot
-        totfree_buf(i)=totfree(i)
+          write (iout,'(2i5,30f10.2)') i,ii,etot,
+     &    (enetb(ijk,i),ijk=1,max_ene)
+#endif
+          Fdimless(i)=beta_h(ib)*etot+entfac(ii)
 #ifdef DEBUG
-        write (iout,*) "fdim calc", i,ii,ib,
-     &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
-     &   entfac(ii),Fdimless(i)
+          write (iout,*) "fdim calc", i,ii,ib,
+     &     1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
+     &    entfac(ii),Fdimless(i)
 #endif
-      enddo   ! i
-
-      do ijk=1,maxconf
-      entfac_buf(ijk)=entfac(ijk)
-      Fdimless_buf(ijk)=Fdimless(ijk)
-      enddo
-      do ijk=0,maxconf
-      totfree_buf(ijk)=totfree(ijk)
-      enddo
-
-
-c      scount_buf=scount(me)
-c      scount_buf2=scount(0)
-
-c      entfac_buf(indstart(me)+1)=entfac(indstart(me)+1)
+          Fdimless_buf(i)=Fdimless(i)
+          totfree_buf(i)=totfree(i)
+        enddo   ! i
 
 #ifdef MPI
 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)"
-      call MPI_Gatherv(Fdimless_buf(1),scount(me),
-     & MPI_REAL,Fdimless(1),
-     & scount(0),idispl(0),MPI_REAL,Master,
-     & MPI_COMM_WORLD, IERROR)
+        call MPI_Gatherv(Fdimless_buf(1),scount(me),
+     &    MPI_REAL,Fdimless(1),
+     &    scount(0),idispl(0),MPI_REAL,Master,
+     &    MPI_COMM_WORLD, IERROR)
 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)"
-      call MPI_Gatherv(totfree_buf(1),scount(me),
-     & MPI_DOUBLE_PRECISION,totfree(1),
-     & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
-     & MPI_COMM_WORLD, IERROR)
+        call MPI_Gatherv(totfree_buf(1),scount(me),
+     &    MPI_DOUBLE_PRECISION,totfree(1),
+     &    scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
+     &    MPI_COMM_WORLD, IERROR)
 c      WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)"
+c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
+        aux=Fdimless(1)
+        do i=1,ncon
+          if (Fdimless(i).lt.aux) aux=Fdimless(i)
+        enddo
+        sumprob=0.0d0
+        do i=1,ncon
+          probab(i)=exp(aux-fdimless(i))
+          sumprob=sumprob+probab(i)
+        enddo
+        do i=1,ncon
+          probab(i)=probab(i)/sumprob
+          totprob(i)=totprob(i)+probab(i)
+c          write (iout,*) ib,i,probab(i),totprob(i)
+        enddo
+      enddo ! ib
+      sumprob=0.0d0
+      do i=1,ncon
+        sumprob=sumprob+totprob(i)
+      enddo
+      write (iout,*) "sumprob",sumprob
+      do i=1,ncon
+        totprob(i)=totprob(i)/sumprob
+      enddo
+      do i=1,ncon
+c        Fdimless(i)=-dlog(Fdimless(i))
+        Fdimless(i)=-alog(totprob(i))
+      enddo
       call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me),
-     & MPI_DOUBLE_PRECISION,entfac(1),
-     & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
-     & MPI_COMM_WORLD, IERROR)
+     &    MPI_DOUBLE_PRECISION,entfac(1),
+     &    scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
+     &    MPI_COMM_WORLD, IERROR)
 c      WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)"
       if (me.eq.Master) then
 c      WRITE (iout,*) "me.eq.Master"
@@ -319,9 +438,16 @@ c     &   " indend",indend(iproc)
 c      enddo
       write (iout,*) "nlist",nlist
 #endif
-      write (iout,'(80(1h-)/i4,a,f6.3,a,f5.1/80(1h-))') 
+      if (cumul_prob) then
+      write (iout,'(80(1h-)/i4,a,f6.3,a,10f6.1/80(1h-))') 
      & nlist," conformations found within",prob_limit,
-     & " probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ib))
+     & " probablity cut_off in the temperature range ",
+     &  (1.0d0/(1.987d-3*beta_h(ib)),ib=1,nT)
+      else
+      write (iout,'(80(1h-)/i4,a,f6.3,a,f6.1/80(1h-))') 
+     & nlist," conformations found within",prob_limit,
+     &" probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ibb))
+      endif
       return
       end
 !--------------------------------------------------

diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F
index 05d6f22..5a42f4b 100644 (file)
--- a/source/cluster/wham/src-M/readrtns.F
+++ b/source/cluster/wham/src-M/readrtns.F
@@ -110,7 +110,9 @@ C long axis of side chain
       lside = index(controlcard,"SIDE").gt.0
       efree = index(controlcard,"EFREE").gt.0
       call readi(controlcard,'NTEMP',nT,1)
-      write (iout,*) "nT",nT
+      cumul_prob=nt.lt.0
+      nT=iabs(nT)
+      write (iout,*) "nT",nT," cumul_prob",cumul_prob
       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
       write (iout,*) "nT",nT
       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)