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,22)
 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       ientmax=entmax-entmin
 855       if (ientmax.gt.2000) ientmax=2000
 856       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 857       call flush(iout)
 858       do t=1,scount(me1)
 859 c        ient=-dlog(entfac(t))-entmin
 860         ient=entfac(t)-entmin
 861         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 862       enddo
 863       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 864      &  MPI_SUM,WHAM_COMM,IERROR)
 865       if (me1.eq.Master) then
 866         write (iout,*) "Entropy histogram"
 867         do i=0,ientmax
 868           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 869         enddo
 870       endif
 871 #else
 872       entmin=1.0d10
 873       entmax=-1.0d10
 874       do t=1,ntot(islice)
 875         ent=entfac(t)
 876         if (ent.lt.entmin) entmin=ent
 877         if (ent.gt.entmax) entmax=ent
 878       enddo
 879       ientmax=-dlog(entmax)-entmin
 880       if (ientmax.gt.2000) ientmax=2000
 881       do t=1,ntot(islice)
 882         ient=entfac(t)-entmin
 883         if (ient.le.2000) histent(ient)=histent(ient)+1
 884       enddo
 885       write (iout,*) "Entropy histogram"
 886       do i=0,ientmax
 887         write(iout,'(2f15.4)') entmin+i,histent(i)
 888       enddo
 889 #endif
 890
 891 #ifdef MPI
 892 c      write (iout,*) "me1",me1," scount",scount(me1)
 893
 894       do iparm=1,nParmSet
 895
 896 #ifdef MPI
 897         do ib=1,nT_h(iparm)
 898           do t=0,tmax
 899             hfin_p(t,ib)=0.0d0
 900           enddo
 901         enddo
 902         do i=1,maxindE
 903           histE_p(i)=0.0d0
 904         enddo
 905 #else
 906         do ib=1,nT_h(iparm)
 907           do t=0,tmax
 908             hfin(t,ib)=0.0d0
 909           enddo
 910         enddo
 911         do i=1,maxindE
 912           histE(i)=0.0d0
 913         enddo
 914 #endif
 915         do ib=1,nT_h(iparm)
 916           do i=0,MaxBinRms
 917             do j=0,MaxBinRgy
 918               hrmsrgy(j,i,ib)=0.0d0
 919 #ifdef MPI
 920               hrmsrgy_p(j,i,ib)=0.0d0
 921 #endif
 922             enddo
 923           enddo
 924         enddo
 925
 926         do t=1,scount(me1)
 927 #else
 928         do t=1,ntot(islice)
 929 #endif
 930           ind = ind_point(t)
 931 #ifdef MPI
 932           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 933 #else
 934           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 935 #endif
 936 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 937           call restore_parm(iparm)
 938           evdw=enetb(21,t,iparm)
 939           evdw_t=enetb(1,t,iparm)
 940 #ifdef SCP14
 941           evdw2_14=enetb(17,t,iparm)
 942           evdw2=enetb(2,t,iparm)+evdw2_14
 943 #else
 944           evdw2=enetb(2,t,iparm)
 945           evdw2_14=0.0d0
 946 #endif
 947 #ifdef SPLITELE
 948           ees=enetb(3,t,iparm)
 949           evdw1=enetb(16,t,iparm)
 950 #else
 951           ees=enetb(3,t,iparm)
 952           evdw1=0.0d0
 953 #endif
 954           ecorr=enetb(4,t,iparm)
 955           ecorr5=enetb(5,t,iparm)
 956           ecorr6=enetb(6,t,iparm)
 957           eel_loc=enetb(7,t,iparm)
 958           eello_turn3=enetb(8,t,iparm)
 959           eello_turn4=enetb(9,t,iparm)
 960           eturn6=enetb(10,t,iparm)
 961           ebe=enetb(11,t,iparm)
 962           escloc=enetb(12,t,iparm)
 963           etors=enetb(13,t,iparm)
 964           etors_d=enetb(14,t,iparm)
 965           ehpb=enetb(15,t,iparm)
 966           estr=enetb(18,t,iparm)
 967           esccor=enetb(19,t,iparm)
 968           edihcnstr=enetb(20,t,iparm)
 969 C          edihcnstr=0.0d0
 970           eliptran=enetb(22,i,iparm)
 971             etube=enetb(25,i,iparm)
 972
 973           do k=0,nGridT
 974             betaT=startGridT+k*delta_T
 975             temper=betaT
 976 c            fT=T0/betaT
 977 c            ft=2*T0/(T0+betaT)
 978             if (rescale_mode.eq.1) then
 979               quot=betaT/T0
 980               quotl=1.0d0
 981               kfacl=1.0d0
 982               do l=1,5
 983                 quotl1=quotl
 984                 quotl=quotl*quot
 985                 kfacl=kfacl*kfac
 986                 denom=kfacl-1.0d0+quotl
 987                 fT(l)=kfacl/denom
 988                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 989                 ftbis(l)=l*kfacl*quotl1*
 990      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 991               enddo
 992 #if defined(FUNCTH)
 993               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 994      &                  320.0d0
 995               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 996              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 997      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 998 #elif defined(FUNCT)
 999               fT(6)=betaT/T0
