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       integer MaxPoint,MaxPointProc
  40       parameter (MaxPoint=MaxStr,
  41      & MaxPointProc=MaxStr_Proc)
  42       double precision finorm_max,potfac,entmin,entmax,expfac,vf
  43       parameter (finorm_max=1.0d0)
  44       integer islice
  45       integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  46       integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
  47      &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
  48       integer htot(0:MaxHdim),histent(0:2000)
  49       double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  50       double precision energia(0:max_ene)
  51 #ifdef MPI
  52       integer tmax_t,upindE_p
  53       double precision fi_p(MaxR,MaxT_h,Max_Parm)
  54       double precision sumW_p(0:nGridT,Max_Parm),
  55      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
  56      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
  57      & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
  58      & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
  59      & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
  60      & sumEbis_p(0:nGridT,Max_Parm)
  61       double precision hfin_p(0:MaxHdim,maxT_h),
  62      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
  63      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  64       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  65       double precision potEmin_t,entmin_p,entmax_p
  66       integer histent_p(0:2000)
  67       logical lprint /.true./
  68 #endif
  69       double precision delta_T /1.0d0/
  70       double precision rgymin,rmsmin,rgymax,rmsmax
  71       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
  72      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
  73      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
  74      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
  75      & weight,econstr
  76       double precision fi(MaxR,maxT_h,Max_Parm),
  77      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
  78      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
  79      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
  80      & potEmin,ent,
  81      & hfin_ent(0:MaxHdim),vmax,aux
  82       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  83      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
  84      &  eplus,eminus,logfac,tanhT,tt
  85       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  86      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  87      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  88      &  eliptran
  89
  90       integer ind_point(maxpoint),upindE,indE
  91       character*16 plik
  92       character*1 licz1
  93       character*2 licz2
  94       character*3 licz3
  95       character*128 nazwa
  96       integer ilen
  97       external ilen
  98
  99       write(licz2,'(bz,i2.2)') islice
 100       nbin1 = 1.0d0/delta
 101       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
 102      &  i2/80(1h-)//)') islice
 103       write (iout,*) "delta",delta," nbin1",nbin1
 104       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 105       call flush(iout)
 106       dmin=0.0d0
 107       tmax=0
 108       potEmin=1.0d10
 109       rgymin=1.0d10
 110       rmsmin=1.0d10
 111       rgymax=0.0d0
 112       rmsmax=0.0d0
 113       do t=0,MaxN
 114         htot(t)=0
 115       enddo
 116 #ifdef MPI
 117       do i=1,scount(me1)
 118 #else
 119       do i=1,ntot(islice)
 120 #endif
 121         do j=1,nParmSet
 122           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 123         enddo
 124         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 125         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 126         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 127         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 128         ind_point(i)=0
 129         do j=nQ,1,-1
 130           ind=(q(j,i)-dmin+1.0d-8)/delta
 131           if (j.eq.1) then
 132             ind_point(i)=ind_point(i)+ind
 133           else
 134             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 135           endif
 136 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 137 c     &      ind_point(i)
 138           call flush(iout)
 139           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 140             write (iout,*) "Error - index exceeds range for point",i,
 141      &      " q=",q(j,i)," ind",ind_point(i)
 142 #ifdef MPI
 143             write (iout,*) "Processor",me1
 144             call flush(iout)
 145             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 146 #endif
 147             stop
 148           endif
 149         enddo ! j
 150         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 151         htot(ind_point(i))=htot(ind_point(i))+1
 152 #ifdef DEBUG
 153         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
 154      &   " htot",htot(ind_point(i))
 155         call flush(iout)
 156 #endif
 157       enddo ! i
 158       call flush(iout)
 159
 160       nbin=nbin1**nQ-1
 161       write (iout,'(a)') "Numbers of counts in Q bins"
 162       do t=0,tmax
 163         if (htot(t).gt.0) then
 164         write (iout,'(i15,$)') t
 165         liczba=t
 166         do j=1,nQ
 167           jj = mod(liczba,nbin1)
 168           liczba=liczba/nbin1
 169           write (iout,'(i5,$)') jj
 170         enddo
 171         write (iout,'(i8)') htot(t)
 172         endif
 173       enddo
 174       do iparm=1,nParmSet
 175       write (iout,'(a,i3)') "Number of data points for parameter set",
 176      & iparm
 177       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
 178      & ib=1,nT_h(iparm))
 179       write (iout,'(i8)') stot(islice)
 180       write (iout,'(a)')
 181       enddo
 182       call flush(iout)
 183
 184 #ifdef MPI
 185       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
 186      &  WHAM_COMM,IERROR)
 187       tmax=tmax_t
 188       call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
 189      &  MPI_MIN,WHAM_COMM,IERROR)
 190       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
 191      &  MPI_MIN,WHAM_COMM,IERROR)
 192       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
 193      &  MPI_MAX,WHAM_COMM,IERROR)
 194       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
 195      &  MPI_MIN,WHAM_COMM,IERROR)
 196       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
 197      &  MPI_MAX,WHAM_COMM,IERROR)
 198       potEmin=potEmin_t/2
 199       rgymin=rgymin_t
 200       rgymax=rgymax_t
 201       rmsmin=rmsmin_t
 202       rmsmax=rmsmax_t
 203       write (iout,*) "potEmin",potEmin
 204 #endif
 205       rmsmin=deltrms*dint(rmsmin/deltrms)
 206       rmsmax=deltrms*dint(rmsmax/deltrms)
 207       rgymin=deltrms*dint(rgymin/deltrgy)
 208       rgymax=deltrms*dint(rgymax/deltrgy)
 209       nbin_rms=(rmsmax-rmsmin)/deltrms
 210       nbin_rgy=(rgymax-rgymin)/deltrgy
 211       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
 212      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 213       nFi=0
 214       do i=1,nParmSet
 215         do j=1,nT_h(i)
 216           nFi=nFi+nR(j,i)
 217         enddo
 218       enddo
 219       write (iout,*) "nFi",nFi
 220 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 221 ! simulations.
 222 #ifdef MPI
 223       do i=1,scount(me1)
 224 #else
 225       do i=1,ntot(islice)
 226 #endif
 227 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 228         do iparm=1,nParmSet
 229 #ifdef DEBUG
 230           write (iout,'(2i5,21f8.2)') i,iparm,
 231      &     (enetb(k,i,iparm),k=1,21)
 232 #endif
 233           call restore_parm(iparm)
 234 #ifdef DEBUG
 235           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 236      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 237      &      wtor_d,wsccor,wbond
 238 #endif
 239           do ib=1,nT_h(iparm)
 240             if (rescale_mode.eq.1) then
 241               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 242               quotl=1.0d0
 243               kfacl=1.0d0
 244               do l=1,5
 245                 quotl1=quotl
 246                 quotl=quotl*quot
 247                 kfacl=kfacl*kfac
 248                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 249               enddo
 250 #if defined(FUNCTH)
 251               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 252               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 253 #elif defined(FUNCT)
 254               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 255 #else
 256               ft(6)=1.0d0
 257 #endif
 258             else if (rescale_mode.eq.2) then
 259               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 260               quotl=1.0d0
 261               do l=1,5
 262                 quotl=quotl*quot
 263                 fT(l)=1.12692801104297249644d0/
 264      &             dlog(dexp(quotl)+dexp(-quotl))
 265               enddo
 266 #if defined(FUNCTH)
 267               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 268               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 269 #elif defined(FUNCT)
 270               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 271 #else
 272               ft(6)=1.0d0
 273 #endif
 274 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 275             else if (rescale_mode.eq.0) then
 276               do l=1,6
 277                 fT(l)=1.0d0
 278               enddo
 279             else
 280               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 281      &          rescale_mode
 282               call flush(iout)
 283               return1
 284             endif
 285             evdw=enetb(1,i,iparm)
 286             evdw_t=enetb(21,i,iparm)
 287 #ifdef SCP14
 288             evdw2_14=enetb(17,i,iparm)
 289             evdw2=enetb(2,i,iparm)+evdw2_14
 290 #else
 291             evdw2=enetb(2,i,iparm)
 292             evdw2_14=0.0d0
 293 #endif
 294 #ifdef SPLITELE
 295             ees=enetb(3,i,iparm)
 296             evdw1=enetb(16,i,iparm)
 297 #else
 298             ees=enetb(3,i,iparm)
 299             evdw1=0.0d0
 300 #endif
 301             ecorr=enetb(4,i,iparm)
 302             ecorr5=enetb(5,i,iparm)
 303             ecorr6=enetb(6,i,iparm)
 304             eel_loc=enetb(7,i,iparm)
 305             eello_turn3=enetb(8,i,iparm)
 306             eello_turn4=enetb(9,i,iparm)
 307             eturn6=enetb(10,i,iparm)
 308             ebe=enetb(11,i,iparm)
 309             escloc=enetb(12,i,iparm)
 310             etors=enetb(13,i,iparm)
 311             etors_d=enetb(14,i,iparm)
 312             ehpb=enetb(15,i,iparm)
 313             estr=enetb(18,i,iparm)
 314             esccor=enetb(19,i,iparm)
 315             edihcnstr=enetb(20,i,iparm)
 316             eliptran=enetb(22,i,iparm)
 317
 318 #ifdef DEBUG
 319             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
 320      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 321      &       etors,etors_d,eello_turn3,eello_turn4,esccor
 322 #endif
 323
 324 #ifdef SPLITELE
 325             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 326      &      +wvdwpp*evdw1
 327      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 328      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 329      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 330      &      +ft(2)*wturn3*eello_turn3
 331      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 332      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 333      &      +wbond*estr+wliptran*eliptran
 334 #else
 335             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 336      &      +ft(1)*welec*(ees+evdw1)
 337      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 338      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 339      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 340      &      +ft(2)*wturn3*eello_turn3
 341      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 342      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 343      &      +wbond*estr+wliptran*eliptran
 344 #endif
 345 #ifdef DEBUG
 346             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 347      &        etot,potEmin
 348 #endif
 349 #ifdef DEBUG
 350             if (iparm.eq.1 .and. ib.eq.1) then
 351             write (iout,*)"Conformation",i
 352             energia(0)=etot
 353             do k=1,max_ene
 354               energia(k)=enetb(k,i,iparm)
 355             enddo
 356             call enerprint(energia(0),fT)
 357             endif
 358 #endif
 359             do kk=1,nR(ib,iparm)
 360               Econstr=0.0d0
 361               do j=1,nQ
 362                 dd = q(j,i)
 363                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 364      &           *(dd-q0(j,kk,ib,iparm))**2
 365               enddo
 366               v(i,kk,ib,iparm)=
 367      &          -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 368 #ifdef DEBUG
 369               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 370      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 371 #endif
 372             enddo ! kk
 373           enddo   ! ib
 374         enddo     ! iparm
 375       enddo       ! i
 376 ! Simple iteration to calculate free energies corresponding to all simulation
 377 ! runs.
 378       do iter=1,maxit
 379
 380 ! Compute new free-energy values corresponding to the righ-hand side of the
 381 ! equation and their derivatives.
 382         write (iout,*) "------------------------fi"
 383 #ifdef MPI
 384         do t=1,scount(me1)
 385 #else
 386         do t=1,ntot(islice)
 387 #endif
 388           vmax=-1.0d+20
 389           do i=1,nParmSet
 390             do k=1,nT_h(i)
 391               do l=1,nR(k,i)
 392                 vf=v(t,l,k,i)+f(l,k,i)
 393                 if (vf.gt.vmax) vmax=vf
 394               enddo
 395             enddo
 396           enddo
 397           denom=0.0d0
 398           do i=1,nParmSet
 399             do k=1,nT_h(i)
 400               do l=1,nR(k,i)
 401                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 402                 if (aux.gt.-200.0d0)
 403      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 404               enddo
 405             enddo
 406           enddo
 407           entfac(t)=-dlog(denom)-vmax
 408 #ifdef DEBUG
 409           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 410 #endif
 411         enddo
 412         do iparm=1,nParmSet
 413           do iib=1,nT_h(iparm)
 414             do ii=1,nR(iib,iparm)
 415 #ifdef MPI
 416               fi_p(ii,iib,iparm)=0.0d0
 417               do t=1,scount(me)
 418                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 419      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
 420 #ifdef DEBUG
 421               write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
 422      &         v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
 423 #endif
 424               enddo
 425 #else
 426               fi(ii,iib,iparm)=0.0d0
 427               do t=1,ntot(islice)
 428                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 429      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
 430               enddo
 431 #endif
 432             enddo ! ii
 433           enddo ! iib
 434         enddo ! iparm
 435
 436 #ifdef MPI
 437 #ifdef DEBUG
 438         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 439         do iparm=1,nParmSet
 440           do ib=1,nT_h(nparmset)
 441             write (iout,*) "iparm",iparm," ib",ib
 442             write (iout,*) "beta=",beta_h(ib,iparm)
 443             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 444           enddo
 445         enddo
 446 #endif
 447         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 448      &   maxR*MaxT_h*nParmSet
 449         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 450      &   " WHAM_COMM",WHAM_COMM
 451         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 452      &   MPI_DOUBLE_PRECISION,
 453      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 454 #ifdef DEBUG
 455         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 456         do iparm=1,nParmSet
 457           write (iout,*) "iparm",iparm
 458           do ib=1,nT_h(iparm)
 459             write (iout,*) "beta=",beta_h(ib,iparm)
 460             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 461           enddo
 462         enddo
 463 #endif
 464         if (me1.eq.Master) then
 465 #endif
 466         avefi=0.0d0
 467         do iparm=1,nParmSet
 468           do ib=1,nT_h(iparm)
 469             do i=1,nR(ib,iparm)
 470               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
 471               avefi=avefi+fi(i,ib,iparm)
 472             enddo
 473           enddo
 474         enddo
 475         avefi=avefi/nFi
 476         do iparm=1,nParmSet
 477           write (iout,*) "Parameter set",iparm
 478           do ib =1,nT_h(iparm)
 479             write (iout,*) "beta=",beta_h(ib,iparm)
 480             do i=1,nR(ib,iparm)
 481               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 482             enddo
 483             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 484             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 485           enddo
 486         enddo
 487
 488 ! Compute the norm of free-energy increments.
 489         finorm=0.0d0
 490         do iparm=1,nParmSet
 491           do ib=1,nT_h(iparm)
 492             do i=1,nR(ib,iparm)
 493               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 494               f(i,ib,iparm)=fi(i,ib,iparm)
 495             enddo
 496           enddo
 497         enddo
 498
 499         write (iout,*) 'Iteration',iter,' finorm',finorm
 500
 501 #ifdef MPI
 502         endif
 503         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 504      &   MPI_DOUBLE_PRECISION,Master,
 505      &   WHAM_COMM,IERROR)
 506         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 507      &   WHAM_COMM,IERROR)
 508 #endif
 509 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 510         if (finorm.lt.fimin) then
 511           write (iout,*) 'Iteration converged'
 512           goto 20
 513         endif
 514
 515       enddo ! iter
 516
 517    20 continue
 518 ! Now, put together the histograms from all simulations, in order to get the
 519 ! unbiased total histogram.
 520 #ifdef MPI
 521       do t=0,tmax
 522         hfin_ent_p(t)=0.0d0
 523       enddo
 524 #else
 525       do t=0,tmax
 526         hfin_ent(t)=0.0d0
 527       enddo
 528 #endif
 529       write (iout,*) "--------------hist"
 530 #ifdef MPI
 531       do iparm=1,nParmSet
 532         do i=0,nGridT
 533           sumW_p(i,iparm)=0.0d0
 534           sumE_p(i,iparm)=0.0d0
 535           sumEbis_p(i,iparm)=0.0d0
 536           sumEsq_p(i,iparm)=0.0d0
 537           do j=1,nQ+2
 538             sumQ_p(j,i,iparm)=0.0d0
 539             sumQsq_p(j,i,iparm)=0.0d0
 540             sumEQ_p(j,i,iparm)=0.0d0
 541           enddo
 542         enddo
 543       enddo
 544       upindE_p=0
 545 #else
 546       do iparm=1,nParmSet
 547         do i=0,nGridT
 548           sumW(i,iparm)=0.0d0
 549           sumE(i,iparm)=0.0d0
 550           sumEbis(i,iparm)=0.0d0
 551           sumEsq(i,iparm)=0.0d0
 552           do j=1,nQ+2
 553             sumQ(j,i,iparm)=0.0d0
 554             sumQsq(j,i,iparm)=0.0d0
 555             sumEQ(j,i,iparm)=0.0d0
 556           enddo
 557         enddo
 558       enddo
 559       upindE=0
 560 #endif
 561 c 8/26/05 entropy distribution
 562 #ifdef MPI
 563       entmin_p=1.0d10
 564       entmax_p=-1.0d10
 565       do t=1,scount(me1)
 566 c        ent=-dlog(entfac(t))
 567         ent=entfac(t)
 568         if (ent.lt.entmin_p) entmin_p=ent
 569         if (ent.gt.entmax_p) entmax_p=ent
 570       enddo
 571       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 572       call flush(iout)
 573       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 574      &  WHAM_COMM,IERROR)
 575       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 576      &  WHAM_COMM,IERROR)
 577       ientmax=entmax-entmin
 578       if (ientmax.gt.2000) ientmax=2000
 579       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 580       call flush(iout)
 581       do t=1,scount(me1)
 582 c        ient=-dlog(entfac(t))-entmin
 583         ient=entfac(t)-entmin
 584         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 585       enddo
 586       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 587      &  MPI_SUM,WHAM_COMM,IERROR)
 588       if (me1.eq.Master) then
 589         write (iout,*) "Entropy histogram"
 590         do i=0,ientmax
 591           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 592         enddo
 593       endif
 594 #else
 595       entmin=1.0d10
 596       entmax=-1.0d10
 597       do t=1,ntot(islice)
 598         ent=entfac(t)
 599         if (ent.lt.entmin) entmin=ent
 600         if (ent.gt.entmax) entmax=ent
 601       enddo
 602       ientmax=-dlog(entmax)-entmin
 603       if (ientmax.gt.2000) ientmax=2000
 604       do t=1,ntot(islice)
 605         ient=entfac(t)-entmin
 606         if (ient.le.2000) histent(ient)=histent(ient)+1
 607       enddo
 608       write (iout,*) "Entropy histogram"
 609       do i=0,ientmax
 610         write(iout,'(2f15.4)') entmin+i,histent(i)
 611       enddo
 612 #endif
 613
 614 #ifdef MPI
 615 c      write (iout,*) "me1",me1," scount",scount(me1)
 616
 617       do iparm=1,nParmSet
 618
 619 #ifdef MPI
 620         do ib=1,nT_h(iparm)
 621           do t=0,tmax
 622             hfin_p(t,ib)=0.0d0
 623           enddo
 624         enddo
 625         do i=1,maxindE
 626           histE_p(i)=0.0d0
 627         enddo
 628 #else
 629         do ib=1,nT_h(iparm)
 630           do t=0,tmax
 631             hfin(t,ib)=0.0d0
 632           enddo
 633         enddo
 634         do i=1,maxindE
 635           histE(i)=0.0d0
 636         enddo
 637 #endif
 638         do ib=1,nT_h(iparm)
 639           do i=0,MaxBinRms
 640             do j=0,MaxBinRgy
 641               hrmsrgy(j,i,ib)=0.0d0
 642 #ifdef MPI
 643               hrmsrgy_p(j,i,ib)=0.0d0
 644 #endif
 645             enddo
 646           enddo
 647         enddo
 648
 649         do t=1,scount(me1)
 650 #else
 651         do t=1,ntot(islice)
 652 #endif
 653           ind = ind_point(t)
 654 #ifdef MPI
 655           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 656 #else
 657           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 658 #endif
 659 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 660           call restore_parm(iparm)
 661           evdw=enetb(21,t,iparm)
 662           evdw_t=enetb(1,t,iparm)
 663 #ifdef SCP14
 664           evdw2_14=enetb(17,t,iparm)
 665           evdw2=enetb(2,t,iparm)+evdw2_14
 666 #else
 667           evdw2=enetb(2,t,iparm)
 668           evdw2_14=0.0d0
 669 #endif
 670 #ifdef SPLITELE
 671           ees=enetb(3,t,iparm)
 672           evdw1=enetb(16,t,iparm)
 673 #else
 674           ees=enetb(3,t,iparm)
 675           evdw1=0.0d0
 676 #endif
 677           ecorr=enetb(4,t,iparm)
 678           ecorr5=enetb(5,t,iparm)
 679           ecorr6=enetb(6,t,iparm)
 680           eel_loc=enetb(7,t,iparm)
 681           eello_turn3=enetb(8,t,iparm)
 682           eello_turn4=enetb(9,t,iparm)
 683           eturn6=enetb(10,t,iparm)
 684           ebe=enetb(11,t,iparm)
 685           escloc=enetb(12,t,iparm)
 686           etors=enetb(13,t,iparm)
 687           etors_d=enetb(14,t,iparm)
 688           ehpb=enetb(15,t,iparm)
 689           estr=enetb(18,t,iparm)
 690           esccor=enetb(19,t,iparm)
 691           edihcnstr=enetb(20,t,iparm)
 692           edihcnstr=0.0d0
 693           do k=0,nGridT
 694             betaT=startGridT+k*delta_T
 695             temper=betaT
 696 c            fT=T0/betaT
 697 c            ft=2*T0/(T0+betaT)
 698             if (rescale_mode.eq.1) then
 699               quot=betaT/T0
 700               quotl=1.0d0
 701               kfacl=1.0d0
 702               do l=1,5
 703                 quotl1=quotl
 704                 quotl=quotl*quot
 705                 kfacl=kfacl*kfac
 706                 denom=kfacl-1.0d0+quotl
 707                 fT(l)=kfacl/denom
 708                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 709                 ftbis(l)=l*kfacl*quotl1*
 710      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 711               enddo
 712 #if defined(FUNCTH)
 713               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 714      &                  320.0d0
 715               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 716              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 717      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 718 #elif defined(FUNCT)
 719               fT(6)=betaT/T0
 720               ftprim(6)=1.0d0/T0
 721               ftbis(6)=0.0d0
 722 #else
 723               fT(6)=1.0d0
 724               ftprim(6)=0.0d0
 725               ftbis(6)=0.0d0
 726 #endif
 727             else if (rescale_mode.eq.2) then
 728               quot=betaT/T0
 729               quotl=1.0d0
 730               do l=1,5
 731                 quotl1=quotl
 732                 quotl=quotl*quot
 733                 eplus=dexp(quotl)
 734                 eminus=dexp(-quotl)
 735                 logfac=1.0d0/dlog(eplus+eminus)
 736                 tanhT=(eplus-eminus)/(eplus+eminus)
 737                 fT(l)=1.12692801104297249644d0*logfac
 738                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 739                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 740      &          2*l*quotl1/T0*logfac*
 741      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 742      &          +ftprim(l)*tanhT)
 743               enddo
 744 #if defined(FUNCTH)
 745               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 746      &                 320.0d0
 747               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 748              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 749      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 750 #elif defined(FUNCT)
 751               fT(6)=betaT/T0
 752               ftprim(6)=1.0d0/T0
 753               ftbis(6)=0.0d0
 754 #else
 755               fT(6)=1.0d0
 756               ftprim(6)=0.0d0
 757               ftbis(6)=0.0d0
 758 #endif
 759             else if (rescale_mode.eq.0) then
 760               do l=1,5
 761                 fT(l)=1.0d0
 762                 ftprim(l)=0.0d0
 763               enddo
 764             else
 765               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 766      &          rescale_mode
 767               call flush(iout)
 768               return1
 769             endif
 770 c            write (iout,*) "ftprim",ftprim
 771 c            write (iout,*) "ftbis",ftbis
 772             betaT=1.0d0/(1.987D-3*betaT)
 773 #ifdef SPLITELE
 774             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 775      &      +wvdwpp*evdw1
 776      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 777      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 778      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 779      &      +ft(2)*wturn3*eello_turn3
 780      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 781      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 782      &      +wbond*estr+wliptran*eliptran
 783             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
 784      &            +ftprim(1)*wtor*etors+
 785      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
 786      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
 787      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
 788      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
 789      &            ftprim(1)*wsccor*esccor
 790             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
 791      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
 792      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
 793      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
 794      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
 795      &            ftbis(1)*wsccor*esccor
 796 #else
 797             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 798      &      +ft(1)*welec*(ees+evdw1)
 799      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 800      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 801      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 802      &      +ft(2)*wturn3*eello_turn3
 803      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 804      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 805      &      +wbond*estr+wliptran*eliptran
 806             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
 807      &           +ftprim(1)*wtor*etors+
 808      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
 809      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
 810      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
 811      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
 812      &            ftprim(1)*wsccor*esccor
 813             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
 814      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
 815      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
 816      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
 817      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
 818      &            ftprim(1)*wsccor*esccor
 819 #endif
 820             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
 821 #ifdef DEBUG
 822             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
 823      &        " etot",etot," entfac",entfac(t),
 824      &        " weight",weight," ebis",ebis
 825 #endif
 826             etot=etot-temper*eprim
 827 #ifdef MPI
 828             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
 829             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
 830             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
 831             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
 832             do j=1,nQ+2
 833               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
 834               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
 835               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
 836      &         +etot*q(j,t)*weight
 837             enddo
 838 #else
 839             sumW(k,iparm)=sumW(k,iparm)+weight
 840             sumE(k,iparm)=sumE(k,iparm)+etot*weight
 841             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
 842             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
 843             do j=1,nQ+2
 844               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
 845               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
 846               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
 847      &         +etot*q(j,t)*weight
 848             enddo
 849 #endif
 850           enddo
 851           indE = aint(potE(t,iparm)-aint(potEmin))
 852           if (indE.ge.0 .and. indE.le.maxinde) then
 853             if (indE.gt.upindE_p) upindE_p=indE
 854             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
 855           endif
 856 #ifdef MPI
 857           do ib=1,nT_h(iparm)
 858             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 859             hfin_p(ind,ib)=hfin_p(ind,ib)+
 860      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 861              if (rmsrgymap) then
 862                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 863                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 864                hrmsrgy_p(indrgy,indrms,ib)=
 865      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
 866              endif
 867           enddo
 868 #else
 869           do ib=1,nT_h(iparm)
 870             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 871             hfin(ind,ib)=hfin(ind,ib)+
 872      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 873              if (rmsrgymap) then
 874                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 875                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 876                hrmsrgy(indrgy,indrms,ib)=
 877      &           hrmsrgy(indrgy,indrms,ib)+expfac
 878              endif
 879           enddo
 880 #endif
 881         enddo ! t
 882         do ib=1,nT_h(iparm)
 883           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
 884      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 885           if (rmsrgymap) then
 886           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
 887      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
 888      &       WHAM_COMM,IERROR)
 889           endif
 890         enddo
 891         call MPI_Reduce(upindE_p,upindE,1,
 892      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
 893         call MPI_Reduce(histE_p(0),histE(0),maxindE,
 894      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 895
 896         if (me1.eq.master) then
 897
 898         if (histout) then
 899
 900         write (iout,'(6x,$)')
 901         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
 902      &   ib=1,nT_h(iparm))
 903         write (iout,*)
 904
 905         write (iout,'(/a)') 'Final histograms'
 906         if (histfile) then
 907           if (nslice.eq.1) then
 908             if (separate_parset) then
 909               write(licz3,"(bz,i3.3)") myparm
 910               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
 911             else
 912               histname=prefix(:ilen(prefix))//'.hist'
 913             endif
 914           else
 915             if (separate_parset) then
 916               write(licz3,"(bz,i3.3)") myparm
 917               histname=prefix(:ilen(prefix))//'_par'//licz3//
 918      &         '_slice_'//licz2//'.hist'
 919             else
 920               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
 921             endif
 922           endif
 923 #if defined(AIX) || defined(PGI)
 924           open (ihist,file=histname,position='append')
 925 #else
 926           open (ihist,file=histname,access='append')
 927 #endif
 928         endif
 929
 930         do t=0,tmax
 931           liczba=t
 932           sumH=0.0d0
 933           do ib=1,nT_h(iparm)
 934             sumH=sumH+hfin(t,ib)
 935           enddo
 936           if (sumH.gt.0.0d0) then
 937             do j=1,nQ
 938               jj = mod(liczba,nbin1)
 939               liczba=liczba/nbin1
 940               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 941               if (histfile)
 942      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 943             enddo
 944             do ib=1,nT_h(iparm)
 945               write (iout,'(e20.10,$)') hfin(t,ib)
 946               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
 947             enddo
 948             write (iout,'(i5)') iparm
 949             if (histfile) write (ihist,'(i5)') iparm
 950           endif
 951         enddo
 952
 953         endif
 954
 955         if (entfile) then
 956           if (nslice.eq.1) then
 957             if (separate_parset) then
 958               write(licz3,"(bz,i3.3)") myparm
 959               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
 960             else
 961               histname=prefix(:ilen(prefix))//'.ent'
 962             endif
 963           else
 964             if (separate_parset) then
 965               write(licz3,"(bz,i3.3)") myparm
 966               histname=prefix(:ilen(prefix))//'par_'//licz3//
 967      &           '_slice_'//licz2//'.ent'
 968             else
 969               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
 970             endif
 971           endif
 972 #if defined(AIX) || defined(PGI)
 973           open (ihist,file=histname,position='append')
 974 #else
 975           open (ihist,file=histname,access='append')
 976 #endif
 977           write (ihist,'(a)') "# Microcanonical entropy"
 978           do i=0,upindE
 979             write (ihist,'(f8.0,$)') dint(potEmin)+i
 980             if (histE(i).gt.0.0e0) then
 981               write (ihist,'(f15.5,$)') dlog(histE(i))
 982             else
 983               write (ihist,'(f15.5,$)') 0.0d0
 984             endif
 985           enddo
 986           write (ihist,*)
 987           close(ihist)
 988         endif
 989         write (iout,*) "Microcanonical entropy"
 990         do i=0,upindE
 991           write (iout,'(f8.0,$)') dint(potEmin)+i
 992           if (histE(i).gt.0.0e0) then
 993             write (iout,'(f15.5,$)') dlog(histE(i))
 994           else
 995             write (iout,'(f15.5,$)') 0.0d0
 996           endif
 997           write (iout,*)
 998         enddo
 999         if (rmsrgymap) then
1000           if (nslice.eq.1) then
1001             if (separate_parset) then
1002               write(licz3,"(bz,i3.3)") myparm
1003               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1004             else
1005               histname=prefix(:ilen(prefix))//'.rmsrgy'
1006             endif
1007           else
1008             if (separate_parset) then
1009               write(licz3,"(bz,i3.3)") myparm
1010               histname=prefix(:ilen(prefix))//'_par'//licz3//
1011      &         '_slice_'//licz2//'.rmsrgy'
1012             else
1013              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1014             endif
1015           endif
1016 #if defined(AIX) || defined(PGI)
1017           open (ihist,file=histname,position='append')
1018 #else
1019           open (ihist,file=histname,access='append')
1020 #endif
1021           do i=0,nbin_rms
1022             do j=0,nbin_rgy
1023               write(ihist,'(2f8.2,$)')
1024      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1025               do ib=1,nT_h(iparm)
1026                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1027                   write(ihist,'(e14.5,$)')
1028      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1029      &              +potEmin
1030                 else
1031                   write(ihist,'(e14.5,$)') 1.0d6
1032                 endif
1033               enddo
1034               write (ihist,'(i2)') iparm
1035             enddo
1036           enddo
1037           close(ihist)
1038         endif
1039         endif
1040       enddo ! iparm
1041 #ifdef MPI
1042       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1043      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1044       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1045      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1046       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1047      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1048       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1049      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1050       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1051      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1052       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1053      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1054      &   WHAM_COMM,IERROR)
1055       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1056      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1057      &   WHAM_COMM,IERROR)
1058       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1059      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1060      &   WHAM_COMM,IERROR)
1061       if (me.eq.master) then
1062 #endif
1063       write (iout,'(/a)') 'Thermal characteristics of folding'
1064       if (nslice.eq.1) then
1065         nazwa=prefix
1066       else
1067         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1068       endif
1069       iln=ilen(nazwa)
1070       if (nparmset.eq.1 .and. .not.separate_parset) then
1071         nazwa=nazwa(:iln)//".thermal"
1072       else if (nparmset.eq.1 .and. separate_parset) then
1073         write(licz3,"(bz,i3.3)") myparm
1074         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1075       endif
1076       do iparm=1,nParmSet
1077       if (nparmset.gt.1) then
1078         write(licz3,"(bz,i3.3)") iparm
1079         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1080       endif
1081       open(34,file=nazwa)
1082       if (separate_parset) then
1083         write (iout,'(a,i3)') "Parameter set",myparm
1084       else
1085         write (iout,'(a,i3)') "Parameter set",iparm
1086       endif
1087       do i=0,NGridT
1088         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1089         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1090      &    sumW(i,iparm)
1091         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1092      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1093         do j=1,nQ+2
1094           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1095           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1096      &     -sumQ(j,i,iparm)**2
1097           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1098      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1099         enddo
1100         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1101      &   (startGridT+i*delta_T))+potEmin
1102         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1103      &   sumW(i,iparm),sumE(i,iparm)
1104         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1105         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1106      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1107         write (iout,*)
1108         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1109      &   sumW(i,iparm),sumE(i,iparm)
1110         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1111         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1112      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1113         write (34,*)
1114         call flush(34)
1115       enddo
1116       close(34)
1117       enddo
1118       if (histout) then
1119       do t=0,tmax
1120         if (hfin_ent(t).gt.0.0d0) then
1121           liczba=t
1122           jj = mod(liczba,nbin1)
1123           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1124      &     hfin_ent(t)
1125           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1126      &      dmin+(jj+0.5d0)*delta,
1127      &     hfin_ent(t)
1128         endif
1129       enddo
1130       if (histfile) close(ihist)
1131       endif
1132
1133 #ifdef ZSCORE
1134 ! Write data for zscore
1135       if (nslice.eq.1) then
1136         zscname=prefix(:ilen(prefix))//".zsc"
1137       else
1138         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1139       endif
1140 #if defined(AIX) || defined(PGI)
1141       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1142 #else
1143       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1144 #endif
1145       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1146       do iparm=1,nParmSet
1147         write (izsc,'("NT=",i1)') nT_h(iparm)
1148       do ib=1,nT_h(iparm)
1149         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1150      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1151         jj = min0(nR(ib,iparm),7)
1152         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1153         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1154         write (izsc,'("&")')
1155         if (nR(ib,iparm).gt.7) then
1156           do ii=8,nR(ib,iparm),9
1157             jj = min0(nR(ib,iparm),ii+8)
1158             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1159             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1160             write (izsc,'("&")')
1161           enddo
1162         endif
1163         write (izsc,'("FI=",$)')
1164         jj=min0(nR(ib,iparm),7)
1165         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1166         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1167         write (izsc,'("&")')
1168         if (nR(ib,iparm).gt.7) then
1169           do ii=8,nR(ib,iparm),9
1170             jj = min0(nR(ib,iparm),ii+8)
1171             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1172             if (jj.eq.nR(ib,iparm)) then
1173               write (izsc,*)
1174             else
1175               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1176               write (izsc,'(t80,"&")')
1177             endif
1178           enddo
1179         endif
1180         do i=1,nR(ib,iparm)
1181           write (izsc,'("KH=",$)')
1182           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1183           write (izsc,'(" Q0=",$)')
1184           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1185           write (izsc,*)
1186         enddo
1187       enddo
1188       enddo
1189       close(izsc)
1190 #endif
1191 #ifdef MPI
1192       endif
1193 #endif
1194
1195       return
1196
1197       end