source/wham/src-M/wham_calc1.F

   1               subroutine WHAM_CALC(islice,*)
   2 ! Weighed Histogram Analysis Method (WHAM) code
   3 ! Written by A. Liwo based on the work of Kumar et al.,
   4 ! J.Comput.Chem., 13, 1011 (1992)
   5 !
   6 ! 2/1/05 Multiple temperatures allowed.
   7 ! 2/2/05 Free energies calculated directly from data points
   8 !  acc. to Eq. (21) of Kumar et al.; final histograms also
   9 !  constructed based on this equation.
  10 ! 2/12/05 Multiple parameter sets included
  11 !
  12 ! 2/2/05 Parallel version
  13       implicit none
  14       include "DIMENSIONS"
  15       include "DIMENSIONS.ZSCOPT"
  16       include "DIMENSIONS.FREE"
  17       integer nGridT
  18       parameter (NGridT=400)
  19       integer MaxBinRms,MaxBinRgy
  20       parameter (MaxBinRms=100,MaxBinRgy=100)
  21       integer MaxHdim
  22 c      parameter (MaxHdim=200000)
  23       parameter (MaxHdim=200)
  24       integer maxinde
  25       parameter (maxinde=200)
  26 #ifdef MPI
  27       include "mpif.h"
  28       include "COMMON.MPI"
  29       integer ierror,errcode,status(MPI_STATUS_SIZE)
  30 #endif
  31       include "COMMON.CONTROL"
  32       include "COMMON.IOUNITS"
  33       include "COMMON.FREE"
  34       include "COMMON.ENERGIES"
  35       include "COMMON.FFIELD"
  36       include "COMMON.SBRIDGE"
  37       include "COMMON.PROT"
  38       include "COMMON.ENEPS"
  39       include "COMMON.SHIELD"
  40       integer MaxPoint,MaxPointProc
  41       parameter (MaxPoint=MaxStr,
  42      & MaxPointProc=MaxStr_Proc)
  43       double precision finorm_max,potfac,entmin,entmax,expfac,vf
  44       parameter (finorm_max=1.0d0)
  45       integer islice
  46       integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  47       integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
  48      &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
  49       integer htot(0:MaxHdim),histent(0:2000)
  50       double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  51       double precision energia(0:max_ene)
  52 #ifdef MPI
  53       integer tmax_t,upindE_p
  54       double precision fi_p(MaxR,MaxT_h,Max_Parm),
  55      &  fi_p_min(MaxR,MaxT_h,Max_Parm)
  56       double precision sumW_p(0:nGridT,Max_Parm),
  57      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
  58      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
  59      & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
  60      & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
  61      & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
  62      & sumEbis_p(0:nGridT,Max_Parm)
  63       double precision hfin_p(0:MaxHdim,maxT_h),
  64      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
  65      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  66       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  67       double precision potEmin_t,entmin_p,entmax_p,
  68      & potEmin_t_all(maxT_h,Max_Parm)
  69       integer histent_p(0:2000)
  70       logical lprint /.true./
  71 #endif
  72       double precision delta_T /1.0d0/
  73       double precision rgymin,rmsmin,rgymax,rmsmax
  74       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
  75      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
  76      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
  77      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
  78      & weight,econstr
  79       double precision fi(MaxR,maxT_h,Max_Parm),
  80      & fi_min(MaxR,maxT_h,Max_Parm),
  81      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
  82      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
  83      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
  84      & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min,
  85      & hfin_ent(0:MaxHdim),vmax,aux,entfac_min
  86       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  87      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
  88      &  eplus,eminus,logfac,tanhT,tt
  89       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  90      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  91      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  92      &  eliptran,etube
  93
  94       integer ind_point(maxpoint),upindE,indE
  95       character*16 plik
  96       character*1 licz1
  97       character*2 licz2
  98       character*3 licz3
  99       character*128 nazwa
 100       integer ilen
 101       external ilen
 102
 103       write(licz2,'(bz,i2.2)') islice
 104       nbin1 = 1.0d0/delta
 105       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
 106      &  i2/80(1h-)//)') islice
 107       write (iout,*) "delta",delta," nbin1",nbin1
 108       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 109       call flush(iout)
 110       dmin=0.0d0
 111       tmax=0
 112       potEmin=1.0d10
 113        do i=1,nParmset
 114          do j=1,nT_h(i)
 115            potEmin_all(j,i)=1.0d10
 116          enddo
 117        enddo
 118       rgymin=1.0d10
 119       rmsmin=1.0d10
 120       rgymax=0.0d0
 121       rmsmax=0.0d0
 122       do t=0,MaxN
 123         htot(t)=0
 124       enddo
 125 C#define DEBUG
 126 #ifdef MPI
 127       do i=1,scount(me1)
 128 #else
 129       do i=1,ntot(islice)
 130 #endif
 131 C        do j=1,nParmSet
 132 C          if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 133 C        enddo
 134         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 135         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 136         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 137         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 138         ind_point(i)=0
 139         do j=nQ,1,-1
 140           ind=(q(j,i)-dmin+1.0d-8)/delta
 141           if (j.eq.1) then
 142             ind_point(i)=ind_point(i)+ind
 143           else
 144             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 145           endif
 146 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 147 c     &      ind_point(i)
 148           call flush(iout)
 149           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 150             write (iout,*) "Error - index exceeds range for point",i,
 151      &      " q=",q(j,i)," ind",ind_point(i)
 152 #ifdef MPI
 153             write (iout,*) "Processor",me1
 154             call flush(iout)
 155             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 156 #endif
 157             stop
 158           endif
 159         enddo ! j
 160         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 161         htot(ind_point(i))=htot(ind_point(i))+1
 162 #ifdef DEBUG
 163         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
 164      &   " htot",htot(ind_point(i))
 165         call flush(iout)
 166 #endif
 167       enddo ! i
 168       call flush(iout)
 169
 170       nbin=nbin1**nQ-1
 171       write (iout,'(a)') "Numbers of counts in Q bins"
 172       do t=0,tmax
 173         if (htot(t).gt.0) then
 174         write (iout,'(i15,$)') t
 175         liczba=t
 176         do j=1,nQ
 177           jj = mod(liczba,nbin1)
 178           liczba=liczba/nbin1
 179           write (iout,'(i5,$)') jj
 180         enddo
 181         write (iout,'(i8)') htot(t)
 182         endif
 183       enddo
 184       do iparm=1,nParmSet
 185       write (iout,'(a,i3)') "Number of data points for parameter set",
 186      & iparm
 187       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
 188      & ib=1,nT_h(iparm))
 189       write (iout,'(i8)') stot(islice)
 190       write (iout,'(a)')
 191       enddo
 192       call flush(iout)
 193
 194 #ifdef MPI
 195       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
 196      &  WHAM_COMM,IERROR)
 197       tmax=tmax_t
 198 C      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
 199 C     &  MPI_MIN,WHAM_COMM,IERROR)
 200       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
 201      &  MPI_MIN,WHAM_COMM,IERROR)
 202       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
 203      &  MPI_MAX,WHAM_COMM,IERROR)
 204       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
 205      &  MPI_MIN,WHAM_COMM,IERROR)
 206       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
 207      &  MPI_MAX,WHAM_COMM,IERROR)
 208 C      potEmin=potEmin_t/2
 209
 210       rgymin=rgymin_t
 211       rgymax=rgymax_t
 212       rmsmin=rmsmin_t
 213       rmsmax=rmsmax_t
 214       write (iout,*) "potEmin",potEmin
 215 #endif
 216       rmsmin=deltrms*dint(rmsmin/deltrms)
 217       rmsmax=deltrms*dint(rmsmax/deltrms)
 218       rgymin=deltrms*dint(rgymin/deltrgy)
 219       rgymax=deltrms*dint(rgymax/deltrgy)
 220       nbin_rms=(rmsmax-rmsmin)/deltrms
 221       nbin_rgy=(rgymax-rgymin)/deltrgy
 222       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
 223      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 224       nFi=0
 225       do i=1,nParmSet
 226         do j=1,nT_h(i)
 227           nFi=nFi+nR(j,i)
 228         enddo
 229       enddo
 230       write (iout,*) "nFi",nFi
 231 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 232 ! simulations.
 233 #ifdef MPI
 234       do i=1,scount(me1)
 235 #else
 236       do i=1,ntot(islice)
 237 #endif
 238 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 239         do iparm=1,nParmSet
 240 #ifdef DEBUG
 241           write (iout,'(2i5,21f8.2)') i,iparm,
 242      &     (enetb(k,i,iparm),k=1,max_ene)
 243 #endif
 244           call restore_parm(iparm)
 245 #ifdef DEBUG
 246           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 247      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 248      &      wtor_d,wsccor,wbond
 249 #endif
 250           do ib=1,nT_h(iparm)
 251             if (rescale_mode.eq.1) then
 252               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 253               quotl=1.0d0
 254               kfacl=1.0d0
 255               do l=1,5
 256                 quotl1=quotl
 257                 quotl=quotl*quot
 258                 kfacl=kfacl*kfac
 259                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 260               enddo
 261 #if defined(FUNCTH)
 262               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 263               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 264 #elif defined(FUNCT)
 265               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 266 #else
 267               ft(6)=1.0d0
 268 #endif
 269             else if (rescale_mode.eq.2) then
 270               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 271               quotl=1.0d0
 272               do l=1,5
 273                 quotl=quotl*quot
 274                 fT(l)=1.12692801104297249644d0/
 275      &             dlog(dexp(quotl)+dexp(-quotl))
 276               enddo
 277 #if defined(FUNCTH)
 278               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 279               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 280 #elif defined(FUNCT)
 281               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 282 #else
 283               ft(6)=1.0d0
 284 #endif
 285 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 286             else if (rescale_mode.eq.0) then
 287               do l=1,6
 288                 fT(l)=1.0d0
 289               enddo
 290             else
 291               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 292      &          rescale_mode
 293               call flush(iout)
 294               return1
 295             endif
 296             evdw=enetb(1,i,iparm)
 297             evdw_t=enetb(21,i,iparm)
 298 #ifdef SCP14
 299             evdw2_14=enetb(17,i,iparm)
 300             evdw2=enetb(2,i,iparm)+evdw2_14
 301 #else
 302             evdw2=enetb(2,i,iparm)
 303             evdw2_14=0.0d0
 304 #endif
 305 #ifdef SPLITELE
 306             ees=enetb(3,i,iparm)
 307             evdw1=enetb(16,i,iparm)
 308 #else
 309             ees=enetb(3,i,iparm)
 310             evdw1=0.0d0
 311 #endif
 312             ecorr=enetb(4,i,iparm)
 313             ecorr5=enetb(5,i,iparm)
 314             ecorr6=enetb(6,i,iparm)
 315             eel_loc=enetb(7,i,iparm)
 316             eello_turn3=enetb(8,i,iparm)
 317             eello_turn4=enetb(9,i,iparm)
 318             eturn6=enetb(10,i,iparm)
 319             ebe=enetb(11,i,iparm)
 320             escloc=enetb(12,i,iparm)
 321             etors=enetb(13,i,iparm)
 322             etors_d=enetb(14,i,iparm)
 323             ehpb=enetb(15,i,iparm)
 324             estr=enetb(18,i,iparm)
 325             esccor=enetb(19,i,iparm)
 326             edihcnstr=enetb(20,i,iparm)
 327             eliptran=enetb(22,i,iparm)
 328             etube=enetb(25,i,iparm)
 329
 330
 331 #ifdef DEBUG
 332             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
 333      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 334      &       etors,etors_d,eello_turn3,eello_turn4,esccor
 335 #endif
 336
 337 #ifdef SPLITELE
 338             if (shield_mode.gt.0) then
 339             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 340      &      +ft(1)*welec*ees
 341      &      +ft(1)*wvdwpp*evdw1
 342      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 343      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 344      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 345      &      +ft(2)*wturn3*eello_turn3
 346      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 347      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 348      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 349              else
 350             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 351      &      +wvdwpp*evdw1
 352      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 353      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 354      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 355      &      +ft(2)*wturn3*eello_turn3
 356      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 357      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 358      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 359              endif
 360 #else
 361       if (shield_mode.gt.0) then
 362             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 363      &      +ft(1)*welec*(ees+evdw1)
 364      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 365      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 366      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 367      &      +ft(2)*wturn3*eello_turn3
 368      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 369      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 370      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 371             else
 372             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 373      &      +ft(1)*welec*(ees+evdw1)
 374      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 375      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 376      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 377      &      +ft(2)*wturn3*eello_turn3
 378      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 379      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 380      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 381            endif
 382
 383 #endif
 384 #ifdef DEBUG
 385             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 386      &        etot,potEmin
 387 #endif
 388 #ifdef DEBUG
 389             if (iparm.eq.1 .and. ib.eq.1) then
 390             write (iout,*)"Conformation",i
 391             energia(0)=etot
 392             do k=1,max_ene
 393               energia(k)=enetb(k,i,iparm)
 394             enddo
 395             call enerprint(energia(0),fT)
 396             endif
 397 #endif
 398             do kk=1,nR(ib,iparm)
 399               Econstr=0.0d0
 400               do j=1,nQ
 401                 dd = q(j,i)
 402                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 403      &           *(dd-q0(j,kk,ib,iparm))**2
 404               enddo
 405               v(i,kk,ib,iparm)=
 406      &          -beta_h(ib,iparm)*(etot+Econstr)
 407 #ifdef DEBUG
 408               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 409      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 410 #endif
 411             enddo ! kk
 412           enddo   ! ib
 413         enddo     ! iparm
 414       enddo       ! i
 415 ! Simple iteration to calculate free energies corresponding to all simulation
 416 ! runs.
 417       do iter=1,maxit
 418
 419 ! Compute new free-energy values corresponding to the righ-hand side of the
 420 ! equation and their derivatives.
 421         write (iout,*) "------------------------fi"
 422         entfac_min=1.0d10
 423
 424 #ifdef MPI
 425         do t=1,scount(me1)
 426 #else
 427         do t=1,ntot(islice)
 428 #endif
 429           vmax=-1.0d+20
 430           do i=1,nParmSet
 431             do k=1,nT_h(i)
 432               do l=1,nR(k,i)
 433                 vf=v(t,l,k,i)+f(l,k,i)
 434                 if (vf.gt.vmax) vmax=vf
 435               enddo
 436             enddo
 437           enddo
 438           denom=0.0d0
 439           do i=1,nParmSet
 440             do k=1,nT_h(i)
 441               do l=1,nR(k,i)
 442                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 443                 if (aux.gt.-200.0d0)
 444      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 445               enddo
 446             enddo
 447           enddo
 448           entfac(t)=-dlog(denom)-vmax
 449           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 450 #ifdef DEBUG
 451           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 452 #endif
 453         enddo
 454
 455         do iparm=1,nParmSet
 456           do iib=1,nT_h(iparm)
 457             do ii=1,nR(iib,iparm)
 458 #ifdef MPI
 459               fi_p_min(ii,iib,iparm)=-1.0d10
 460               do t=1,scount(me)
 461                 aux=v(t,ii,iib,iparm)+entfac(t)
 462                 if (aux.gt.fi_p_min(ii,iib,iparm))
 463      &                                   fi_p_min(ii,iib,iparm)=aux
 464               enddo
 465 #else
 466               do t=1,ntot(islice)
 467                 aux=v(t,ii,iib,iparm)+entfac(t)
 468                 if (aux.gt.fi_min(ii,iib,iparm))
 469      &                                     fi_min(ii,iib,iparm)=aux
 470               enddo
 471 #endif
 472             enddo ! ii
 473           enddo ! iib
 474         enddo ! iparm
 475 #ifdef MPI
 476 #ifdef DEBUG
 477         write (iout,*) "fi_min before AllReduce"
 478         do i=1,nParmSet
 479           do j=1,nT_h(i)
 480             write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
 481           enddo
 482         enddo
 483 #endif
 484         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
 485      &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 486 #ifdef DEBUG
 487         write (iout,*) "fi_min after AllReduce"
 488         do i=1,nParmSet
 489           do j=1,nT_h(i)
 490             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 491           enddo
 492         enddo
 493 #endif
 494 #endif
 495         do iparm=1,nParmSet
 496           do iib=1,nT_h(iparm)
 497             do ii=1,nR(iib,iparm)
 498 #ifdef MPI
 499               fi_p(ii,iib,iparm)=0.0d0
 500               do t=1,scount(me)
 501                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 502      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 503 #ifdef DEBUG
 504               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
 505      &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
 506      &         fi_p(ii,iib,iparm)
 507 #endif
 508               enddo
 509 #else
 510               fi(ii,iib,iparm)=0.0d0
 511               do t=1,ntot(islice)
 512                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 513      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 514               enddo
 515 #endif
 516             enddo ! ii
 517           enddo ! iib
 518         enddo ! iparm
 519 #ifdef MPI
 520 #ifdef DEBUG
 521         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 522         do iparm=1,nParmSet
 523           do ib=1,nT_h(nparmset)
 524             write (iout,*) "iparm",iparm," ib",ib
 525             write (iout,*) "beta=",beta_h(ib,iparm)
 526             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 527           enddo
 528         enddo
 529 #endif
 530         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 531      &   maxR*MaxT_h*nParmSet
 532         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 533      &   " WHAM_COMM",WHAM_COMM
 534         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 535      &   MPI_DOUBLE_PRECISION,
 536      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 537 #ifdef DEBUG
 538         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 539         do iparm=1,nParmSet
 540           write (iout,*) "iparm",iparm
 541           do ib=1,nT_h(iparm)
 542             write (iout,*) "beta=",beta_h(ib,iparm)
 543             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 544           enddo
 545         enddo
 546 #endif
 547         if (me1.eq.Master) then
 548 #endif
 549         avefi=0.0d0
 550         do iparm=1,nParmSet
 551           do ib=1,nT_h(iparm)
 552             do i=1,nR(ib,iparm)
 553               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 554               avefi=avefi+fi(i,ib,iparm)
 555             enddo
 556           enddo
 557         enddo
 558         avefi=avefi/nFi
 559         do iparm=1,nParmSet
 560           write (iout,*) "Parameter set",iparm
 561           do ib =1,nT_h(iparm)
 562             write (iout,*) "beta=",beta_h(ib,iparm)
 563             do i=1,nR(ib,iparm)
 564               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 565             enddo
 566             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 567             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 568           enddo
 569         enddo
 570
 571 ! Compute the norm of free-energy increments.
 572         finorm=0.0d0
 573         do iparm=1,nParmSet
 574           do ib=1,nT_h(iparm)
 575             do i=1,nR(ib,iparm)
 576               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 577               f(i,ib,iparm)=fi(i,ib,iparm)
 578             enddo
 579           enddo
 580         enddo
 581
 582         write (iout,*) 'Iteration',iter,' finorm',finorm
 583
 584 #ifdef MPI
 585         endif
 586         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 587      &   MPI_DOUBLE_PRECISION,Master,
 588      &   WHAM_COMM,IERROR)
 589         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 590      &   WHAM_COMM,IERROR)
 591 #endif
 592 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 593         if (finorm.lt.fimin) then
 594           write (iout,*) 'Iteration converged'
 595           goto 20
 596         endif
 597
 598       enddo ! iter
 599
 600    20 continue
 601
 602 #ifdef MPI
 603       do i=1,scount(me1)
 604 #else
 605       do i=1,ntot(islice)
 606 #endif
 607 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 608         do iparm=1,nParmSet
 609 #ifdef DEBUG
 610           write (iout,'(2i5,21f8.2)') i,iparm,
 611      &     (enetb(k,i,iparm),k=1,21)
 612 #endif
 613           call restore_parm(iparm)
 614 #ifdef DEBUG
 615           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 616      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 617      &      wtor_d,wsccor,wbond
 618 #endif
 619           do ib=1,nT_h(iparm)
 620             if (rescale_mode.eq.1) then
 621               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 622               quotl=1.0d0
 623               kfacl=1.0d0
 624               do l=1,5
 625                 quotl1=quotl
 626                 quotl=quotl*quot
 627                 kfacl=kfacl*kfac
 628                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 629               enddo
 630 #if defined(FUNCTH)
 631               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 632               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 633 #elif defined(FUNCT)
 634               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 635 #else
 636               ft(6)=1.0d0
 637 #endif
 638             else if (rescale_mode.eq.2) then
 639               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 640               quotl=1.0d0
 641               do l=1,5
 642                 quotl=quotl*quot
 643                 fT(l)=1.12692801104297249644d0/
 644      &             dlog(dexp(quotl)+dexp(-quotl))
 645               enddo
 646 #if defined(FUNCTH)
 647               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 648               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 649 #elif defined(FUNCT)
 650               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 651 #else
 652               ft(6)=1.0d0
 653 #endif
 654 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 655             else if (rescale_mode.eq.0) then
 656               do l=1,6
 657                 fT(l)=1.0d0
 658               enddo
 659             else
 660               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 661      &          rescale_mode
 662               call flush(iout)
 663               return1
 664             endif
 665           evdw=enetb(21,i,iparm)
 666           evdw_t=enetb(1,i,iparm)
 667 #ifdef SCP14
 668           evdw2_14=enetb(17,i,iparm)
 669           evdw2=enetb(2,i,iparm)+evdw2_14
 670 #else
 671           evdw2=enetb(2,i,iparm)
 672           evdw2_14=0.0d0
 673 #endif
 674 #ifdef SPLITELE
 675           ees=enetb(3,i,iparm)
 676           evdw1=enetb(16,i,iparm)
 677 #else
 678           ees=enetb(3,i,iparm)
 679           evdw1=0.0d0
 680 #endif
 681           ecorr=enetb(4,i,iparm)
 682           ecorr5=enetb(5,i,iparm)
 683           ecorr6=enetb(6,i,iparm)
 684           eel_loc=enetb(7,i,iparm)
 685           eello_turn3=enetb(8,i,iparm)
 686           eello_turn4=enetb(9,i,iparm)
 687           eturn6=enetb(10,i,iparm)
 688           ebe=enetb(11,i,iparm)
 689           escloc=enetb(12,i,iparm)
 690           etors=enetb(13,i,iparm)
 691           etors_d=enetb(14,i,iparm)
 692           ehpb=enetb(15,i,iparm)
 693           estr=enetb(18,i,iparm)
 694           esccor=enetb(19,i,iparm)
 695           edihcnstr=enetb(20,i,iparm)
 696 C          edihcnstr=0.0d0
 697           eliptran=enetb(22,i,iparm)
 698             etube=enetb(25,i,iparm)
 699
 700 #ifdef SPLITELE
 701       if (shield_mode.gt.0) then
 702             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 703      &      +ft(1)*welec*ees
 704      &      +ft(1)*wvdwpp*evdw1
 705      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 706      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 707      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 708      &      +ft(2)*wturn3*eello_turn3
 709      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 710      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 711      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 712         else
 713             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 714      &      +wvdwpp*evdw1
 715      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 716      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 717      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 718      &      +ft(2)*wturn3*eello_turn3
 719      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 720      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 721      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 722       endif
 723 #else
 724       if (shield_mode.gt.0) then
 725             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 726      &      +ft(1)*welec*(ees+evdw1)
 727      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 728      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 729      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 730      &      +ft(2)*wturn3*eello_turn3
 731      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 732      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 733      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 734        else
 735             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 736      &      +ft(1)*welec*(ees+evdw1)
 737      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 738      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 739      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 740      &      +ft(2)*wturn3*eello_turn3
 741      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 742      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 743      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 744        endif
 745
 746 #endif
 747             etot=etot-entfac(i)/beta_h(ib,iparm)
 748             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 749
 750           enddo   ! ib
 751         enddo     ! iparm
 752       enddo       ! i
 753
 754
 755 #ifdef DEBUG
 756       write (iout,*) "The potEmin array before reduction"
 757       do i=1,nParmSet
 758         write (iout,*) "Parameter set",i
 759         do j=1,nT_h(i)
 760           write (iout,*) j,PotEmin_all(j,i)
 761         enddo
 762       enddo
 763       write (iout,*) "potEmin_min",potEmin_min
 764 #endif
 765 #ifdef MPI
 766 C Determine the minimum energes for all parameter sets and temperatures
 767       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
 768      &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 769       do i=1,nParmSet
 770         do j=1,nT_h(i)
 771           potEmin_all(j,i)=potEmin_t_all(j,i)
 772         enddo
 773       enddo
 774 #endif
 775       potEmin_min=potEmin_all(1,1)
 776       do i=1,nParmSet
 777         do j=1,nT_h(i)
 778           if (potEmin_all(j,i).lt.potEmin_min)
 779      &             potEmin_min=potEmin_all(j,i)
 780         enddo
 781       enddo
 782 #ifdef DEBUG
 783       write (iout,*) "The potEmin array"
 784       do i=1,nParmSet
 785         write (iout,*) "Parameter set",i
 786         do j=1,nT_h(i)
 787           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
 788      &      PotEmin_all(j,i)
 789         enddo
 790       enddo
 791       write (iout,*) "potEmin_min",potEmin_min
 792 #endif
 793
 794
 795 ! Now, put together the histograms from all simulations, in order to get the
 796 ! unbiased total histogram.
 797 #ifdef MPI
 798       do t=0,tmax
 799         hfin_ent_p(t)=0.0d0
 800       enddo
 801 #else
 802       do t=0,tmax
 803         hfin_ent(t)=0.0d0
 804       enddo
 805 #endif
 806       write (iout,*) "--------------hist"
 807 #ifdef MPI
 808       do iparm=1,nParmSet
 809         do i=0,nGridT
 810           sumW_p(i,iparm)=0.0d0
 811           sumE_p(i,iparm)=0.0d0
 812           sumEbis_p(i,iparm)=0.0d0
 813           sumEsq_p(i,iparm)=0.0d0
 814           do j=1,nQ+2
 815             sumQ_p(j,i,iparm)=0.0d0
 816             sumQsq_p(j,i,iparm)=0.0d0
 817             sumEQ_p(j,i,iparm)=0.0d0
 818           enddo
 819         enddo
 820       enddo
 821       upindE_p=0
 822 #else
 823       do iparm=1,nParmSet
 824         do i=0,nGridT
 825           sumW(i,iparm)=0.0d0
 826           sumE(i,iparm)=0.0d0
 827           sumEbis(i,iparm)=0.0d0
 828           sumEsq(i,iparm)=0.0d0
 829           do j=1,nQ+2
 830             sumQ(j,i,iparm)=0.0d0
 831             sumQsq(j,i,iparm)=0.0d0
 832             sumEQ(j,i,iparm)=0.0d0
 833           enddo
 834         enddo
 835       enddo
 836       upindE=0
 837 #endif
 838 c 8/26/05 entropy distribution
 839 #ifdef MPI
 840       entmin_p=1.0d10
 841       entmax_p=-1.0d10
 842       do t=1,scount(me1)
 843 c        ent=-dlog(entfac(t))
 844         ent=entfac(t)
 845         if (ent.lt.entmin_p) entmin_p=ent
 846         if (ent.gt.entmax_p) entmax_p=ent
 847       enddo
 848       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 849       call flush(iout)
 850       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 851      &  WHAM_COMM,IERROR)
 852       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 853      &  WHAM_COMM,IERROR)
 854 C      ientmax=entmax-entmin
 855 C      if (ientmax.gt.2000) ientmax=2000
 856       if ((-dlog(entmax)-entmin).lt.2000.0d0) then
 857       ientmax=-dlog(entmax)-entmin
 858       else
 859        ientmax=2000
 860       endif
 861       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 862       call flush(iout)
 863       do t=1,scount(me1)
 864 c        ient=-dlog(entfac(t))-entmin
 865         ient=entfac(t)-entmin
 866       write (iout,*) "ient",ient,entfac(t),entmin
 867 C        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 868       enddo
 869 C      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 870 C     &  MPI_SUM,WHAM_COMM,IERROR)
 871 C      if (me1.eq.Master) then
 872 C        write (iout,*) "Entropy histogram"
 873 C        do i=0,ientmax
 874 C          write(iout,'(f15.4,i10)') entmin+i,histent(i)
 875 C        enddo
 876 C      endif
 877 #else
 878       entmin=1.0d10
 879       entmax=-1.0d10
 880       do t=1,ntot(islice)
 881         ent=entfac(t)
 882         if (ent.lt.entmin) entmin=ent
 883         if (ent.gt.entmax) entmax=ent
 884       enddo
 885       if ((-dlog(entmax)-entmin).lt.2000.0d0) then
 886       ientmax=-dlog(entmax)-entmin
 887       else
 888        ientmax=2000
 889       endif
 890 C      do t=1,ntot(islice)
 891 C        ient=entfac(t)-entmin
 892 C        if (ient.le.2000) histent(ient)=histent(ient)+1
 893 C      enddo
 894 C      write (iout,*) "Entropy histogram"
 895 C      do i=0,ientmax
 896 C        write(iout,'(2f15.4)') entmin+i,histent(i)
 897 C      enddo
 898 #endif
 899
 900 #ifdef MPI
 901 c      write (iout,*) "me1",me1," scount",scount(me1)
 902
 903       do iparm=1,nParmSet
 904
 905 #ifdef MPI
 906         do ib=1,nT_h(iparm)
 907           do t=0,tmax
 908             hfin_p(t,ib)=0.0d0
 909           enddo
 910         enddo
 911         do i=1,maxindE
 912           histE_p(i)=0.0d0
 913         enddo
 914 #else
 915         do ib=1,nT_h(iparm)
 916           do t=0,tmax
 917             hfin(t,ib)=0.0d0
 918           enddo
 919         enddo
 920         do i=1,maxindE
 921           histE(i)=0.0d0
 922         enddo
 923 #endif
 924         do ib=1,nT_h(iparm)
 925           do i=0,MaxBinRms
 926             do j=0,MaxBinRgy
 927               hrmsrgy(j,i,ib)=0.0d0
 928 #ifdef MPI
 929               hrmsrgy_p(j,i,ib)=0.0d0
 930 #endif
 931             enddo
 932           enddo
 933         enddo
 934
 935         do t=1,scount(me1)
 936 #else
 937         do t=1,ntot(islice)
 938 #endif
 939           ind = ind_point(t)
 940 #ifdef MPI
 941           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 942 #else
 943           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 944 #endif
 945 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 946           call restore_parm(iparm)
 947           evdw=enetb(21,t,iparm)
 948           evdw_t=enetb(1,t,iparm)
 949 #ifdef SCP14
 950           evdw2_14=enetb(17,t,iparm)
 951           evdw2=enetb(2,t,iparm)+evdw2_14
 952 #else
 953           evdw2=enetb(2,t,iparm)
 954           evdw2_14=0.0d0
 955 #endif
 956 #ifdef SPLITELE
 957           ees=enetb(3,t,iparm)
 958           evdw1=enetb(16,t,iparm)
 959 #else
 960           ees=enetb(3,t,iparm)
 961           evdw1=0.0d0
 962 #endif
 963           ecorr=enetb(4,t,iparm)
 964           ecorr5=enetb(5,t,iparm)
 965           ecorr6=enetb(6,t,iparm)
 966           eel_loc=enetb(7,t,iparm)
 967           eello_turn3=enetb(8,t,iparm)
 968           eello_turn4=enetb(9,t,iparm)
 969           eturn6=enetb(10,t,iparm)
 970           ebe=enetb(11,t,iparm)
 971           escloc=enetb(12,t,iparm)
 972           etors=enetb(13,t,iparm)
 973           etors_d=enetb(14,t,iparm)
 974           ehpb=enetb(15,t,iparm)
 975           estr=enetb(18,t,iparm)
 976           esccor=enetb(19,t,iparm)
 977           edihcnstr=enetb(20,t,iparm)
 978 C          edihcnstr=0.0d0
 979           eliptran=enetb(22,t,iparm)
 980             etube=enetb(25,t,iparm)
 981
 982           do k=0,nGridT
 983             betaT=startGridT+k*delta_T
 984             temper=betaT
 985 c            fT=T0/betaT
 986 c            ft=2*T0/(T0+betaT)
 987             if (rescale_mode.eq.1) then
 988               quot=betaT/T0
 989               quotl=1.0d0
 990               kfacl=1.0d0
 991               do l=1,5
 992                 quotl1=quotl
 993                 quotl=quotl*quot
 994                 kfacl=kfacl*kfac
 995                 denom=kfacl-1.0d0+quotl
 996                 fT(l)=kfacl/denom
 997                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 998                 ftbis(l)=l*kfacl*quotl1*
 999      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1000               enddo
