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