source/wham/src-M-Pawel/wham_calc1.F.safe

   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                 write (iout,*) "i,k,l",i,k,l," f",f(l,k,i),
 400      &            " v",v(t,l,k,i)," vmax",vmax,"snk",snk(l,k,i,islice),
 401      &            " aux",aux
 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 #define DEBUG
 409 #ifdef DEBUG
 410           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 411 #endif
 412         enddo
 413         do iparm=1,nParmSet
 414           do iib=1,nT_h(iparm)
 415             do ii=1,nR(iib,iparm)
 416 #ifdef MPI
 417               fi_p(ii,iib,iparm)=0.0d0
 418               do t=1,scount(me)
 419                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 420      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
 421 #ifdef DEBUG
 422               write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
 423      &         v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
 424 #endif
 425               enddo
 426 #else
 427               fi(ii,iib,iparm)=0.0d0
 428               do t=1,ntot(islice)
 429                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 430      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
 431               enddo
 432 #endif
 433             enddo ! ii
 434           enddo ! iib
 435         enddo ! iparm
 436
 437 #ifdef MPI
 438 #ifdef DEBUG
 439         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 440         do iparm=1,nParmSet
 441           do ib=1,nT_h(nparmset)
 442             write (iout,*) "iparm",iparm," ib",ib
 443             write (iout,*) "beta=",beta_h(ib,iparm)
 444             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 445           enddo
 446         enddo
 447 #endif
 448         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 449      &   maxR*MaxT_h*nParmSet
 450         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 451      &   " WHAM_COMM",WHAM_COMM
 452         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 453      &   MPI_DOUBLE_PRECISION,
 454      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 455 #ifdef DEBUG
 456         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 457         do iparm=1,nParmSet
 458           write (iout,*) "iparm",iparm
 459           do ib=1,nT_h(iparm)
 460             write (iout,*) "beta=",beta_h(ib,iparm)
 461             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 462           enddo
 463         enddo
 464 #endif
 465 #undef DEBUG
 466         if (me1.eq.Master) then
 467 #endif
 468         avefi=0.0d0
 469         do iparm=1,nParmSet
 470           do ib=1,nT_h(iparm)
 471             do i=1,nR(ib,iparm)
 472               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
 473               avefi=avefi+fi(i,ib,iparm)
 474             enddo
 475           enddo
 476         enddo
 477         avefi=avefi/nFi
 478         do iparm=1,nParmSet
 479           write (iout,*) "Parameter set",iparm
 480           do ib =1,nT_h(iparm)
 481             write (iout,*) "beta=",beta_h(ib,iparm)
 482             do i=1,nR(ib,iparm)
 483               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 484             enddo
 485             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 486             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 487           enddo
 488         enddo
 489
 490 ! Compute the norm of free-energy increments.
 491         finorm=0.0d0
 492         do iparm=1,nParmSet
 493           do ib=1,nT_h(iparm)
 494             do i=1,nR(ib,iparm)
 495               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 496               f(i,ib,iparm)=fi(i,ib,iparm)
 497             enddo
 498           enddo
 499         enddo
 500
 501         write (iout,*) 'Iteration',iter,' finorm',finorm
 502
 503 #ifdef MPI
 504         endif
 505         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 506      &   MPI_DOUBLE_PRECISION,Master,
 507      &   WHAM_COMM,IERROR)
 508         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 509      &   WHAM_COMM,IERROR)
 510 #endif
 511 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 512         if (finorm.lt.fimin) then
 513           write (iout,*) 'Iteration converged'
 514           goto 20
 515         endif
 516
 517       enddo ! iter
 518
 519    20 continue
 520 ! Now, put together the histograms from all simulations, in order to get the
 521 ! unbiased total histogram.
 522 #ifdef MPI
 523       do t=0,tmax
 524         hfin_ent_p(t)=0.0d0
 525       enddo
 526 #else
 527       do t=0,tmax
 528         hfin_ent(t)=0.0d0
 529       enddo
 530 #endif
 531       write (iout,*) "--------------hist"
 532 #ifdef MPI
 533       do iparm=1,nParmSet
 534         do i=0,nGridT
 535           sumW_p(i,iparm)=0.0d0
 536           sumE_p(i,iparm)=0.0d0
 537           sumEbis_p(i,iparm)=0.0d0
 538           sumEsq_p(i,iparm)=0.0d0
 539           do j=1,nQ+2
 540             sumQ_p(j,i,iparm)=0.0d0
 541             sumQsq_p(j,i,iparm)=0.0d0
 542             sumEQ_p(j,i,iparm)=0.0d0
 543           enddo
 544         enddo
 545       enddo
 546       upindE_p=0
 547 #else
 548       do iparm=1,nParmSet
 549         do i=0,nGridT
 550           sumW(i,iparm)=0.0d0
 551           sumE(i,iparm)=0.0d0
 552           sumEbis(i,iparm)=0.0d0
 553           sumEsq(i,iparm)=0.0d0
 554           do j=1,nQ+2
 555             sumQ(j,i,iparm)=0.0d0
 556             sumQsq(j,i,iparm)=0.0d0
 557             sumEQ(j,i,iparm)=0.0d0
 558           enddo
 559         enddo
 560       enddo
 561       upindE=0
 562 #endif
 563 c 8/26/05 entropy distribution
 564 #ifdef MPI
 565       entmin_p=1.0d10
 566       entmax_p=-1.0d10
 567       do t=1,scount(me1)
 568 c        ent=-dlog(entfac(t))
 569         ent=entfac(t)
 570         if (ent.lt.entmin_p) entmin_p=ent
 571         if (ent.gt.entmax_p) entmax_p=ent
 572       enddo
 573       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 574       call flush(iout)
 575       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 576      &  WHAM_COMM,IERROR)
 577       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 578      &  WHAM_COMM,IERROR)
 579       ientmax=entmax-entmin
 580       if (ientmax.gt.2000) ientmax=2000
 581       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 582       call flush(iout)
 583       do t=1,scount(me1)
 584 c        ient=-dlog(entfac(t))-entmin
 585         ient=entfac(t)-entmin
 586         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 587       enddo
 588       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 589      &  MPI_SUM,WHAM_COMM,IERROR)
 590       if (me1.eq.Master) then
 591         write (iout,*) "Entropy histogram"
 592         do i=0,ientmax
 593           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 594         enddo
 595       endif
 596 #else
 597       entmin=1.0d10
 598       entmax=-1.0d10
 599       do t=1,ntot(islice)
 600         ent=entfac(t)
 601         if (ent.lt.entmin) entmin=ent
 602         if (ent.gt.entmax) entmax=ent
 603       enddo
 604       ientmax=-dlog(entmax)-entmin
 605       if (ientmax.gt.2000) ientmax=2000
 606       do t=1,ntot(islice)
 607         ient=entfac(t)-entmin
 608         if (ient.le.2000) histent(ient)=histent(ient)+1
 609       enddo
 610       write (iout,*) "Entropy histogram"
 611       do i=0,ientmax
 612         write(iout,'(2f15.4)') entmin+i,histent(i)
 613       enddo
 614 #endif
 615
 616 #ifdef MPI
 617 c      write (iout,*) "me1",me1," scount",scount(me1)
 618
 619       do iparm=1,nParmSet
 620
 621 #ifdef MPI
 622         do ib=1,nT_h(iparm)
 623           do t=0,tmax
 624             hfin_p(t,ib)=0.0d0
 625           enddo
 626         enddo
 627         do i=1,maxindE
 628           histE_p(i)=0.0d0
 629         enddo
 630 #else
 631         do ib=1,nT_h(iparm)
 632           do t=0,tmax
 633             hfin(t,ib)=0.0d0
 634           enddo
 635         enddo
 636         do i=1,maxindE
 637           histE(i)=0.0d0
 638         enddo
 639 #endif
 640         do ib=1,nT_h(iparm)
 641           do i=0,MaxBinRms
 642             do j=0,MaxBinRgy
 643               hrmsrgy(j,i,ib)=0.0d0
 644 #ifdef MPI
 645               hrmsrgy_p(j,i,ib)=0.0d0
 646 #endif
 647             enddo
 648           enddo
 649         enddo
 650
 651         do t=1,scount(me1)
 652 #else
 653         do t=1,ntot(islice)
 654 #endif
 655           ind = ind_point(t)
 656 #ifdef MPI
 657           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 658 #else
 659           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 660 #endif
 661 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 662           call restore_parm(iparm)
 663           evdw=enetb(21,t,iparm)
 664           evdw_t=enetb(1,t,iparm)
 665 #ifdef SCP14
 666           evdw2_14=enetb(17,t,iparm)
 667           evdw2=enetb(2,t,iparm)+evdw2_14
 668 #else
 669           evdw2=enetb(2,t,iparm)
 670           evdw2_14=0.0d0
 671 #endif
 672 #ifdef SPLITELE
 673           ees=enetb(3,t,iparm)
 674           evdw1=enetb(16,t,iparm)
 675 #else
 676           ees=enetb(3,t,iparm)
 677           evdw1=0.0d0
 678 #endif
 679           ecorr=enetb(4,t,iparm)
 680           ecorr5=enetb(5,t,iparm)
 681           ecorr6=enetb(6,t,iparm)
 682           eel_loc=enetb(7,t,iparm)
 683           eello_turn3=enetb(8,t,iparm)
 684           eello_turn4=enetb(9,t,iparm)
 685           eturn6=enetb(10,t,iparm)
 686           ebe=enetb(11,t,iparm)
 687           escloc=enetb(12,t,iparm)
 688           etors=enetb(13,t,iparm)
 689           etors_d=enetb(14,t,iparm)
 690           ehpb=enetb(15,t,iparm)
 691           estr=enetb(18,t,iparm)
 692           esccor=enetb(19,t,iparm)
 693           edihcnstr=enetb(20,t,iparm)
 694           edihcnstr=0.0d0
 695           do k=0,nGridT
 696             betaT=startGridT+k*delta_T
 697             temper=betaT
 698 c            fT=T0/betaT
 699 c            ft=2*T0/(T0+betaT)
 700             if (rescale_mode.eq.1) then
 701               quot=betaT/T0
 702               quotl=1.0d0
 703               kfacl=1.0d0
 704               do l=1,5
 705                 quotl1=quotl
 706                 quotl=quotl*quot
 707                 kfacl=kfacl*kfac
 708                 denom=kfacl-1.0d0+quotl
 709                 fT(l)=kfacl/denom
 710                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 711                 ftbis(l)=l*kfacl*quotl1*
 712      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 713               enddo
 714 #if defined(FUNCTH)
 715               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 716      &                  320.0d0
 717               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 718              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 719      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 720 #elif defined(FUNCT)
 721               fT(6)=betaT/T0
 722               ftprim(6)=1.0d0/T0
 723               ftbis(6)=0.0d0
 724 #else
 725               fT(6)=1.0d0
 726               ftprim(6)=0.0d0
 727               ftbis(6)=0.0d0
 728 #endif
 729             else if (rescale_mode.eq.2) then
 730               quot=betaT/T0
 731               quotl=1.0d0
 732               do l=1,5
 733                 quotl1=quotl
 734                 quotl=quotl*quot
 735                 eplus=dexp(quotl)
 736                 eminus=dexp(-quotl)
 737                 logfac=1.0d0/dlog(eplus+eminus)
 738                 tanhT=(eplus-eminus)/(eplus+eminus)
 739                 fT(l)=1.12692801104297249644d0*logfac
 740                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 741                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 742      &          2*l*quotl1/T0*logfac*
 743      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 744      &          +ftprim(l)*tanhT)
 745               enddo
 746 #if defined(FUNCTH)
 747               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 748      &                 320.0d0
 749               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 750              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 751      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 752 #elif defined(FUNCT)
 753               fT(6)=betaT/T0
 754               ftprim(6)=1.0d0/T0
 755               ftbis(6)=0.0d0
 756 #else
 757               fT(6)=1.0d0
 758               ftprim(6)=0.0d0
 759               ftbis(6)=0.0d0
 760 #endif
 761             else if (rescale_mode.eq.0) then
 762               do l=1,5
 763                 fT(l)=1.0d0
 764                 ftprim(l)=0.0d0
 765               enddo
 766             else
 767               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 768      &          rescale_mode
 769               call flush(iout)
 770               return1
 771             endif
 772 c            write (iout,*) "ftprim",ftprim
 773 c            write (iout,*) "ftbis",ftbis
 774             betaT=1.0d0/(1.987D-3*betaT)
 775 #ifdef SPLITELE
 776             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 777      &      +wvdwpp*evdw1
 778      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 779      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 780      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 781      &      +ft(2)*wturn3*eello_turn3
 782      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 783      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 784      &      +wbond*estr
 785             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
 786      &            +ftprim(1)*wtor*etors+
 787      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
 788      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
 789      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
 790      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
 791      &            ftprim(1)*wsccor*esccor
 792             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
 793      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
 794      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
 795      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
 796      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
 797      &            ftbis(1)*wsccor*esccor
 798 #else
 799             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 800      &      +ft(1)*welec*(ees+evdw1)
 801      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 802      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 803      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 804      &      +ft(2)*wturn3*eello_turn3
 805      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 806      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 807      &      +wbond*estr
 808             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
 809      &           +ftprim(1)*wtor*etors+
 810      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
 811      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
 812      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
 813      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
 814      &            ftprim(1)*wsccor*esccor
 815             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
 816      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
 817      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
 818      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
 819      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
 820      &            ftprim(1)*wsccor*esccor
 821 #endif
 822             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
 823 #ifdef DEBUG
 824             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
 825      &        " etot",etot," entfac",entfac(t),
 826      &        " weight",weight," ebis",ebis
 827 #endif
 828             etot=etot-temper*eprim
 829 #ifdef MPI
 830             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
 831             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
 832             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
 833             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
 834             do j=1,nQ+2
 835               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
 836               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
 837               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
 838      &         +etot*q(j,t)*weight
 839             enddo
 840 #else
 841             sumW(k,iparm)=sumW(k,iparm)+weight
 842             sumE(k,iparm)=sumE(k,iparm)+etot*weight
 843             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
 844             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
 845             do j=1,nQ+2
 846               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
 847               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
 848               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
 849      &         +etot*q(j,t)*weight
 850             enddo
 851 #endif
 852           enddo
 853           indE = aint(potE(t,iparm)-aint(potEmin))
 854           if (indE.ge.0 .and. indE.le.maxinde) then
 855             if (indE.gt.upindE_p) upindE_p=indE
 856             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
 857           endif
 858 #ifdef MPI
 859           do ib=1,nT_h(iparm)
 860             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 861             hfin_p(ind,ib)=hfin_p(ind,ib)+
 862      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 863              if (rmsrgymap) then
 864                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 865                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 866                hrmsrgy_p(indrgy,indrms,ib)=
 867      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
 868              endif
 869           enddo
 870 #else
 871           do ib=1,nT_h(iparm)
 872             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 873             hfin(ind,ib)=hfin(ind,ib)+
 874      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
 875              if (rmsrgymap) then
 876                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
 877                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
 878                hrmsrgy(indrgy,indrms,ib)=
 879      &           hrmsrgy(indrgy,indrms,ib)+expfac
 880              endif
 881           enddo
 882 #endif
 883         enddo ! t
 884         do ib=1,nT_h(iparm)
 885           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
 886      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 887           if (rmsrgymap) then
 888           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
 889      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
 890      &       WHAM_COMM,IERROR)
 891           endif
 892         enddo
 893         call MPI_Reduce(upindE_p,upindE,1,
 894      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
 895         call MPI_Reduce(histE_p(0),histE(0),maxindE,
 896      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
 897
 898         if (me1.eq.master) then
 899
 900         if (histout) then
 901
 902         write (iout,'(6x,$)')
 903         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
 904      &   ib=1,nT_h(iparm))
 905         write (iout,*)
 906
 907         write (iout,'(/a)') 'Final histograms'
 908         if (histfile) then
 909           if (nslice.eq.1) then
 910             if (separate_parset) then
 911               write(licz3,"(bz,i3.3)") myparm
 912               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
 913             else
 914               histname=prefix(:ilen(prefix))//'.hist'
 915             endif
 916           else
 917             if (separate_parset) then
 918               write(licz3,"(bz,i3.3)") myparm
 919               histname=prefix(:ilen(prefix))//'_par'//licz3//
 920      &         '_slice_'//licz2//'.hist'
 921             else
 922               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
 923             endif
 924           endif
 925 #if defined(AIX) || defined(PGI)
 926           open (ihist,file=histname,position='append')
 927 #else
 928           open (ihist,file=histname,access='append')
 929 #endif
 930         endif
 931
 932         do t=0,tmax
 933           liczba=t
 934           sumH=0.0d0
 935           do ib=1,nT_h(iparm)
 936             sumH=sumH+hfin(t,ib)
 937           enddo
 938           if (sumH.gt.0.0d0) then
 939             do j=1,nQ
 940               jj = mod(liczba,nbin1)
 941               liczba=liczba/nbin1
 942               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 943               if (histfile)
 944      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
 945             enddo
 946             do ib=1,nT_h(iparm)
 947               write (iout,'(e20.10,$)') hfin(t,ib)
 948               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
 949             enddo
 950             write (iout,'(i5)') iparm
 951             if (histfile) write (ihist,'(i5)') iparm
 952           endif
 953         enddo
 954
 955         endif
 956
 957         if (entfile) then
 958           if (nslice.eq.1) then
 959             if (separate_parset) then
 960               write(licz3,"(bz,i3.3)") myparm
 961               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
 962             else
 963               histname=prefix(:ilen(prefix))//'.ent'
 964             endif
 965           else
 966             if (separate_parset) then
 967               write(licz3,"(bz,i3.3)") myparm
 968               histname=prefix(:ilen(prefix))//'par_'//licz3//
 969      &           '_slice_'//licz2//'.ent'
 970             else
 971               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
 972             endif
 973           endif
 974 #if defined(AIX) || defined(PGI)
 975           open (ihist,file=histname,position='append')
 976 #else
 977           open (ihist,file=histname,access='append')
 978 #endif
 979           write (ihist,'(a)') "# Microcanonical entropy"
 980           do i=0,upindE
 981             write (ihist,'(f8.0,$)') dint(potEmin)+i
 982             if (histE(i).gt.0.0e0) then
 983               write (ihist,'(f15.5,$)') dlog(histE(i))
 984             else
 985               write (ihist,'(f15.5,$)') 0.0d0
 986             endif
 987           enddo
 988           write (ihist,*)
 989           close(ihist)
 990         endif
 991         write (iout,*) "Microcanonical entropy"
 992         do i=0,upindE
 993           write (iout,'(f8.0,$)') dint(potEmin)+i
 994           if (histE(i).gt.0.0e0) then
 995             write (iout,'(f15.5,$)') dlog(histE(i))
 996           else
 997             write (iout,'(f15.5,$)') 0.0d0
 998           endif
 999           write (iout,*)
1000         enddo
1001         if (rmsrgymap) then
1002           if (nslice.eq.1) then
1003             if (separate_parset) then
1004               write(licz3,"(bz,i3.3)") myparm
1005               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1006             else
1007               histname=prefix(:ilen(prefix))//'.rmsrgy'
1008             endif
1009           else
1010             if (separate_parset) then
1011               write(licz3,"(bz,i3.3)") myparm
1012               histname=prefix(:ilen(prefix))//'_par'//licz3//
1013      &         '_slice_'//licz2//'.rmsrgy'
1014             else
1015              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1016             endif
1017           endif
1018 #if defined(AIX) || defined(PGI)
1019           open (ihist,file=histname,position='append')
1020 #else
1021           open (ihist,file=histname,access='append')
1022 #endif
1023           do i=0,nbin_rms
1024             do j=0,nbin_rgy
1025               write(ihist,'(2f8.2,$)')
1026      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1027               do ib=1,nT_h(iparm)
1028                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1029                   write(ihist,'(e14.5,$)')
1030      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1031      &              +potEmin
1032                 else
1033                   write(ihist,'(e14.5,$)') 1.0d6
1034                 endif
1035               enddo
1036               write (ihist,'(i2)') iparm
1037             enddo
1038           enddo
1039           close(ihist)
1040         endif
1041         endif
1042       enddo ! iparm
1043 #ifdef MPI
1044       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1045      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1046       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1047      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1048       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1049      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1050       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1051      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1052       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1053      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1054       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1055      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1056      &   WHAM_COMM,IERROR)
1057       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1058      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1059      &   WHAM_COMM,IERROR)
1060       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1061      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1062      &   WHAM_COMM,IERROR)
1063       if (me.eq.master) then
1064 #endif
1065       write (iout,'(/a)') 'Thermal characteristics of folding'
1066       if (nslice.eq.1) then
1067         nazwa=prefix
1068       else
1069         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1070       endif
1071       iln=ilen(nazwa)
1072       if (nparmset.eq.1 .and. .not.separate_parset) then
1073         nazwa=nazwa(:iln)//".thermal"
1074       else if (nparmset.eq.1 .and. separate_parset) then
1075         write(licz3,"(bz,i3.3)") myparm
1076         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1077       endif
1078       do iparm=1,nParmSet
1079       if (nparmset.gt.1) then
1080         write(licz3,"(bz,i3.3)") iparm
1081         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1082       endif
1083       open(34,file=nazwa)
1084       if (separate_parset) then
1085         write (iout,'(a,i3)') "Parameter set",myparm
1086       else
1087         write (iout,'(a,i3)') "Parameter set",iparm
1088       endif
1089       do i=0,NGridT
1090         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1091         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1092      &    sumW(i,iparm)
1093         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1094      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1095         do j=1,nQ+2
1096           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1097           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1098      &     -sumQ(j,i,iparm)**2
1099           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1100      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1101         enddo
1102         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1103      &   (startGridT+i*delta_T))+potEmin
1104         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1105      &   sumW(i,iparm),sumE(i,iparm)
1106         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1107         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1108      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1109         write (iout,*)
1110         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1111      &   sumW(i,iparm),sumE(i,iparm)
1112         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1113         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1114      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1115         write (34,*)
1116         call flush(34)
1117       enddo
1118       close(34)
1119       enddo
1120       if (histout) then
1121       do t=0,tmax
1122         if (hfin_ent(t).gt.0.0d0) then
1123           liczba=t
1124           jj = mod(liczba,nbin1)
1125           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1126      &     hfin_ent(t)
1127           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1128      &      dmin+(jj+0.5d0)*delta,
1129      &     hfin_ent(t)
1130         endif
1131       enddo
1132       if (histfile) close(ihist)
1133       endif
1134
1135 #ifdef ZSCORE
1136 ! Write data for zscore
1137       if (nslice.eq.1) then
1138         zscname=prefix(:ilen(prefix))//".zsc"
1139       else
1140         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1141       endif
1142 #if defined(AIX) || defined(PGI)
1143       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1144 #else
1145       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1146 #endif
1147       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1148       do iparm=1,nParmSet
1149         write (izsc,'("NT=",i1)') nT_h(iparm)
1150       do ib=1,nT_h(iparm)
1151         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1152      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1153         jj = min0(nR(ib,iparm),7)
1154         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1155         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1156         write (izsc,'("&")')
1157         if (nR(ib,iparm).gt.7) then
1158           do ii=8,nR(ib,iparm),9
1159             jj = min0(nR(ib,iparm),ii+8)
1160             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1161             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1162             write (izsc,'("&")')
1163           enddo
1164         endif
1165         write (izsc,'("FI=",$)')
1166         jj=min0(nR(ib,iparm),7)
1167         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1168         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1169         write (izsc,'("&")')
1170         if (nR(ib,iparm).gt.7) then
1171           do ii=8,nR(ib,iparm),9
1172             jj = min0(nR(ib,iparm),ii+8)
1173             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1174             if (jj.eq.nR(ib,iparm)) then
1175               write (izsc,*)
1176             else
1177               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1178               write (izsc,'(t80,"&")')
1179             endif
1180           enddo
1181         endif
1182         do i=1,nR(ib,iparm)
1183           write (izsc,'("KH=",$)')
1184           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1185           write (izsc,'(" Q0=",$)')
1186           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1187           write (izsc,*)
1188         enddo
1189       enddo
1190       enddo
1191       close(izsc)
1192 #endif
1193 #ifdef MPI
1194       endif
1195 #endif
1196
1197       return
1198
1199       end