+
+C Determine the minimum free energies
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
+ do iparm=1,nParmSet
+#ifdef DEBUG
+ write (iout,'(2i5,21f8.2)') i,iparm,
+ & (enetb(k,i,iparm),k=1,21)
+#endif
+ call restore_parm(iparm)
+#ifdef DEBUG
+ write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
+ & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
+ & wtor_d,wsccor,wbond
+#endif
+ do ib=1,nT_h(iparm)
+ if (rescale_mode.eq.1) then
+ quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+ 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)
+ 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)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/
+ & dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+#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)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_mode.eq.0) then
+ do l=1,6
+ fT(l)=1.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ 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 DEBUG
+ write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
+ & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
+ & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr
+#endif
+
+#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
+c write (iout,*) "i",i," ib",ib,
+c & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
+c & " entfac",entfac(i)
+ etot=etot-entfac(i)/beta_h(ib,iparm)
+ if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
+c write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
+ enddo ! ib
+ enddo ! iparm
+ enddo ! i
+#ifdef DEBUG
+ write (iout,*) "The potEmin array before reduction"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#ifdef MPI
+C Determine the minimum energes for all parameter sets and temperatures
+ call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
+ & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=potEmin_t_all(j,i)
+ enddo
+ enddo
+#endif
+ potEmin_min=potEmin_all(1,1)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ if (potEmin_all(j,i).lt.potEmin_min)
+ & potEmin_min=potEmin_all(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "The potEmin array"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#undef DEBUG
+