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
  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
 329 #ifdef DEBUG
 330             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
 331      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 332      &       etors,etors_d,eello_turn3,eello_turn4,esccor
 333 #endif
 334
 335 #ifdef SPLITELE
 336             if (shield_mode.gt.0) then
 337             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 338      &      +ft(1)*welec*ees
 339      &      +ft(1)*wvdwpp*evdw1
 340      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 341      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 342      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 343      &      +ft(2)*wturn3*eello_turn3
 344      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 345      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 346      &      +wbond*estr+wliptran*eliptran
 347              else
 348             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 349      &      +wvdwpp*evdw1
 350      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 351      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 352      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 353      &      +ft(2)*wturn3*eello_turn3
 354      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 355      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 356      &      +wbond*estr+wliptran*eliptran
 357              endif
 358 #else
 359       if (shield_mode.gt.0) then
 360             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 361      &      +ft(1)*welec*(ees+evdw1)
 362      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 363      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 364      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 365      &      +ft(2)*wturn3*eello_turn3
 366      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 367      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 368      &      +wbond*estr+wliptran*eliptran
 369             else
 370             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 371      &      +ft(1)*welec*(ees+evdw1)
 372      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 373      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 374      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 375      &      +ft(2)*wturn3*eello_turn3
 376      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 377      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 378      &      +wbond*estr+wliptran*eliptran
 379            endif
 380
 381 #endif
 382 #ifdef DEBUG
 383             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 384      &        etot,potEmin
 385 #endif
 386 #ifdef DEBUG
 387             if (iparm.eq.1 .and. ib.eq.1) then
 388             write (iout,*)"Conformation",i
 389             energia(0)=etot
 390             do k=1,max_ene
 391               energia(k)=enetb(k,i,iparm)
 392             enddo
 393             call enerprint(energia(0),fT)
 394             endif
 395 #endif
 396             do kk=1,nR(ib,iparm)
 397               Econstr=0.0d0
 398               do j=1,nQ
 399                 dd = q(j,i)
 400                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 401      &           *(dd-q0(j,kk,ib,iparm))**2
 402               enddo
 403               v(i,kk,ib,iparm)=
 404      &          -beta_h(ib,iparm)*(etot+Econstr)
 405 #ifdef DEBUG
 406               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 407      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 408 #endif
 409             enddo ! kk
 410           enddo   ! ib
 411         enddo     ! iparm
 412       enddo       ! i
 413 ! Simple iteration to calculate free energies corresponding to all simulation
 414 ! runs.
 415       do iter=1,maxit
 416
 417 ! Compute new free-energy values corresponding to the righ-hand side of the
 418 ! equation and their derivatives.
 419         write (iout,*) "------------------------fi"
 420         entfac_min=1.0d10
 421
 422 #ifdef MPI
 423         do t=1,scount(me1)
 424 #else
 425         do t=1,ntot(islice)
 426 #endif
 427           vmax=-1.0d+20
 428           do i=1,nParmSet
 429             do k=1,nT_h(i)
 430               do l=1,nR(k,i)
 431                 vf=v(t,l,k,i)+f(l,k,i)
 432                 if (vf.gt.vmax) vmax=vf
 433               enddo
 434             enddo
 435           enddo
 436           denom=0.0d0
 437           do i=1,nParmSet
 438             do k=1,nT_h(i)
 439               do l=1,nR(k,i)
 440                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 441                 if (aux.gt.-200.0d0)
 442      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 443               enddo
 444             enddo
 445           enddo
 446           entfac(t)=-dlog(denom)-vmax
 447           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 448 #ifdef DEBUG
 449           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 450 #endif
 451         enddo
 452
 453         do iparm=1,nParmSet
 454           do iib=1,nT_h(iparm)
 455             do ii=1,nR(iib,iparm)
 456 #ifdef MPI
 457               fi_p_min(ii,iib,iparm)=-1.0d10
 458               do t=1,scount(me)
 459                 aux=v(t,ii,iib,iparm)+entfac(t)
 460                 if (aux.gt.fi_p_min(ii,iib,iparm))
 461      &                                   fi_p_min(ii,iib,iparm)=aux
 462               enddo
 463 #else
 464               do t=1,ntot(islice)
 465                 aux=v(t,ii,iib,iparm)+entfac(t)
 466                 if (aux.gt.fi_min(ii,iib,iparm))
 467      &                                     fi_min(ii,iib,iparm)=aux
 468               enddo
 469 #endif
 470             enddo ! ii
 471           enddo ! iib
 472         enddo ! iparm
 473 #ifdef MPI
 474 #ifdef DEBUG
 475         write (iout,*) "fi_min before AllReduce"
 476         do i=1,nParmSet
 477           do j=1,nT_h(i)
 478             write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
 479           enddo
 480         enddo
 481 #endif
 482         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
 483      &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 484 #ifdef DEBUG
 485         write (iout,*) "fi_min after AllReduce"
 486         do i=1,nParmSet
 487           do j=1,nT_h(i)
 488             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 489           enddo
 490         enddo
 491 #endif
 492 #endif
 493         do iparm=1,nParmSet
 494           do iib=1,nT_h(iparm)
 495             do ii=1,nR(iib,iparm)
 496 #ifdef MPI
 497               fi_p(ii,iib,iparm)=0.0d0
 498               do t=1,scount(me)
 499                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 500      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 501 #ifdef DEBUG
 502               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
 503      &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
 504      &         fi_p(ii,iib,iparm)
 505 #endif
 506               enddo
 507 #else
 508               fi(ii,iib,iparm)=0.0d0
 509               do t=1,ntot(islice)
 510                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 511      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 512               enddo
 513 #endif
 514             enddo ! ii
 515           enddo ! iib
 516         enddo ! iparm
 517 #ifdef MPI
 518 #ifdef DEBUG
 519         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 520         do iparm=1,nParmSet
 521           do ib=1,nT_h(nparmset)
 522             write (iout,*) "iparm",iparm," ib",ib
 523             write (iout,*) "beta=",beta_h(ib,iparm)
 524             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 525           enddo
 526         enddo
 527 #endif
 528         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 529      &   maxR*MaxT_h*nParmSet
 530         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 531      &   " WHAM_COMM",WHAM_COMM
 532         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 533      &   MPI_DOUBLE_PRECISION,
 534      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 535 #ifdef DEBUG
 536         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 537         do iparm=1,nParmSet
 538           write (iout,*) "iparm",iparm
 539           do ib=1,nT_h(iparm)
 540             write (iout,*) "beta=",beta_h(ib,iparm)
 541             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 542           enddo
 543         enddo
 544 #endif
 545         if (me1.eq.Master) then
 546 #endif
 547         avefi=0.0d0
 548         do iparm=1,nParmSet
 549           do ib=1,nT_h(iparm)
 550             do i=1,nR(ib,iparm)
 551               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 552               avefi=avefi+fi(i,ib,iparm)
 553             enddo
 554           enddo
 555         enddo
 556         avefi=avefi/nFi
 557         do iparm=1,nParmSet
 558           write (iout,*) "Parameter set",iparm
 559           do ib =1,nT_h(iparm)
 560             write (iout,*) "beta=",beta_h(ib,iparm)
 561             do i=1,nR(ib,iparm)
 562               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 563             enddo
 564             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 565             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 566           enddo
 567         enddo
 568
 569 ! Compute the norm of free-energy increments.
 570         finorm=0.0d0
 571         do iparm=1,nParmSet
 572           do ib=1,nT_h(iparm)
 573             do i=1,nR(ib,iparm)
 574               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 575               f(i,ib,iparm)=fi(i,ib,iparm)
 576             enddo
 577           enddo
 578         enddo
 579
 580         write (iout,*) 'Iteration',iter,' finorm',finorm
 581
 582 #ifdef MPI
 583         endif
 584         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 585      &   MPI_DOUBLE_PRECISION,Master,
 586      &   WHAM_COMM,IERROR)
 587         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 588      &   WHAM_COMM,IERROR)
 589 #endif
 590 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 591         if (finorm.lt.fimin) then
 592           write (iout,*) 'Iteration converged'
 593           goto 20
 594         endif
 595
 596       enddo ! iter
 597
 598    20 continue
 599
 600 #ifdef MPI
 601       do i=1,scount(me1)
 602 #else
 603       do i=1,ntot(islice)
 604 #endif
 605 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 606         do iparm=1,nParmSet
 607 #ifdef DEBUG
 608           write (iout,'(2i5,21f8.2)') i,iparm,
 609      &     (enetb(k,i,iparm),k=1,21)
 610 #endif
 611           call restore_parm(iparm)
 612 #ifdef DEBUG
 613           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 614      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 615      &      wtor_d,wsccor,wbond
 616 #endif
 617           do ib=1,nT_h(iparm)
 618             if (rescale_mode.eq.1) then
 619               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 620               quotl=1.0d0
 621               kfacl=1.0d0
 622               do l=1,5
 623                 quotl1=quotl
 624                 quotl=quotl*quot
 625                 kfacl=kfacl*kfac
 626                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 627               enddo
 628 #if defined(FUNCTH)
 629               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 630               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 631 #elif defined(FUNCT)
 632               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 633 #else
 634               ft(6)=1.0d0
 635 #endif
 636             else if (rescale_mode.eq.2) then
 637               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 638               quotl=1.0d0
 639               do l=1,5
 640                 quotl=quotl*quot
 641                 fT(l)=1.12692801104297249644d0/
 642      &             dlog(dexp(quotl)+dexp(-quotl))
 643               enddo
 644 #if defined(FUNCTH)
 645               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 646               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 647 #elif defined(FUNCT)
 648               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 649 #else
 650               ft(6)=1.0d0
 651 #endif
 652 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 653             else if (rescale_mode.eq.0) then
 654               do l=1,6
 655                 fT(l)=1.0d0
 656               enddo
 657             else
 658               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 659      &          rescale_mode
 660               call flush(iout)
 661               return1
 662             endif
 663           evdw=enetb(21,i,iparm)
 664           evdw_t=enetb(1,i,iparm)
 665 #ifdef SCP14
 666           evdw2_14=enetb(17,i,iparm)
 667           evdw2=enetb(2,i,iparm)+evdw2_14
 668 #else
 669           evdw2=enetb(2,i,iparm)
 670           evdw2_14=0.0d0
 671 #endif
 672 #ifdef SPLITELE
 673           ees=enetb(3,i,iparm)
 674           evdw1=enetb(16,i,iparm)
 675 #else
 676           ees=enetb(3,i,iparm)
 677           evdw1=0.0d0
 678 #endif
 679           ecorr=enetb(4,i,iparm)
 680           ecorr5=enetb(5,i,iparm)
 681           ecorr6=enetb(6,i,iparm)
 682           eel_loc=enetb(7,i,iparm)
 683           eello_turn3=enetb(8,i,iparm)
 684           eello_turn4=enetb(9,i,iparm)
 685           eturn6=enetb(10,i,iparm)
 686           ebe=enetb(11,i,iparm)
 687           escloc=enetb(12,i,iparm)
 688           etors=enetb(13,i,iparm)
 689           etors_d=enetb(14,i,iparm)
 690           ehpb=enetb(15,i,iparm)
 691           estr=enetb(18,i,iparm)
 692           esccor=enetb(19,i,iparm)
 693           edihcnstr=enetb(20,i,iparm)
 694 C          edihcnstr=0.0d0
 695           eliptran=enetb(22,i,iparm)
 696 #ifdef SPLITELE
 697       if (shield_mode.gt.0) then
 698             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 699      &      +ft(1)*welec*ees
 700      &      +ft(1)*wvdwpp*evdw1
 701      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 702      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 703      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 704      &      +ft(2)*wturn3*eello_turn3
 705      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 706      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 707      &      +wbond*estr+wliptran*eliptran
 708         else
 709             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 710      &      +wvdwpp*evdw1
 711      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 712      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 713      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 714      &      +ft(2)*wturn3*eello_turn3
 715      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 716      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 717      &      +wbond*estr+wliptran*eliptran
 718       endif
 719 #else
 720       if (shield_mode.gt.0) then
 721             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 722      &      +ft(1)*welec*(ees+evdw1)
 723      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 724      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 725      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 726      &      +ft(2)*wturn3*eello_turn3
 727      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 728      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 729      &      +wbond*estr+wliptran*eliptran
 730        else
 731             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 732      &      +ft(1)*welec*(ees+evdw1)
 733      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 734      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 735      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 736      &      +ft(2)*wturn3*eello_turn3
 737      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 738      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 739      &      +wbond*estr+wliptran*eliptran
 740        endif
 741
 742 #endif
 743             etot=etot-entfac(i)/beta_h(ib,iparm)
 744             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 745
 746           enddo   ! ib
 747         enddo     ! iparm
 748       enddo       ! i
 749
 750
 751 #ifdef DEBUG
 752       write (iout,*) "The potEmin array before reduction"
 753       do i=1,nParmSet
 754         write (iout,*) "Parameter set",i
 755         do j=1,nT_h(i)
 756           write (iout,*) j,PotEmin_all(j,i)
 757         enddo
 758       enddo
 759       write (iout,*) "potEmin_min",potEmin_min
 760 #endif
 761 #ifdef MPI
 762 C Determine the minimum energes for all parameter sets and temperatures
 763       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
 764      &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 765       do i=1,nParmSet
 766         do j=1,nT_h(i)
 767           potEmin_all(j,i)=potEmin_t_all(j,i)
 768         enddo
 769       enddo
 770 #endif
 771       potEmin_min=potEmin_all(1,1)
 772       do i=1,nParmSet
 773         do j=1,nT_h(i)
 774           if (potEmin_all(j,i).lt.potEmin_min)
 775      &             potEmin_min=potEmin_all(j,i)
 776         enddo
 777       enddo
 778 #ifdef DEBUG
 779       write (iout,*) "The potEmin array"
 780       do i=1,nParmSet
 781         write (iout,*) "Parameter set",i
 782         do j=1,nT_h(i)
 783           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
 784      &      PotEmin_all(j,i)
 785         enddo
 786       enddo
 787       write (iout,*) "potEmin_min",potEmin_min
 788 #endif
 789
 790
 791 ! Now, put together the histograms from all simulations, in order to get the
 792 ! unbiased total histogram.
 793 #ifdef MPI
 794       do t=0,tmax
 795         hfin_ent_p(t)=0.0d0
 796       enddo
 797 #else
 798       do t=0,tmax
 799         hfin_ent(t)=0.0d0
 800       enddo
 801 #endif
 802       write (iout,*) "--------------hist"
 803 #ifdef MPI
 804       do iparm=1,nParmSet
 805         do i=0,nGridT
 806           sumW_p(i,iparm)=0.0d0
 807           sumE_p(i,iparm)=0.0d0
 808           sumEbis_p(i,iparm)=0.0d0
 809           sumEsq_p(i,iparm)=0.0d0
 810           do j=1,nQ+2
 811             sumQ_p(j,i,iparm)=0.0d0
 812             sumQsq_p(j,i,iparm)=0.0d0
 813             sumEQ_p(j,i,iparm)=0.0d0
 814           enddo
 815         enddo
 816       enddo
 817       upindE_p=0
 818 #else
 819       do iparm=1,nParmSet
 820         do i=0,nGridT
 821           sumW(i,iparm)=0.0d0
 822           sumE(i,iparm)=0.0d0
 823           sumEbis(i,iparm)=0.0d0
 824           sumEsq(i,iparm)=0.0d0
 825           do j=1,nQ+2
 826             sumQ(j,i,iparm)=0.0d0
 827             sumQsq(j,i,iparm)=0.0d0
 828             sumEQ(j,i,iparm)=0.0d0
 829           enddo
 830         enddo
 831       enddo
 832       upindE=0
 833 #endif
 834 c 8/26/05 entropy distribution
 835 #ifdef MPI
 836       entmin_p=1.0d10
 837       entmax_p=-1.0d10
 838       do t=1,scount(me1)
 839 c        ent=-dlog(entfac(t))
 840         ent=entfac(t)
 841         if (ent.lt.entmin_p) entmin_p=ent
 842         if (ent.gt.entmax_p) entmax_p=ent
 843       enddo
 844       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 845       call flush(iout)
 846       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 847      &  WHAM_COMM,IERROR)
 848       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 849      &  WHAM_COMM,IERROR)
 850       ientmax=entmax-entmin
 851       if (ientmax.gt.2000) ientmax=2000
 852       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 853       call flush(iout)
 854       do t=1,scount(me1)
 855 c        ient=-dlog(entfac(t))-entmin
 856         ient=entfac(t)-entmin
 857         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 858       enddo
 859       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 860      &  MPI_SUM,WHAM_COMM,IERROR)
 861       if (me1.eq.Master) then
 862         write (iout,*) "Entropy histogram"
 863         do i=0,ientmax
 864           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 865         enddo
 866       endif
 867 #else
 868       entmin=1.0d10
 869       entmax=-1.0d10
 870       do t=1,ntot(islice)
 871         ent=entfac(t)
 872         if (ent.lt.entmin) entmin=ent
 873         if (ent.gt.entmax) entmax=ent
 874       enddo
 875       ientmax=-dlog(entmax)-entmin
 876       if (ientmax.gt.2000) ientmax=2000
 877       do t=1,ntot(islice)
 878         ient=entfac(t)-entmin
 879         if (ient.le.2000) histent(ient)=histent(ient)+1
 880       enddo
 881       write (iout,*) "Entropy histogram"
 882       do i=0,ientmax
 883         write(iout,'(2f15.4)') entmin+i,histent(i)
 884       enddo
 885 #endif
 886
 887 #ifdef MPI
 888 c      write (iout,*) "me1",me1," scount",scount(me1)
 889
 890       do iparm=1,nParmSet
 891
 892 #ifdef MPI
 893         do ib=1,nT_h(iparm)
 894           do t=0,tmax
 895             hfin_p(t,ib)=0.0d0
 896           enddo
 897         enddo
 898         do i=1,maxindE
 899           histE_p(i)=0.0d0
 900         enddo
 901 #else
 902         do ib=1,nT_h(iparm)
 903           do t=0,tmax
 904             hfin(t,ib)=0.0d0
 905           enddo
 906         enddo
 907         do i=1,maxindE
 908           histE(i)=0.0d0
 909         enddo
 910 #endif
 911         do ib=1,nT_h(iparm)
 912           do i=0,MaxBinRms
 913             do j=0,MaxBinRgy
 914               hrmsrgy(j,i,ib)=0.0d0
 915 #ifdef MPI
 916               hrmsrgy_p(j,i,ib)=0.0d0
 917 #endif
 918             enddo
 919           enddo
 920         enddo
 921
 922         do t=1,scount(me1)
 923 #else
 924         do t=1,ntot(islice)
 925 #endif
 926           ind = ind_point(t)
 927 #ifdef MPI
 928           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 929 #else
 930           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 931 #endif
 932 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 933           call restore_parm(iparm)
 934           evdw=enetb(21,t,iparm)
 935           evdw_t=enetb(1,t,iparm)
 936 #ifdef SCP14
 937           evdw2_14=enetb(17,t,iparm)
 938           evdw2=enetb(2,t,iparm)+evdw2_14
 939 #else
 940           evdw2=enetb(2,t,iparm)
 941           evdw2_14=0.0d0
 942 #endif
 943 #ifdef SPLITELE
 944           ees=enetb(3,t,iparm)
 945           evdw1=enetb(16,t,iparm)
 946 #else
 947           ees=enetb(3,t,iparm)
 948           evdw1=0.0d0
 949 #endif
 950           ecorr=enetb(4,t,iparm)
 951           ecorr5=enetb(5,t,iparm)
 952           ecorr6=enetb(6,t,iparm)
 953           eel_loc=enetb(7,t,iparm)
 954           eello_turn3=enetb(8,t,iparm)
 955           eello_turn4=enetb(9,t,iparm)
 956           eturn6=enetb(10,t,iparm)
 957           ebe=enetb(11,t,iparm)
 958           escloc=enetb(12,t,iparm)
 959           etors=enetb(13,t,iparm)
 960           etors_d=enetb(14,t,iparm)
 961           ehpb=enetb(15,t,iparm)
 962           estr=enetb(18,t,iparm)
 963           esccor=enetb(19,t,iparm)
 964           edihcnstr=enetb(20,t,iparm)
 965 C          edihcnstr=0.0d0
 966           eliptran=enetb(22,i,iparm)
 967
 968           do k=0,nGridT
 969             betaT=startGridT+k*delta_T
 970             temper=betaT
 971 c            fT=T0/betaT
 972 c            ft=2*T0/(T0+betaT)
 973             if (rescale_mode.eq.1) then
 974               quot=betaT/T0
 975               quotl=1.0d0
 976               kfacl=1.0d0
 977               do l=1,5
 978                 quotl1=quotl
 979                 quotl=quotl*quot
 980                 kfacl=kfacl*kfac
 981                 denom=kfacl-1.0d0+quotl
 982                 fT(l)=kfacl/denom
 983                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 984                 ftbis(l)=l*kfacl*quotl1*
 985      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 986               enddo
 987 #if defined(FUNCTH)
 988               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 989      &                  320.0d0
 990               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 991              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 992      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 993 #elif defined(FUNCT)
 994               fT(6)=betaT/T0
 995               ftprim(6)=1.0d0/T0
 996               ftbis(6)=0.0d0
 997 #else
 998               fT(6)=1.0d0
 999               ftprim(6)=0.0d0
