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