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