1000               ftprim(6)=1.0d0/T0
1001               ftbis(6)=0.0d0
1002 #else
1003               fT(6)=1.0d0
1004               ftprim(6)=0.0d0
1005               ftbis(6)=0.0d0
1006 #endif
1007             else if (rescale_mode.eq.2) then
1008               quot=betaT/T0
1009               quotl=1.0d0
1010               do l=1,5
1011                 quotl1=quotl
1012                 quotl=quotl*quot
1013                 eplus=dexp(quotl)
1014                 eminus=dexp(-quotl)
1015                 logfac=1.0d0/dlog(eplus+eminus)
1016                 tanhT=(eplus-eminus)/(eplus+eminus)
1017                 fT(l)=1.12692801104297249644d0*logfac
1018                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1019                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1020      &          2*l*quotl1/T0*logfac*
1021      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1022      &          +ftprim(l)*tanhT)
1023               enddo
1024 #if defined(FUNCTH)
1025               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1026      &                 320.0d0
1027               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1028              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1029      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1030 #elif defined(FUNCT)
1031               fT(6)=betaT/T0
1032               ftprim(6)=1.0d0/T0
1033               ftbis(6)=0.0d0
1034 #else
1035               fT(6)=1.0d0
1036               ftprim(6)=0.0d0
1037               ftbis(6)=0.0d0
1038 #endif
1039             else if (rescale_mode.eq.0) then
1040               do l=1,5
1041                 fT(l)=1.0d0
1042                 ftprim(l)=0.0d0
1043               enddo
1044             else
1045               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1046      &          rescale_mode
1047               call flush(iout)
1048               return1
1049             endif
1050 c            write (iout,*) "ftprim",ftprim
1051 c            write (iout,*) "ftbis",ftbis
1052             betaT=1.0d0/(1.987D-3*betaT)
1053             if (betaT.ge.beta_h(1,iparm)) then
1054               potEmin=potEmin_all(1,iparm)+
1055      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1056      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1057      &         (1.0/betaT-1.0/beta_h(1,iparm))
1058 #ifdef DEBUG
1059               write(iout,*) "first",temper,potEmin
1060 #endif
1061             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1062               potEmin=potEmin_all(nT_h(iparm),iparm)+
1063      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1064      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1065      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1066 #ifdef DEBUG
1067               write (iout,*) "last",temper,potEmin
1068 #endif
1069             else
1070               do l=1,nT_h(iparm)-1
1071                 if (betaT.le.beta_h(l,iparm) .and.
1072      &              betaT.gt.beta_h(l+1,iparm)) then
1073                   potEmin=potEmin_all(l,iparm)
1074 #ifdef DEBUG
1075                   write (iout,*) "l",l,
1076      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1077      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1078 #endif
1079                   exit
1080                 endif
1081               enddo
1082             endif
1083 #ifdef SPLITELE
1084       if (shield_mode.gt.0) then
1085             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1086      &      +ft(1)*welec*ees
1087      &      +ft(1)*wvdwpp*evdw1
1088      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1089      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1090      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1091      &      +ft(2)*wturn3*eello_turn3
1092      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1093      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1094      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1095             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1096 C     &            +ftprim(6)*evdw_t
1097      &            +ftprim(1)*wscp*evdw2
1098      &            +ftprim(1)*welec*ees
1099      &            +ftprim(1)*wvdwpp*evdw1
1100      &            +ftprim(1)*wtor*etors+
1101      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1102      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1103      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1104      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1105      &            ftprim(1)*wsccor*esccor
1106             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1107      &            +ftbis(1)*wscp*evdw2+
1108      &            ftbis(1)*welec*ees
1109      &            +ftbis(1)*wvdwpp*evdw
1110      &            +ftbis(1)*wtor*etors+
1111      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1112      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1113      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1114      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1115      &            ftbis(1)*wsccor*esccor
1116         else
1117             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1118      &      +wvdwpp*evdw1
1119      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1120      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1121      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1122      &      +ft(2)*wturn3*eello_turn3
1123      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1124      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1125      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1126             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1127      &            +ftprim(1)*wtor*etors+
1128      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1129      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1130      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1131      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1132      &            ftprim(1)*wsccor*esccor
1133             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1134      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1135      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1136      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1137      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1138      &            ftbis(1)*wsccor*esccor
1139       endif
1140 #else
1141       if (shield_mode.gt.0) then
1142             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1143      &      +ft(1)*welec*(ees+evdw1)
1144      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1145      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1146      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1147      &      +ft(2)*wturn3*eello_turn3
1148      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1149      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1150      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1151             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1152      &           +ftprim(1)*welec*(ees+evdw1)
1153      &           +ftprim(1)*wtor*etors+
1154      &            ftprim(1)*wscp*evdw2+
1155      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1156      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1157      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1158      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1159      &            ftprim(1)*wsccor*esccor
1160             ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1161      &           +ftbis(1)*wscp*evdw2
1162      &           +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1163      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1164      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1165      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1166      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1167      &            ftbis(1)*wsccor*esccor
1168        else
1169             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1170      &      +ft(1)*welec*(ees+evdw1)
1171      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1172      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1173      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1174      &      +ft(2)*wturn3*eello_turn3
1175      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1176      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1177      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1178             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1179      &           +ftprim(1)*wtor*etors+
1180      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1181      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1182      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1183      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1184      &            ftprim(1)*wsccor*esccor
1185             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1186      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1187      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1188      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1189      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1190      &            ftbis(1)*wsccor*esccor
1191
1192        endif
1193
1194 #endif
1195             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1196 #ifdef DEBUG
1197             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1198      &        " etot",etot," entfac",entfac(t),
1199      &        " weight",weight," ebis",ebis
1200 #endif
1201             etot=etot-temper*eprim
1202 #ifdef MPI
1203             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1204             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1205             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1206             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1207             do j=1,nQ+2
1208               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1209               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1210               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1211      &         +etot*q(j,t)*weight
1212             enddo
1213 #else
1214             sumW(k,iparm)=sumW(k,iparm)+weight
1215             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1216             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1217             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1218             do j=1,nQ+2
1219               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1220               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1221               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1222      &         +etot*q(j,t)*weight
1223             enddo
1224 #endif
1225           enddo
1226           indE = aint(potE(t,iparm)-aint(potEmin))
1227           if (indE.ge.0 .and. indE.le.maxinde) then
1228             if (indE.gt.upindE_p) upindE_p=indE
1229             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1230           endif
1231 #ifdef MPI
1232           do ib=1,nT_h(iparm)
1233             potEmin=potEmin_all(ib,iparm)
1234             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1235             hfin_p(ind,ib)=hfin_p(ind,ib)+
1236      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1237              if (rmsrgymap) then
1238                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1239                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1240                hrmsrgy_p(indrgy,indrms,ib)=
1241      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1242              endif
1243           enddo
1244 #else
1245           do ib=1,nT_h(iparm)
1246             potEmin=potEmin_all(ib,iparm)
1247             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1248             hfin(ind,ib)=hfin(ind,ib)+
1249      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1250              if (rmsrgymap) then
1251                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1252                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1253                hrmsrgy(indrgy,indrms,ib)=
1254      &           hrmsrgy(indrgy,indrms,ib)+expfac
1255              endif
1256           enddo
1257 #endif
1258         enddo ! t
1259         do ib=1,nT_h(iparm)
1260           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1261      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1262           if (rmsrgymap) then
1263           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1264      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1265      &       WHAM_COMM,IERROR)
1266           endif
1267         enddo
1268         call MPI_Reduce(upindE_p,upindE,1,
1269      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1270         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1271      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1272
1273         if (me1.eq.master) then
1274
1275         if (histout) then
1276
1277         write (iout,'(6x,$)')
1278         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1279      &   ib=1,nT_h(iparm))
1280         write (iout,*)
1281
1282         write (iout,'(/a)') 'Final histograms'
1283         if (histfile) then
1284           if (nslice.eq.1) then
1285             if (separate_parset) then
1286               write(licz3,"(bz,i3.3)") myparm
1287               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1288             else
1289               histname=prefix(:ilen(prefix))//'.hist'
1290             endif
1291           else
1292             if (separate_parset) then
1293               write(licz3,"(bz,i3.3)") myparm
1294               histname=prefix(:ilen(prefix))//'_par'//licz3//
1295      &         '_slice_'//licz2//'.hist'
1296             else
1297               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1298             endif
1299           endif
1300 #if defined(AIX) || defined(PGI)
1301           open (ihist,file=histname,position='append')
1302 #else
1303           open (ihist,file=histname,access='append')
1304 #endif
1305         endif
1306
1307         do t=0,tmax
1308           liczba=t
1309           sumH=0.0d0
1310           do ib=1,nT_h(iparm)
1311             sumH=sumH+hfin(t,ib)
1312           enddo
1313           if (sumH.gt.0.0d0) then
1314             do j=1,nQ
1315               jj = mod(liczba,nbin1)
1316               liczba=liczba/nbin1
1317               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1318               if (histfile)
1319      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1320             enddo
1321             do ib=1,nT_h(iparm)
1322               write (iout,'(e20.10,$)') hfin(t,ib)
1323               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1324             enddo
1325             write (iout,'(i5)') iparm
1326             if (histfile) write (ihist,'(i5)') iparm
1327           endif
1328         enddo
1329
1330         endif
1331
1332         if (entfile) then
1333           if (nslice.eq.1) then
1334             if (separate_parset) then
1335               write(licz3,"(bz,i3.3)") myparm
1336               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1337             else
1338               histname=prefix(:ilen(prefix))//'.ent'
1339             endif
1340           else
1341             if (separate_parset) then
1342               write(licz3,"(bz,i3.3)") myparm
1343               histname=prefix(:ilen(prefix))//'par_'//licz3//
1344      &           '_slice_'//licz2//'.ent'
1345             else
1346               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1347             endif
1348           endif
1349 #if defined(AIX) || defined(PGI)
1350           open (ihist,file=histname,position='append')
1351 #else
1352           open (ihist,file=histname,access='append')
1353 #endif
1354           write (ihist,'(a)') "# Microcanonical entropy"
1355           do i=0,upindE
1356             write (ihist,'(f8.0,$)') dint(potEmin)+i
1357             if (histE(i).gt.0.0e0) then
1358               write (ihist,'(f15.5,$)') dlog(histE(i))
1359             else
1360               write (ihist,'(f15.5,$)') 0.0d0
1361             endif
1362           enddo
1363           write (ihist,*)
1364           close(ihist)
1365         endif
1366         write (iout,*) "Microcanonical entropy"
1367         do i=0,upindE
1368           write (iout,'(f8.0,$)') dint(potEmin)+i
1369           if (histE(i).gt.0.0e0) then
1370             write (iout,'(f15.5,$)') dlog(histE(i))
1371           else
1372             write (iout,'(f15.5,$)') 0.0d0
1373           endif
1374           write (iout,*)
1375         enddo
1376         if (rmsrgymap) then
1377           if (nslice.eq.1) then
1378             if (separate_parset) then
1379               write(licz3,"(bz,i3.3)") myparm
1380               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1381             else
1382               histname=prefix(:ilen(prefix))//'.rmsrgy'
1383             endif
1384           else
1385             if (separate_parset) then
1386               write(licz3,"(bz,i3.3)") myparm
1387               histname=prefix(:ilen(prefix))//'_par'//licz3//
1388      &         '_slice_'//licz2//'.rmsrgy'
1389             else
1390              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1391             endif
1392           endif
1393 #if defined(AIX) || defined(PGI)
1394           open (ihist,file=histname,position='append')
1395 #else
1396           open (ihist,file=histname,access='append')
1397 #endif
1398           do i=0,nbin_rms
1399             do j=0,nbin_rgy
1400               write(ihist,'(2f8.2,$)')
1401      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1402               do ib=1,nT_h(iparm)
1403                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1404                   write(ihist,'(e14.5,$)')
1405      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1406      &              +potEmin
1407                 else
1408                   write(ihist,'(e14.5,$)') 1.0d6
1409                 endif
1410               enddo
1411               write (ihist,'(i2)') iparm
1412             enddo
1413           enddo
1414           close(ihist)
1415         endif
1416         endif
1417       enddo ! iparm
1418 #ifdef MPI
1419       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1420      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1421       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1422      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1423       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1424      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1425       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1426      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1427       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1428      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1429       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1430      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1431      &   WHAM_COMM,IERROR)
1432       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1433      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1434      &   WHAM_COMM,IERROR)
1435       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1436      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1437      &   WHAM_COMM,IERROR)
1438       if (me.eq.master) then
1439 #endif
1440       write (iout,'(/a)') 'Thermal characteristics of folding'
1441       if (nslice.eq.1) then
1442         nazwa=prefix
1443       else
1444         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1445       endif
1446       iln=ilen(nazwa)
1447       if (nparmset.eq.1 .and. .not.separate_parset) then
1448         nazwa=nazwa(:iln)//".thermal"
1449       else if (nparmset.eq.1 .and. separate_parset) then
1450         write(licz3,"(bz,i3.3)") myparm
1451         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1452       endif
1453       do iparm=1,nParmSet
1454       if (nparmset.gt.1) then
1455         write(licz3,"(bz,i3.3)") iparm
1456         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1457       endif
1458       open(34,file=nazwa)
1459       if (separate_parset) then
1460         write (iout,'(a,i3)') "Parameter set",myparm
1461       else
1462         write (iout,'(a,i3)') "Parameter set",iparm
1463       endif
1464       do i=0,NGridT
1465         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1466         if (betaT.ge.beta_h(1,iparm)) then
1467           potEmin=potEmin_all(1,iparm)
1468         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1469           potEmin=potEmin_all(nT_h(iparm),iparm)
1470         else
1471           do l=1,nT_h(iparm)-1
1472             if (betaT.le.beta_h(l,iparm) .and.
1473      &          betaT.gt.beta_h(l+1,iparm)) then
1474               potEmin=potEmin_all(l,iparm)
1475               exit
1476             endif
1477           enddo
1478         endif
1479
1480         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1481         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1482      &    sumW(i,iparm)
1483         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1484      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1485         do j=1,nQ+2
1486           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1487           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1488      &     -sumQ(j,i,iparm)**2
1489           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1490      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1491         enddo
1492         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1493      &   (startGridT+i*delta_T))+potEmin
1494         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1495      &   sumW(i,iparm),sumE(i,iparm)
1496         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1497         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1498      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1499         write (iout,*)
1500         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1501      &   sumW(i,iparm),sumE(i,iparm)
1502         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1503         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1504      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1505         write (34,*)
1506         call flush(34)
1507       enddo
1508       close(34)
1509       enddo
1510       if (histout) then
1511       do t=0,tmax
1512         if (hfin_ent(t).gt.0.0d0) then
1513           liczba=t
1514           jj = mod(liczba,nbin1)
1515           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1516      &     hfin_ent(t)
1517           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1518      &      dmin+(jj+0.5d0)*delta,
1519      &     hfin_ent(t)
1520         endif
1521       enddo
1522       if (histfile) close(ihist)
1523       endif
1524
1525 #ifdef ZSCORE
1526 ! Write data for zscore
1527       if (nslice.eq.1) then
1528         zscname=prefix(:ilen(prefix))//".zsc"
1529       else
1530         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1531       endif
1532 #if defined(AIX) || defined(PGI)
1533       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1534 #else
1535       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1536 #endif
1537       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1538       do iparm=1,nParmSet
1539         write (izsc,'("NT=",i1)') nT_h(iparm)
1540       do ib=1,nT_h(iparm)
1541         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1542      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1543         jj = min0(nR(ib,iparm),7)
1544         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1545         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1546         write (izsc,'("&")')
1547         if (nR(ib,iparm).gt.7) then
1548           do ii=8,nR(ib,iparm),9
1549             jj = min0(nR(ib,iparm),ii+8)
1550             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1551             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1552             write (izsc,'("&")')
1553           enddo
1554         endif
1555         write (izsc,'("FI=",$)')
1556         jj=min0(nR(ib,iparm),7)
1557         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1558         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1559         write (izsc,'("&")')
1560         if (nR(ib,iparm).gt.7) then
1561           do ii=8,nR(ib,iparm),9
1562             jj = min0(nR(ib,iparm),ii+8)
1563             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1564             if (jj.eq.nR(ib,iparm)) then
1565               write (izsc,*)
1566             else
1567               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1568               write (izsc,'(t80,"&")')
1569             endif
1570           enddo
1571         endif
1572         do i=1,nR(ib,iparm)
1573           write (izsc,'("KH=",$)')
1574           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1575           write (izsc,'(" Q0=",$)')
1576           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1577           write (izsc,*)
1578         enddo
1579       enddo
1580       enddo
1581       close(izsc)
1582 #endif
1583 #ifdef MPI
1584       endif
1585 #endif
1586
1587       return
1588 C#undef DEBUG
1589       end