1000               ftbis(6)=0.0d0
1001 #endif
1002             else if (rescale_mode.eq.2) then
1003               quot=betaT/T0
1004               quotl=1.0d0
1005               do l=1,5
1006                 quotl1=quotl
1007                 quotl=quotl*quot
1008                 eplus=dexp(quotl)
1009                 eminus=dexp(-quotl)
1010                 logfac=1.0d0/dlog(eplus+eminus)
1011                 tanhT=(eplus-eminus)/(eplus+eminus)
1012                 fT(l)=1.12692801104297249644d0*logfac
1013                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1014                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1015      &          2*l*quotl1/T0*logfac*
1016      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1017      &          +ftprim(l)*tanhT)
1018               enddo
1019 #if defined(FUNCTH)
1020               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1021      &                 320.0d0
1022               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1023              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1024      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1025 #elif defined(FUNCT)
1026               fT(6)=betaT/T0
1027               ftprim(6)=1.0d0/T0
1028               ftbis(6)=0.0d0
1029 #else
1030               fT(6)=1.0d0
1031               ftprim(6)=0.0d0
1032               ftbis(6)=0.0d0
1033 #endif
1034             else if (rescale_mode.eq.0) then
1035               do l=1,5
1036                 fT(l)=1.0d0
1037                 ftprim(l)=0.0d0
1038               enddo
1039             else
1040               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1041      &          rescale_mode
1042               call flush(iout)
1043               return1
1044             endif
1045 c            write (iout,*) "ftprim",ftprim
1046 c            write (iout,*) "ftbis",ftbis
1047             betaT=1.0d0/(1.987D-3*betaT)
1048             if (betaT.ge.beta_h(1,iparm)) then
1049               potEmin=potEmin_all(1,iparm)+
1050      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1051      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1052      &         (1.0/betaT-1.0/beta_h(1,iparm))
1053 #ifdef DEBUG
1054               write(iout,*) "first",temper,potEmin
1055 #endif
1056             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1057               potEmin=potEmin_all(nT_h(iparm),iparm)+
1058      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1059      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1060      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1061 #ifdef DEBUG
1062               write (iout,*) "last",temper,potEmin
1063 #endif
1064             else
1065               do l=1,nT_h(iparm)-1
1066                 if (betaT.le.beta_h(l,iparm) .and.
1067      &              betaT.gt.beta_h(l+1,iparm)) then
1068                   potEmin=potEmin_all(l,iparm)
1069 #ifdef DEBUG
1070                   write (iout,*) "l",l,
1071      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1072      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1073 #endif
1074                   exit
1075                 endif
1076               enddo
1077             endif
1078 #ifdef SPLITELE
1079       if (shield_mode.gt.0) then
1080             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1081      &      +ft(1)*welec*ees
1082      &      +ft(1)*wvdwpp*evdw1
1083      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1084      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1085      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1086      &      +ft(2)*wturn3*eello_turn3
1087      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1088      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1089      &      +wbond*estr+wliptran*eliptran
1090             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1091 C     &            +ftprim(6)*evdw_t
1092      &            +ftprim(1)*wscp*evdw2
1093      &            +ftprim(1)*welec*ees
1094      &            +ftprim(1)*wvdwpp*evdw1
1095      &            +ftprim(1)*wtor*etors+
1096      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1097      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1098      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1099      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1100      &            ftprim(1)*wsccor*esccor
1101             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1102      &            +ftbis(1)*wscp*evdw2+
1103      &            ftbis(1)*welec*ees
1104      &            +ftbis(1)*wvdwpp*evdw
1105      &            +ftbis(1)*wtor*etors+
1106      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1107      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1108      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1109      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1110      &            ftbis(1)*wsccor*esccor
1111         else
1112             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1113      &      +wvdwpp*evdw1
1114      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1115      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1116      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1117      &      +ft(2)*wturn3*eello_turn3
1118      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1119      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1120      &      +wbond*estr+wliptran*eliptran
1121             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1122      &            +ftprim(1)*wtor*etors+
1123      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1124      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1125      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1126      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1127      &            ftprim(1)*wsccor*esccor
1128             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1129      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1130      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1131      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1132      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1133      &            ftbis(1)*wsccor*esccor
1134       endif
1135 #else
1136       if (shield_mode.gt.0) then
1137             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1138      &      +ft(1)*welec*(ees+evdw1)
1139      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1140      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1141      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1142      &      +ft(2)*wturn3*eello_turn3
1143      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1144      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1145      &      +wbond*estr+wliptran*eliptran
1146             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1147      &           +ftprim(1)*welec*(ees+evdw1)
1148      &           +ftprim(1)*wtor*etors+
1149      &            ftprim(1)*wscp*evdw2+
1150      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1151      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1152      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1153      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1154      &            ftprim(1)*wsccor*esccor
1155             ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1156      &           +ftbis(1)*wscp*evdw2
1157      &           +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1158      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1159      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1160      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1161      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1162      &            ftbis(1)*wsccor*esccor
1163        else
1164             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1165      &      +ft(1)*welec*(ees+evdw1)
1166      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1167      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1168      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1169      &      +ft(2)*wturn3*eello_turn3
1170      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1171      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1172      &      +wbond*estr+wliptran*eliptran
1173             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1174      &           +ftprim(1)*wtor*etors+
1175      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1176      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1177      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1178      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1179      &            ftprim(1)*wsccor*esccor
1180             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1181      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1182      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1183      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1184      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1185      &            ftbis(1)*wsccor*esccor
1186
1187        endif
1188
1189 #endif
1190             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1191 #ifdef DEBUG
1192             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1193      &        " etot",etot," entfac",entfac(t),
1194      &        " weight",weight," ebis",ebis
1195 #endif
1196             etot=etot-temper*eprim
1197 #ifdef MPI
1198             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1199             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1200             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1201             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1202             do j=1,nQ+2
1203               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1204               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1205               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1206      &         +etot*q(j,t)*weight
1207             enddo
1208 #else
1209             sumW(k,iparm)=sumW(k,iparm)+weight
1210             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1211             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1212             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1213             do j=1,nQ+2
1214               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1215               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1216               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1217      &         +etot*q(j,t)*weight
1218             enddo
1219 #endif
1220           enddo
1221           indE = aint(potE(t,iparm)-aint(potEmin))
1222           if (indE.ge.0 .and. indE.le.maxinde) then
1223             if (indE.gt.upindE_p) upindE_p=indE
1224             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1225           endif
1226 #ifdef MPI
1227           do ib=1,nT_h(iparm)
1228             potEmin=potEmin_all(ib,iparm)
1229             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1230             hfin_p(ind,ib)=hfin_p(ind,ib)+
1231      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1232              if (rmsrgymap) then
1233                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1234                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1235                hrmsrgy_p(indrgy,indrms,ib)=
1236      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1237              endif
1238           enddo
1239 #else
1240           do ib=1,nT_h(iparm)
1241             potEmin=potEmin_all(ib,iparm)
1242             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1243             hfin(ind,ib)=hfin(ind,ib)+
1244      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1245              if (rmsrgymap) then
1246                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1247                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1248                hrmsrgy(indrgy,indrms,ib)=
1249      &           hrmsrgy(indrgy,indrms,ib)+expfac
1250              endif
1251           enddo
1252 #endif
1253         enddo ! t
1254         do ib=1,nT_h(iparm)
1255           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1256      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1257           if (rmsrgymap) then
1258           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1259      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1260      &       WHAM_COMM,IERROR)
1261           endif
1262         enddo
1263         call MPI_Reduce(upindE_p,upindE,1,
1264      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1265         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1266      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1267
1268         if (me1.eq.master) then
1269
1270         if (histout) then
1271
1272         write (iout,'(6x,$)')
1273         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1274      &   ib=1,nT_h(iparm))
1275         write (iout,*)
1276
1277         write (iout,'(/a)') 'Final histograms'
1278         if (histfile) then
1279           if (nslice.eq.1) then
1280             if (separate_parset) then
1281               write(licz3,"(bz,i3.3)") myparm
1282               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1283             else
1284               histname=prefix(:ilen(prefix))//'.hist'
1285             endif
1286           else
1287             if (separate_parset) then
1288               write(licz3,"(bz,i3.3)") myparm
1289               histname=prefix(:ilen(prefix))//'_par'//licz3//
1290      &         '_slice_'//licz2//'.hist'
1291             else
1292               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1293             endif
1294           endif
1295 #if defined(AIX) || defined(PGI)
1296           open (ihist,file=histname,position='append')
1297 #else
1298           open (ihist,file=histname,access='append')
1299 #endif
1300         endif
1301
1302         do t=0,tmax
1303           liczba=t
1304           sumH=0.0d0
1305           do ib=1,nT_h(iparm)
1306             sumH=sumH+hfin(t,ib)
1307           enddo
1308           if (sumH.gt.0.0d0) then
1309             do j=1,nQ
1310               jj = mod(liczba,nbin1)
1311               liczba=liczba/nbin1
1312               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1313               if (histfile)
1314      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1315             enddo
1316             do ib=1,nT_h(iparm)
1317               write (iout,'(e20.10,$)') hfin(t,ib)
1318               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1319             enddo
1320             write (iout,'(i5)') iparm
1321             if (histfile) write (ihist,'(i5)') iparm
1322           endif
1323         enddo
1324
1325         endif
1326
1327         if (entfile) then
1328           if (nslice.eq.1) then
1329             if (separate_parset) then
1330               write(licz3,"(bz,i3.3)") myparm
1331               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1332             else
1333               histname=prefix(:ilen(prefix))//'.ent'
1334             endif
1335           else
1336             if (separate_parset) then
1337               write(licz3,"(bz,i3.3)") myparm
1338               histname=prefix(:ilen(prefix))//'par_'//licz3//
1339      &           '_slice_'//licz2//'.ent'
1340             else
1341               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1342             endif
1343           endif
1344 #if defined(AIX) || defined(PGI)
1345           open (ihist,file=histname,position='append')
1346 #else
1347           open (ihist,file=histname,access='append')
1348 #endif
1349           write (ihist,'(a)') "# Microcanonical entropy"
1350           do i=0,upindE
1351             write (ihist,'(f8.0,$)') dint(potEmin)+i
1352             if (histE(i).gt.0.0e0) then
1353               write (ihist,'(f15.5,$)') dlog(histE(i))
1354             else
1355               write (ihist,'(f15.5,$)') 0.0d0
1356             endif
1357           enddo
1358           write (ihist,*)
1359           close(ihist)
1360         endif
1361         write (iout,*) "Microcanonical entropy"
1362         do i=0,upindE
1363           write (iout,'(f8.0,$)') dint(potEmin)+i
1364           if (histE(i).gt.0.0e0) then
1365             write (iout,'(f15.5,$)') dlog(histE(i))
1366           else
1367             write (iout,'(f15.5,$)') 0.0d0
1368           endif
1369           write (iout,*)
1370         enddo
1371         if (rmsrgymap) then
1372           if (nslice.eq.1) then
1373             if (separate_parset) then
1374               write(licz3,"(bz,i3.3)") myparm
1375               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1376             else
1377               histname=prefix(:ilen(prefix))//'.rmsrgy'
1378             endif
1379           else
1380             if (separate_parset) then
1381               write(licz3,"(bz,i3.3)") myparm
1382               histname=prefix(:ilen(prefix))//'_par'//licz3//
1383      &         '_slice_'//licz2//'.rmsrgy'
1384             else
1385              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1386             endif
1387           endif
1388 #if defined(AIX) || defined(PGI)
1389           open (ihist,file=histname,position='append')
1390 #else
1391           open (ihist,file=histname,access='append')
1392 #endif
1393           do i=0,nbin_rms
1394             do j=0,nbin_rgy
1395               write(ihist,'(2f8.2,$)')
1396      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1397               do ib=1,nT_h(iparm)
1398                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1399                   write(ihist,'(e14.5,$)')
1400      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1401      &              +potEmin
1402                 else
1403                   write(ihist,'(e14.5,$)') 1.0d6
1404                 endif
1405               enddo
1406               write (ihist,'(i2)') iparm
1407             enddo
1408           enddo
1409           close(ihist)
1410         endif
1411         endif
1412       enddo ! iparm
1413 #ifdef MPI
1414       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1415      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1416       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1417      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1418       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1419      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1420       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1421      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1422       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1423      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1424       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1425      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1426      &   WHAM_COMM,IERROR)
1427       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1428      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1429      &   WHAM_COMM,IERROR)
1430       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1431      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1432      &   WHAM_COMM,IERROR)
1433       if (me.eq.master) then
1434 #endif
1435       write (iout,'(/a)') 'Thermal characteristics of folding'
1436       if (nslice.eq.1) then
1437         nazwa=prefix
1438       else
1439         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1440       endif
1441       iln=ilen(nazwa)
1442       if (nparmset.eq.1 .and. .not.separate_parset) then
1443         nazwa=nazwa(:iln)//".thermal"
1444       else if (nparmset.eq.1 .and. separate_parset) then
1445         write(licz3,"(bz,i3.3)") myparm
1446         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1447       endif
1448       do iparm=1,nParmSet
1449       if (nparmset.gt.1) then
1450         write(licz3,"(bz,i3.3)") iparm
1451         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1452       endif
1453       open(34,file=nazwa)
1454       if (separate_parset) then
1455         write (iout,'(a,i3)') "Parameter set",myparm
1456       else
1457         write (iout,'(a,i3)') "Parameter set",iparm
1458       endif
1459       do i=0,NGridT
1460         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1461         if (betaT.ge.beta_h(1,iparm)) then
1462           potEmin=potEmin_all(1,iparm)
1463         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1464           potEmin=potEmin_all(nT_h(iparm),iparm)
1465         else
1466           do l=1,nT_h(iparm)-1
1467             if (betaT.le.beta_h(l,iparm) .and.
1468      &          betaT.gt.beta_h(l+1,iparm)) then
1469               potEmin=potEmin_all(l,iparm)
1470               exit
1471             endif
1472           enddo
1473         endif
1474
1475         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1476         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1477      &    sumW(i,iparm)
1478         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1479      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1480         do j=1,nQ+2
1481           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1482           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1483      &     -sumQ(j,i,iparm)**2
1484           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1485      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1486         enddo
1487         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1488      &   (startGridT+i*delta_T))+potEmin
1489         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1490      &   sumW(i,iparm),sumE(i,iparm)
1491         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1492         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1493      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1494         write (iout,*)
1495         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1496      &   sumW(i,iparm),sumE(i,iparm)
1497         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1498         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1499      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1500         write (34,*)
1501         call flush(34)
1502       enddo
1503       close(34)
1504       enddo
1505       if (histout) then
1506       do t=0,tmax
1507         if (hfin_ent(t).gt.0.0d0) then
1508           liczba=t
1509           jj = mod(liczba,nbin1)
1510           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1511      &     hfin_ent(t)
1512           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1513      &      dmin+(jj+0.5d0)*delta,
1514      &     hfin_ent(t)
1515         endif
1516       enddo
1517       if (histfile) close(ihist)
1518       endif
1519
1520 #ifdef ZSCORE
1521 ! Write data for zscore
1522       if (nslice.eq.1) then
1523         zscname=prefix(:ilen(prefix))//".zsc"
1524       else
1525         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1526       endif
1527 #if defined(AIX) || defined(PGI)
1528       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1529 #else
1530       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1531 #endif
1532       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1533       do iparm=1,nParmSet
1534         write (izsc,'("NT=",i1)') nT_h(iparm)
1535       do ib=1,nT_h(iparm)
1536         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1537      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1538         jj = min0(nR(ib,iparm),7)
1539         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1540         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1541         write (izsc,'("&")')
1542         if (nR(ib,iparm).gt.7) then
1543           do ii=8,nR(ib,iparm),9
1544             jj = min0(nR(ib,iparm),ii+8)
1545             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1546             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1547             write (izsc,'("&")')
1548           enddo
1549         endif
1550         write (izsc,'("FI=",$)')
1551         jj=min0(nR(ib,iparm),7)
1552         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1553         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1554         write (izsc,'("&")')
1555         if (nR(ib,iparm).gt.7) then
1556           do ii=8,nR(ib,iparm),9
1557             jj = min0(nR(ib,iparm),ii+8)
1558             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1559             if (jj.eq.nR(ib,iparm)) then
1560               write (izsc,*)
1561             else
1562               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1563               write (izsc,'(t80,"&")')
1564             endif
1565           enddo
1566         endif
1567         do i=1,nR(ib,iparm)
1568           write (izsc,'("KH=",$)')
1569           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1570           write (izsc,'(" Q0=",$)')
1571           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1572           write (izsc,*)
1573         enddo
1574       enddo
1575       enddo
1576       close(izsc)
1577 #endif
1578 #ifdef MPI
1579       endif
1580 #endif
1581
1582       return
1583 C#undef DEBUG
1584       end