1001 #if defined(FUNCTH)
1002               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1003      &                  320.0d0
1004               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1005              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1006      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1007 #elif defined(FUNCT)
1008               fT(6)=betaT/T0
1009               ftprim(6)=1.0d0/T0
1010               ftbis(6)=0.0d0
1011 #else
1012               fT(6)=1.0d0
1013               ftprim(6)=0.0d0
1014               ftbis(6)=0.0d0
1015 #endif
1016             else if (rescale_mode.eq.2) then
1017               quot=betaT/T0
1018               quotl=1.0d0
1019               do l=1,5
1020                 quotl1=quotl
1021                 quotl=quotl*quot
1022                 eplus=dexp(quotl)
1023                 eminus=dexp(-quotl)
1024                 logfac=1.0d0/dlog(eplus+eminus)
1025                 tanhT=(eplus-eminus)/(eplus+eminus)
1026                 fT(l)=1.12692801104297249644d0*logfac
1027                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1028                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1029      &          2*l*quotl1/T0*logfac*
1030      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1031      &          +ftprim(l)*tanhT)
1032               enddo
1033 #if defined(FUNCTH)
1034               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1035      &                 320.0d0
1036               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1037              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1038      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1039 #elif defined(FUNCT)
1040               fT(6)=betaT/T0
1041               ftprim(6)=1.0d0/T0
1042               ftbis(6)=0.0d0
1043 #else
1044               fT(6)=1.0d0
1045               ftprim(6)=0.0d0
1046               ftbis(6)=0.0d0
1047 #endif
1048             else if (rescale_mode.eq.0) then
1049               do l=1,5
1050                 fT(l)=1.0d0
1051                 ftprim(l)=0.0d0
1052               enddo
1053             else
1054               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1055      &          rescale_mode
1056               call flush(iout)
1057               return1
1058             endif
1059 c            write (iout,*) "ftprim",ftprim
1060 c            write (iout,*) "ftbis",ftbis
1061             betaT=1.0d0/(1.987D-3*betaT)
1062             if (betaT.ge.beta_h(1,iparm)) then
1063               potEmin=potEmin_all(1,iparm)+
1064      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1065      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1066      &         (1.0/betaT-1.0/beta_h(1,iparm))
1067 #ifdef DEBUG
1068               write(iout,*) "first",temper,potEmin
1069 #endif
1070             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1071               potEmin=potEmin_all(nT_h(iparm),iparm)+
1072      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1073      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1074      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1075 #ifdef DEBUG
1076               write (iout,*) "last",temper,potEmin
1077 #endif
1078             else
1079               do l=1,nT_h(iparm)-1
1080                 if (betaT.le.beta_h(l,iparm) .and.
1081      &              betaT.gt.beta_h(l+1,iparm)) then
1082                   potEmin=potEmin_all(l,iparm)
1083 #ifdef DEBUG
1084                   write (iout,*) "l",l,
1085      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1086      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1087 #endif
1088                   exit
1089                 endif
1090               enddo
1091             endif
1092 #ifdef SPLITELE
1093       if (shield_mode.gt.0) then
1094             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1095      &      +ft(1)*welec*ees
1096      &      +ft(1)*wvdwpp*evdw1
1097      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1098      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1099      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1100      &      +ft(2)*wturn3*eello_turn3
1101      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1102      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1103      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1104             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1105 C     &            +ftprim(6)*evdw_t
1106      &            +ftprim(1)*wscp*evdw2
1107      &            +ftprim(1)*welec*ees
1108      &            +ftprim(1)*wvdwpp*evdw1
1109      &            +ftprim(1)*wtor*etors+
1110      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1111      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1112      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1113      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1114      &            ftprim(1)*wsccor*esccor
1115             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1116      &            +ftbis(1)*wscp*evdw2+
1117      &            ftbis(1)*welec*ees
1118      &            +ftbis(1)*wvdwpp*evdw
1119      &            +ftbis(1)*wtor*etors+
1120      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1121      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1122      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1123      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1124      &            ftbis(1)*wsccor*esccor
1125         else
1126             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1127      &      +wvdwpp*evdw1
1128      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1129      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1130      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1131      &      +ft(2)*wturn3*eello_turn3
1132      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1133      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1134      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1135             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1136      &            +ftprim(1)*wtor*etors+
1137      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1138      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1139      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1140      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1141      &            ftprim(1)*wsccor*esccor
1142             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1143      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1144      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1145      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1146      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1147      &            ftbis(1)*wsccor*esccor
1148       endif
1149 #else
1150       if (shield_mode.gt.0) then
1151             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1152      &      +ft(1)*welec*(ees+evdw1)
1153      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1154      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1155      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1156      &      +ft(2)*wturn3*eello_turn3
1157      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1158      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1159      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1160             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1161      &           +ftprim(1)*welec*(ees+evdw1)
1162      &           +ftprim(1)*wtor*etors+
1163      &            ftprim(1)*wscp*evdw2+
1164      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1165      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1166      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1167      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1168      &            ftprim(1)*wsccor*esccor
1169             ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1170      &           +ftbis(1)*wscp*evdw2
1171      &           +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1172      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1173      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1174      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1175      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1176      &            ftbis(1)*wsccor*esccor
1177        else
1178             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1179      &      +ft(1)*welec*(ees+evdw1)
1180      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1181      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1182      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1183      &      +ft(2)*wturn3*eello_turn3
1184      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1185      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1186      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1187             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1188      &           +ftprim(1)*wtor*etors+
1189      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1190      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1191      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1192      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1193      &            ftprim(1)*wsccor*esccor
1194             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1195      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1196      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1197      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1198      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1199      &            ftbis(1)*wsccor*esccor
1200
1201        endif
1202
1203 #endif
1204
1205             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1206 #ifdef DEBUG
1207             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1208      &        " etot",etot," entfac",entfac(t),
1209      &        " weight",weight," ebis",ebis
1210 #endif
1211             etot=etot-temper*eprim
1212 #ifdef MPI
1213             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1214             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1215             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1216             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1217             do j=1,nQ+2
1218               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1219               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1220               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1221      &         +etot*q(j,t)*weight
1222             enddo
1223 #else
1224             sumW(k,iparm)=sumW(k,iparm)+weight
1225             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1226             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1227             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1228             do j=1,nQ+2
1229               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1230               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1231               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1232      &         +etot*q(j,t)*weight
1233             enddo
1234 #endif
1235           enddo
1236           indE = aint(potE(t,iparm)-aint(potEmin))
1237           if (indE.ge.0 .and. indE.le.maxinde) then
1238             if (indE.gt.upindE_p) upindE_p=indE
1239             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1240           endif
1241 #ifdef MPI
1242           do ib=1,nT_h(iparm)
1243             potEmin=potEmin_all(ib,iparm)
1244             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1245             hfin_p(ind,ib)=hfin_p(ind,ib)+
1246      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1247              if (rmsrgymap) then
1248                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1249                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1250                hrmsrgy_p(indrgy,indrms,ib)=
1251      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1252              endif
1253           enddo
1254 #else
1255           do ib=1,nT_h(iparm)
1256             potEmin=potEmin_all(ib,iparm)
1257             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1258             hfin(ind,ib)=hfin(ind,ib)+
1259      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1260              if (rmsrgymap) then
1261                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1262                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1263                hrmsrgy(indrgy,indrms,ib)=
1264      &           hrmsrgy(indrgy,indrms,ib)+expfac
1265              endif
1266           enddo
1267 #endif
1268         enddo ! t
1269         do ib=1,nT_h(iparm)
1270           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1271      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1272           if (rmsrgymap) then
1273           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1274      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1275      &       WHAM_COMM,IERROR)
1276           endif
1277         enddo
1278         call MPI_Reduce(upindE_p,upindE,1,
1279      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1280         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1281      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1282
1283         if (me1.eq.master) then
1284
1285         if (histout) then
1286
1287         write (iout,'(6x,$)')
1288         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1289      &   ib=1,nT_h(iparm))
1290         write (iout,*)
1291
1292         write (iout,'(/a)') 'Final histograms'
1293         if (histfile) then
1294           if (nslice.eq.1) then
1295             if (separate_parset) then
1296               write(licz3,"(bz,i3.3)") myparm
1297               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1298             else
1299               histname=prefix(:ilen(prefix))//'.hist'
1300             endif
1301           else
1302             if (separate_parset) then
1303               write(licz3,"(bz,i3.3)") myparm
1304               histname=prefix(:ilen(prefix))//'_par'//licz3//
1305      &         '_slice_'//licz2//'.hist'
1306             else
1307               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1308             endif
1309           endif
1310 #if defined(AIX) || defined(PGI)
1311           open (ihist,file=histname,position='append')
1312 #else
1313           open (ihist,file=histname,access='append')
1314 #endif
1315         endif
1316
1317         do t=0,tmax
1318           liczba=t
1319           sumH=0.0d0
1320           do ib=1,nT_h(iparm)
1321             sumH=sumH+hfin(t,ib)
1322           enddo
1323           if (sumH.gt.0.0d0) then
1324             do j=1,nQ
1325               jj = mod(liczba,nbin1)
1326               liczba=liczba/nbin1
1327               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1328               if (histfile)
1329      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1330             enddo
1331             do ib=1,nT_h(iparm)
1332               write (iout,'(e20.10,$)') hfin(t,ib)
1333               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1334             enddo
1335             write (iout,'(i5)') iparm
1336             if (histfile) write (ihist,'(i5)') iparm
1337           endif
1338         enddo
1339
1340         endif
1341
1342         if (entfile) then
1343           if (nslice.eq.1) then
1344             if (separate_parset) then
1345               write(licz3,"(bz,i3.3)") myparm
1346               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1347             else
1348               histname=prefix(:ilen(prefix))//'.ent'
1349             endif
1350           else
1351             if (separate_parset) then
1352               write(licz3,"(bz,i3.3)") myparm
1353               histname=prefix(:ilen(prefix))//'par_'//licz3//
1354      &           '_slice_'//licz2//'.ent'
1355             else
1356               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1357             endif
1358           endif
1359 #if defined(AIX) || defined(PGI)
1360           open (ihist,file=histname,position='append')
1361 #else
1362           open (ihist,file=histname,access='append')
1363 #endif
1364           write (ihist,'(a)') "# Microcanonical entropy"
1365           do i=0,upindE
1366             write (ihist,'(f8.0,$)') dint(potEmin)+i
1367             if (histE(i).gt.0.0e0) then
1368               write (ihist,'(f15.5,$)') dlog(histE(i))
1369             else
1370               write (ihist,'(f15.5,$)') 0.0d0
1371             endif
1372           enddo
1373           write (ihist,*)
1374           close(ihist)
1375         endif
1376         write (iout,*) "Microcanonical entropy"
1377         do i=0,upindE
1378           write (iout,'(f8.0,$)') dint(potEmin)+i
1379           if (histE(i).gt.0.0e0) then
1380             write (iout,'(f15.5,$)') dlog(histE(i))
1381           else
1382             write (iout,'(f15.5,$)') 0.0d0
1383           endif
1384           write (iout,*)
1385         enddo
1386         if (rmsrgymap) then
1387           if (nslice.eq.1) then
1388             if (separate_parset) then
1389               write(licz3,"(bz,i3.3)") myparm
1390               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1391             else
1392               histname=prefix(:ilen(prefix))//'.rmsrgy'
1393             endif
1394           else
1395             if (separate_parset) then
1396               write(licz3,"(bz,i3.3)") myparm
1397               histname=prefix(:ilen(prefix))//'_par'//licz3//
1398      &         '_slice_'//licz2//'.rmsrgy'
1399             else
1400              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1401             endif
1402           endif
1403 #if defined(AIX) || defined(PGI)
1404           open (ihist,file=histname,position='append')
1405 #else
1406           open (ihist,file=histname,access='append')
1407 #endif
1408           do i=0,nbin_rms
1409             do j=0,nbin_rgy
1410               write(ihist,'(2f8.2,$)')
1411      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1412               do ib=1,nT_h(iparm)
1413                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1414                   write(ihist,'(e14.5,$)')
1415      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1416      &              +potEmin
1417                 else
1418                   write(ihist,'(e14.5,$)') 1.0d6
1419                 endif
1420               enddo
1421               write (ihist,'(i2)') iparm
1422             enddo
1423           enddo
1424           close(ihist)
1425         endif
1426         endif
1427       enddo ! iparm
1428 #ifdef MPI
1429       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1430      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1431       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1432      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1433       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1434      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1435       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1436      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1437       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1438      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1439       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1440      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1441      &   WHAM_COMM,IERROR)
1442       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1443      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1444      &   WHAM_COMM,IERROR)
1445       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1446      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1447      &   WHAM_COMM,IERROR)
1448       if (me.eq.master) then
1449 #endif
1450       write (iout,'(/a)') 'Thermal characteristics of folding'
1451       if (nslice.eq.1) then
1452         nazwa=prefix
1453       else
1454         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1455       endif
1456       iln=ilen(nazwa)
1457       if (nparmset.eq.1 .and. .not.separate_parset) then
1458         nazwa=nazwa(:iln)//".thermal"
1459       else if (nparmset.eq.1 .and. separate_parset) then
1460         write(licz3,"(bz,i3.3)") myparm
1461         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1462       endif
1463       do iparm=1,nParmSet
1464       if (nparmset.gt.1) then
1465         write(licz3,"(bz,i3.3)") iparm
1466         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1467       endif
1468       open(34,file=nazwa)
1469       if (separate_parset) then
1470         write (iout,'(a,i3)') "Parameter set",myparm
1471       else
1472         write (iout,'(a,i3)') "Parameter set",iparm
1473       endif
1474       do i=0,NGridT
1475         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1476         if (betaT.ge.beta_h(1,iparm)) then
1477           potEmin=potEmin_all(1,iparm)
1478         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1479           potEmin=potEmin_all(nT_h(iparm),iparm)
1480         else
1481           do l=1,nT_h(iparm)-1
1482             if (betaT.le.beta_h(l,iparm) .and.
1483      &          betaT.gt.beta_h(l+1,iparm)) then
1484               potEmin=potEmin_all(l,iparm)
1485               exit
1486             endif
1487           enddo
1488         endif
1489
1490         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1491         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1492      &    sumW(i,iparm)
1493         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1494      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1495         do j=1,nQ+2
1496           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1497           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1498      &     -sumQ(j,i,iparm)**2
1499           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1500      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1501         enddo
1502         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1503      &   (startGridT+i*delta_T))+potEmin
1504         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1505      &   sumW(i,iparm),sumE(i,iparm)
1506         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1507         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1508      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1509         write (iout,*)
1510         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1511      &   sumW(i,iparm),sumE(i,iparm)
1512         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1513         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1514      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1515         write (34,*)
1516         call flush(34)
1517       enddo
1518       close(34)
1519       enddo
1520       if (histout) then
1521       do t=0,tmax
1522         if (hfin_ent(t).gt.0.0d0) then
1523           liczba=t
1524           jj = mod(liczba,nbin1)
1525           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1526      &     hfin_ent(t)
1527           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1528      &      dmin+(jj+0.5d0)*delta,
1529      &     hfin_ent(t)
1530         endif
1531       enddo
1532       if (histfile) close(ihist)
1533       endif
1534
1535 #ifdef ZSCORE
1536 ! Write data for zscore
1537       if (nslice.eq.1) then
1538         zscname=prefix(:ilen(prefix))//".zsc"
1539       else
1540         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1541       endif
1542 #if defined(AIX) || defined(PGI)
1543       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1544 #else
1545       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1546 #endif
1547       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1548       do iparm=1,nParmSet
1549         write (izsc,'("NT=",i1)') nT_h(iparm)
1550       do ib=1,nT_h(iparm)
1551         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1552      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1553         jj = min0(nR(ib,iparm),7)
1554         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1555         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1556         write (izsc,'("&")')
1557         if (nR(ib,iparm).gt.7) then
1558           do ii=8,nR(ib,iparm),9
1559             jj = min0(nR(ib,iparm),ii+8)
1560             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1561             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1562             write (izsc,'("&")')
1563           enddo
1564         endif
1565         write (izsc,'("FI=",$)')
1566         jj=min0(nR(ib,iparm),7)
1567         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1568         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1569         write (izsc,'("&")')
1570         if (nR(ib,iparm).gt.7) then
1571           do ii=8,nR(ib,iparm),9
1572             jj = min0(nR(ib,iparm),ii+8)
1573             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1574             if (jj.eq.nR(ib,iparm)) then
1575               write (izsc,*)
1576             else
1577               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1578               write (izsc,'(t80,"&")')
1579             endif
1580           enddo
1581         endif
1582         do i=1,nR(ib,iparm)
1583           write (izsc,'("KH=",$)')
1584           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1585           write (izsc,'(" Q0=",$)')
1586           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1587           write (izsc,*)
1588         enddo
1589       enddo
1590       enddo
1591       close(izsc)
1592 #endif
1593 #ifdef MPI
1594       endif
1595 #endif
1596
1597       return
1598 C#undef DEBUG
1599       end