wham with mpich2 (no aliasing) and new dimenstions
[unres.git] / source / cluster / wham / src / probabl.F
index 64e8c50..2f34f61 100644 (file)
@@ -15,6 +15,7 @@
       include "COMMON.SBRIDGE"
       include "COMMON.CHAIN"
       include "COMMON.CLUSTER"
+      include "COMMON.NAMES"
       real*4 csingle(3,maxres2)
       double precision fT(5),fTprim(5),fTbis(5),quot,quotl1,quotl,kfacl,
      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/
      &      ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
      &      eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
       integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon
-      double precision qfree,sumprob,eini,efree,rmsdev
+      integer ires
+      double precision qfree,sumprob,eini,efree,rmsdev,ehomology_constr,
+     &      edfadis,edfator,edfanei,edfabet
       character*80 bxname
       character*2 licz1
       character*5 ctemper
       integer ilen
       external ilen
-      real*4 Fdimless(maxconf)
+      real*4 Fdimless(maxconf),Fdimless_(maxconf)
       double precision energia(0:max_ene)
+      double precision totfree_(maxconf),entfac_(maxconf)
       do i=1,ncon
         list_conf(i)=i
       enddo
 c      do i=1,ncon
 c        write (iout,*) i,list_conf(i)
 c      enddo
+c      do i=1,ncon
+c        write(iout,*) "entrop before", entfac(i),i
+c      enddo
+
 #ifdef MPI
       write (iout,*) me," indstart",indstart(me)," indend",indend(me)
       call daread_ccoords(indstart(me),indend(me))
 #endif
+c      do i=1,ncon
+c        write(iout,*) "entrop after", entfac(i),i
+c      enddo
+
 c      write (iout,*) "ncon",ncon
+
       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)
@@ -71,7 +84,11 @@ c      enddo
      &             dlog(dexp(quotl)+dexp(-quotl))
               enddo
 c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
-              call flush(iout)
+c#ifdef AIX
+c      call flush_(iout)
+c#else
+c      call flush(iout)
+c#endif
             endif
 
 #ifdef MPI
@@ -82,7 +99,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft
         ii=i
 #endif
 c        write (iout,*) "i",i," ii",ii
-        call flush(iout)
+c        call flush(iout)
         if (ib.eq.1) then
           do j=1,nres
             do k=1,3
@@ -103,7 +120,18 @@ c        write (iout,*) "i",i," ii",ii
           call etotal(energia(0),fT)
           totfree(i)=energia(0)
 #ifdef DEBUG
-          write (iout,*) i," energia",(energia(j),j=0,19)
+          write (iout,*) "Conformation",i
+          call enerprint(energia(0),ft)
+      write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+     & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+      do ires=1,nres
+        write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+     &    restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+     &    (c(j,ires+nres),j=1,3)
+      enddo
+
+c          call intout
+c          call flush(iout)
 #endif
           do k=1,max_ene
             enetb(k,i)=energia(k)
@@ -127,6 +155,7 @@ c        write (iout,*) "i",i," ii",ii
         ecorr=enetb(4,i)
         ecorr5=enetb(5,i)
         ecorr6=enetb(6,i)
+cc        if (wcorr6.eq.0) ecorr6=0.0d0
         eel_loc=enetb(7,i)
         eello_turn3=enetb(8,i)
         eello_turn4=enetb(9,i)
@@ -139,47 +168,63 @@ c        write (iout,*) "i",i," ii",ii
         estr=enetb(18,i)
         esccor=enetb(19,i)
         edihcnstr=enetb(20,i)
+        ehomology_constr=enetb(21,i)
+        edfadis=enetb(22,i)
+        edfator=enetb(23,i)
+        edfanei=enetb(24,i)
+        edfabet=enetb(25,i)
 #ifdef SPLITELE
         etot=wsc*evdw+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
+     &  +wstrain*ehpb+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
+     &  +wbond*estr+ehomology_constr+wdfa_dist*edfadis
+     &  +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 #else
         etot=wsc*evdw+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
+     &  +wstrain*ehpb+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
+     &  +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor+wdfa_dist*edfadis
+     &  +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
+     &  +wbond*estr+ehomology_constr
 #endif
+#ifdef MPI
+        Fdimless_(i)=beta_h(ib)*etot+entfac(ii)
+        totfree_(i)=etot
+#else
         Fdimless(i)=beta_h(ib)*etot+entfac(ii)
         totfree(i)=etot
+#endif
 #ifdef DEBUG
-        write (iout,*) i,ii,ib,
+        
+        write (iout,*) "etrop", i,ii,ib,
      &   1.0d0/(1.987d-3*beta_h(ib)),totfree(i),
      &   entfac(ii),Fdimless(i)
 #endif
       enddo   ! i
 #ifdef MPI
-      call MPI_Gatherv(Fdimless(1),scount(me),
+      call MPI_Gatherv(Fdimless_(1),scount(me),
      & MPI_REAL,Fdimless(1),
      & scount(0),idispl(0),MPI_REAL,Master,
      & MPI_COMM_WORLD, IERROR)
-      call MPI_Gatherv(totfree(1),scount(me),
+      call MPI_Gatherv(totfree_(1),scount(me),
      & MPI_DOUBLE_PRECISION,totfree(1),
      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
      & MPI_COMM_WORLD, IERROR)
       call MPI_Gatherv(entfac(indstart(me)+1),scount(me),
-     & MPI_DOUBLE_PRECISION,entfac(1),
+     & MPI_DOUBLE_PRECISION,entfac_(1),
      & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,
      & MPI_COMM_WORLD, IERROR)
       if (me.eq.Master) then
+        do i=1,ncon
+          entfac(i)=entfac_(i)
+        enddo
 #endif
 #ifdef DEBUG
         write (iout,*) "The FDIMLESS array before sorting"
@@ -199,19 +244,21 @@ c        write (iout,*) "i",i," ii",ii
         enddo
         qfree=0.0d0
         do i=1,ncon
-          qfree=qfree+exp(-fdimless(i)+fdimless(1))
+          qfree=qfree+dexp(dble(-fdimless(i)+fdimless(1)))
         enddo
 c        write (iout,*) "qfree",qfree
         nlist=1
         sumprob=0.0
         do i=1,min0(ncon,maxstr_proc)-1 
-          sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
+          sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree 
+c#define DEBUG
 #ifdef DEBUG
-          write (iout,*) i,ib,beta_h(ib),
+          write (iout,*) 'i=',i,ib,beta_h(ib),
      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
      &     totfree(list_conf(i)),
      &     -entfac(list_conf(i)),fdimless(i),sumprob
 #endif
+c#undef DEBUG
           if (sumprob.gt.prob_limit) goto 122
 c          if (sumprob.gt.1.00d0) goto 122
           nlist=nlist+1