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       include "COMMON.SHIELD"
  40       integer MaxPoint,MaxPointProc
  41       parameter (MaxPoint=MaxStr,
  42      & MaxPointProc=MaxStr_Proc)
  43       double precision finorm_max,potfac,entmin,entmax,expfac,vf
  44       parameter (finorm_max=1.0d0)
  45       integer islice
  46       integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  47       integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
  48      &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
  49       integer htot(0:MaxHdim),histent(0:2000)
  50       double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  51       double precision energia(0:max_ene)
  52 #ifdef MPI
  53       integer tmax_t,upindE_p
  54       double precision fi_p(MaxR,MaxT_h,Max_Parm),
  55      &  fi_p_min(MaxR,MaxT_h,Max_Parm)
  56       double precision sumW_p(0:nGridT,Max_Parm),
  57      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
  58      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
  59      & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
  60      & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
  61      & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
  62      & sumEbis_p(0:nGridT,Max_Parm)
  63       double precision hfin_p(0:MaxHdim,maxT_h),
  64      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
  65      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  66       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  67       double precision potEmin_t,entmin_p,entmax_p,
  68      & potEmin_t_all(maxT_h,Max_Parm)
  69       integer histent_p(0:2000)
  70       logical lprint /.true./
  71 #endif
  72       double precision delta_T /1.0d0/
  73       double precision rgymin,rmsmin,rgymax,rmsmax
  74       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
  75      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
  76      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
  77      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
  78      & weight,econstr
  79       double precision fi(MaxR,maxT_h,Max_Parm),
  80      & fi_min(MaxR,maxT_h,Max_Parm),
  81      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
  82      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
  83      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
  84      & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min,
  85      & hfin_ent(0:MaxHdim),vmax,aux,entfac_min
  86       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  87      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
  88      &  eplus,eminus,logfac,tanhT,tt
  89       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  90      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  91      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  92      &  eliptran,etube
  93
  94       integer ind_point(maxpoint),upindE,indE
  95       character*16 plik
  96       character*1 licz1
  97       character*2 licz2
  98       character*3 licz3
  99       character*128 nazwa
 100       integer ilen
 101       external ilen
 102
 103       write(licz2,'(bz,i2.2)') islice
 104       nbin1 = 1.0d0/delta
 105       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
 106      &  i2/80(1h-)//)') islice
 107       write (iout,*) "delta",delta," nbin1",nbin1
 108       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 109       call flush(iout)
 110       dmin=0.0d0
 111       tmax=0
 112       potEmin=1.0d10
 113        do i=1,nParmset
 114          do j=1,nT_h(i)
 115            potEmin_all(j,i)=1.0d10
 116          enddo
 117        enddo
 118       rgymin=1.0d10
 119       rmsmin=1.0d10
 120       rgymax=0.0d0
 121       rmsmax=0.0d0
 122       do t=0,MaxN
 123         htot(t)=0
 124       enddo
 125 C#define DEBUG
 126 #ifdef MPI
 127       do i=1,scount(me1)
 128 #else
 129       do i=1,ntot(islice)
 130 #endif
 131 C        do j=1,nParmSet
 132 C          if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 133 C        enddo
 134         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 135         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 136         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 137         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 138         ind_point(i)=0
 139         do j=nQ,1,-1
 140           ind=(q(j,i)-dmin+1.0d-8)/delta
 141           if (j.eq.1) then
 142             ind_point(i)=ind_point(i)+ind
 143           else
 144             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 145           endif
 146 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 147 c     &      ind_point(i)
 148           call flush(iout)
 149           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 150             write (iout,*) "Error - index exceeds range for point",i,
 151      &      " q=",q(j,i)," ind",ind_point(i)
 152 #ifdef MPI
 153             write (iout,*) "Processor",me1
 154             call flush(iout)
 155             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 156 #endif
 157             stop
 158           endif
 159         enddo ! j
 160         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 161         htot(ind_point(i))=htot(ind_point(i))+1
 162 #ifdef DEBUG
 163         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
 164      &   " htot",htot(ind_point(i))
 165         call flush(iout)
 166 #endif
 167       enddo ! i
 168       call flush(iout)
 169
 170       nbin=nbin1**nQ-1
 171       write (iout,'(a)') "Numbers of counts in Q bins"
 172       do t=0,tmax
 173         if (htot(t).gt.0) then
 174         write (iout,'(i15,$)') t
 175         liczba=t
 176         do j=1,nQ
 177           jj = mod(liczba,nbin1)
 178           liczba=liczba/nbin1
 179           write (iout,'(i5,$)') jj
 180         enddo
 181         write (iout,'(i8)') htot(t)
 182         endif
 183       enddo
 184       do iparm=1,nParmSet
 185       write (iout,'(a,i3)') "Number of data points for parameter set",
 186      & iparm
 187       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
 188      & ib=1,nT_h(iparm))
 189       write (iout,'(i8)') stot(islice)
 190       write (iout,'(a)')
 191       enddo
 192       call flush(iout)
 193
 194 #ifdef MPI
 195       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
 196      &  WHAM_COMM,IERROR)
 197       tmax=tmax_t
 198 C      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
 199 C     &  MPI_MIN,WHAM_COMM,IERROR)
 200       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
 201      &  MPI_MIN,WHAM_COMM,IERROR)
 202       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
 203      &  MPI_MAX,WHAM_COMM,IERROR)
 204       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
 205      &  MPI_MIN,WHAM_COMM,IERROR)
 206       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
 207      &  MPI_MAX,WHAM_COMM,IERROR)
 208 C      potEmin=potEmin_t/2
 209
 210       rgymin=rgymin_t
 211       rgymax=rgymax_t
 212       rmsmin=rmsmin_t
 213       rmsmax=rmsmax_t
 214       write (iout,*) "potEmin",potEmin
 215 #endif
 216       rmsmin=deltrms*dint(rmsmin/deltrms)
 217       rmsmax=deltrms*dint(rmsmax/deltrms)
 218       rgymin=deltrms*dint(rgymin/deltrgy)
 219       rgymax=deltrms*dint(rgymax/deltrgy)
 220       nbin_rms=(rmsmax-rmsmin)/deltrms
 221       nbin_rgy=(rgymax-rgymin)/deltrgy
 222       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
 223      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 224       nFi=0
 225       do i=1,nParmSet
 226         do j=1,nT_h(i)
 227           nFi=nFi+nR(j,i)
 228         enddo
 229       enddo
 230       write (iout,*) "nFi",nFi
 231 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 232 ! simulations.
 233 #ifdef MPI
 234       do i=1,scount(me1)
 235 #else
 236       do i=1,ntot(islice)
 237 #endif
 238 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 239         do iparm=1,nParmSet
 240 #ifdef DEBUG
 241           write (iout,'(2i5,21f8.2)') i,iparm,
 242      &     (enetb(k,i,iparm),k=1,max_ene)
 243 #endif
 244           call restore_parm(iparm)
 245 #ifdef DEBUG
 246           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 247      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 248      &      wtor_d,wsccor,wbond
 249 #endif
 250           do ib=1,nT_h(iparm)
 251             if (rescale_mode.eq.1) then
 252               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 253               quotl=1.0d0
 254               kfacl=1.0d0
 255               do l=1,5
 256                 quotl1=quotl
 257                 quotl=quotl*quot
 258                 kfacl=kfacl*kfac
 259                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 260               enddo
 261 #if defined(FUNCTH)
 262               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 263               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 264 #elif defined(FUNCT)
 265               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 266 #else
 267               ft(6)=1.0d0
 268 #endif
 269             else if (rescale_mode.eq.2) then
 270               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 271               quotl=1.0d0
 272               do l=1,5
 273                 quotl=quotl*quot
 274                 fT(l)=1.12692801104297249644d0/
 275      &             dlog(dexp(quotl)+dexp(-quotl))
 276               enddo
 277 #if defined(FUNCTH)
 278               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 279               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 280 #elif defined(FUNCT)
 281               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 282 #else
 283               ft(6)=1.0d0
 284 #endif
 285 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 286             else if (rescale_mode.eq.0) then
 287               do l=1,6
 288                 fT(l)=1.0d0
 289               enddo
 290             else
 291               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 292      &          rescale_mode
 293               call flush(iout)
 294               return1
 295             endif
 296             evdw=enetb(1,i,iparm)
 297             evdw_t=enetb(21,i,iparm)
 298 #ifdef SCP14
 299             evdw2_14=enetb(17,i,iparm)
 300             evdw2=enetb(2,i,iparm)+evdw2_14
 301 #else
 302             evdw2=enetb(2,i,iparm)
 303             evdw2_14=0.0d0
 304 #endif
 305 #ifdef SPLITELE
 306             ees=enetb(3,i,iparm)
 307             evdw1=enetb(16,i,iparm)
 308 #else
 309             ees=enetb(3,i,iparm)
 310             evdw1=0.0d0
 311 #endif
 312             ecorr=enetb(4,i,iparm)
 313             ecorr5=enetb(5,i,iparm)
 314             ecorr6=enetb(6,i,iparm)
 315             eel_loc=enetb(7,i,iparm)
 316             eello_turn3=enetb(8,i,iparm)
 317             eello_turn4=enetb(9,i,iparm)
 318             eturn6=enetb(10,i,iparm)
 319             ebe=enetb(11,i,iparm)
 320             escloc=enetb(12,i,iparm)
 321             etors=enetb(13,i,iparm)
 322             etors_d=enetb(14,i,iparm)
 323             ehpb=enetb(15,i,iparm)
 324             estr=enetb(18,i,iparm)
 325             esccor=enetb(19,i,iparm)
 326             edihcnstr=enetb(20,i,iparm)
 327             eliptran=enetb(22,i,iparm)
 328             etube=enetb(25,i,iparm)
 329
 330
 331 #ifdef DEBUG
 332             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
 333      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 334      &       etors,etors_d,eello_turn3,eello_turn4,esccor
 335 #endif
 336
 337 #ifdef SPLITELE
 338             if (shield_mode.gt.0) then
 339             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 340      &      +ft(1)*welec*ees
 341      &      +ft(1)*wvdwpp*evdw1
 342      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 343      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 344      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 345      &      +ft(2)*wturn3*eello_turn3
 346      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 347      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 348      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 349              else
 350             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 351      &      +wvdwpp*evdw1
 352      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 353      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 354      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 355      &      +ft(2)*wturn3*eello_turn3
 356      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 357      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 358      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 359              endif
 360 #else
 361       if (shield_mode.gt.0) then
 362             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 363      &      +ft(1)*welec*(ees+evdw1)
 364      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 365      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 366      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 367      &      +ft(2)*wturn3*eello_turn3
 368      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 369      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 370      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 371             else
 372             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 373      &      +ft(1)*welec*(ees+evdw1)
 374      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 375      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 376      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 377      &      +ft(2)*wturn3*eello_turn3
 378      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 379      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 380      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 381            endif
 382
 383 #endif
 384 #ifdef DEBUG
 385             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 386      &        etot,potEmin
 387 #endif
 388 #ifdef DEBUG
 389             if (iparm.eq.1 .and. ib.eq.1) then
 390             write (iout,*)"Conformation",i
 391             energia(0)=etot
 392             do k=1,max_ene
 393               energia(k)=enetb(k,i,iparm)
 394             enddo
 395             call enerprint(energia(0),fT)
 396             endif
 397 #endif
 398             do kk=1,nR(ib,iparm)
 399               Econstr=0.0d0
 400               do j=1,nQ
 401                 dd = q(j,i)
 402                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 403      &           *(dd-q0(j,kk,ib,iparm))**2
 404               enddo
 405               v(i,kk,ib,iparm)=
 406      &          -beta_h(ib,iparm)*(etot+Econstr)
 407 #ifdef DEBUG
 408               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 409      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 410 #endif
 411             enddo ! kk
 412           enddo   ! ib
 413         enddo     ! iparm
 414       enddo       ! i
 415 ! Simple iteration to calculate free energies corresponding to all simulation
 416 ! runs.
 417       do iter=1,maxit
 418
 419 ! Compute new free-energy values corresponding to the righ-hand side of the
 420 ! equation and their derivatives.
 421         write (iout,*) "------------------------fi"
 422         entfac_min=1.0d10
 423
 424 #ifdef MPI
 425         do t=1,scount(me1)
 426 #else
 427         do t=1,ntot(islice)
 428 #endif
 429           vmax=-1.0d+20
 430           do i=1,nParmSet
 431             do k=1,nT_h(i)
 432               do l=1,nR(k,i)
 433                 vf=v(t,l,k,i)+f(l,k,i)
 434                 if (vf.gt.vmax) vmax=vf
 435               enddo
 436             enddo
 437           enddo
 438           denom=0.0d0
 439           do i=1,nParmSet
 440             do k=1,nT_h(i)
 441               do l=1,nR(k,i)
 442                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 443                 if (aux.gt.-200.0d0)
 444      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 445               enddo
 446             enddo
 447           enddo
 448           entfac(t)=-dlog(denom)-vmax
 449           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 450 #ifdef DEBUG
 451           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 452 #endif
 453         enddo
 454
 455         do iparm=1,nParmSet
 456           do iib=1,nT_h(iparm)
 457             do ii=1,nR(iib,iparm)
 458 #ifdef MPI
 459               fi_p_min(ii,iib,iparm)=-1.0d10
 460               do t=1,scount(me)
 461                 aux=v(t,ii,iib,iparm)+entfac(t)
 462                 if (aux.gt.fi_p_min(ii,iib,iparm))
 463      &                                   fi_p_min(ii,iib,iparm)=aux
 464               enddo
 465 #else
 466               do t=1,ntot(islice)
 467                 aux=v(t,ii,iib,iparm)+entfac(t)
 468                 if (aux.gt.fi_min(ii,iib,iparm))
 469      &                                     fi_min(ii,iib,iparm)=aux
 470               enddo
 471 #endif
 472             enddo ! ii
 473           enddo ! iib
 474         enddo ! iparm
 475 #ifdef MPI
 476 #ifdef DEBUG
 477         write (iout,*) "fi_min before AllReduce"
 478         do i=1,nParmSet
 479           do j=1,nT_h(i)
 480             write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
 481           enddo
 482         enddo
 483 #endif
 484         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
 485      &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 486 #ifdef DEBUG
 487         write (iout,*) "fi_min after AllReduce"
 488         do i=1,nParmSet
 489           do j=1,nT_h(i)
 490             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 491           enddo
 492         enddo
 493 #endif
 494 #endif
 495         do iparm=1,nParmSet
 496           do iib=1,nT_h(iparm)
 497             do ii=1,nR(iib,iparm)
 498 #ifdef MPI
 499               fi_p(ii,iib,iparm)=0.0d0
 500               do t=1,scount(me)
 501                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 502      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 503 #ifdef DEBUG
 504               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
 505      &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
 506      &         fi_p(ii,iib,iparm)
 507 #endif
 508               enddo
 509 #else
 510               fi(ii,iib,iparm)=0.0d0
 511               do t=1,ntot(islice)
 512                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 513      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 514               enddo
 515 #endif
 516             enddo ! ii
 517           enddo ! iib
 518         enddo ! iparm
 519 #ifdef MPI
 520 #ifdef DEBUG
 521         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 522         do iparm=1,nParmSet
 523           do ib=1,nT_h(nparmset)
 524             write (iout,*) "iparm",iparm," ib",ib
 525             write (iout,*) "beta=",beta_h(ib,iparm)
 526             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 527           enddo
 528         enddo
 529 #endif
 530         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 531      &   maxR*MaxT_h*nParmSet
 532         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 533      &   " WHAM_COMM",WHAM_COMM
 534         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 535      &   MPI_DOUBLE_PRECISION,
 536      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 537 #ifdef DEBUG
 538         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 539         do iparm=1,nParmSet
 540           write (iout,*) "iparm",iparm
 541           do ib=1,nT_h(iparm)
 542             write (iout,*) "beta=",beta_h(ib,iparm)
 543             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 544           enddo
 545         enddo
 546 #endif
 547         if (me1.eq.Master) then
 548 #endif
 549         avefi=0.0d0
 550         do iparm=1,nParmSet
 551           do ib=1,nT_h(iparm)
 552             do i=1,nR(ib,iparm)
 553               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 554               avefi=avefi+fi(i,ib,iparm)
 555             enddo
 556           enddo
 557         enddo
 558         avefi=avefi/nFi
 559         do iparm=1,nParmSet
 560           write (iout,*) "Parameter set",iparm
 561           do ib =1,nT_h(iparm)
 562             write (iout,*) "beta=",beta_h(ib,iparm)
 563             do i=1,nR(ib,iparm)
 564               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 565             enddo
 566             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 567             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 568           enddo
 569         enddo
 570
 571 ! Compute the norm of free-energy increments.
 572         finorm=0.0d0
 573         do iparm=1,nParmSet
 574           do ib=1,nT_h(iparm)
 575             do i=1,nR(ib,iparm)
 576               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 577               f(i,ib,iparm)=fi(i,ib,iparm)
 578             enddo
 579           enddo
 580         enddo
 581
 582         write (iout,*) 'Iteration',iter,' finorm',finorm
 583
 584 #ifdef MPI
 585         endif
 586         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 587      &   MPI_DOUBLE_PRECISION,Master,
 588      &   WHAM_COMM,IERROR)
 589         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 590      &   WHAM_COMM,IERROR)
 591 #endif
 592 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 593         if (finorm.lt.fimin) then
 594           write (iout,*) 'Iteration converged'
 595           goto 20
 596         endif
 597
 598       enddo ! iter
 599
 600    20 continue
 601
 602 #ifdef MPI
 603       do i=1,scount(me1)
 604 #else
 605       do i=1,ntot(islice)
 606 #endif
 607 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 608         do iparm=1,nParmSet
 609 #ifdef DEBUG
 610           write (iout,'(2i5,21f8.2)') i,iparm,
 611      &     (enetb(k,i,iparm),k=1,21)
 612 #endif
 613           call restore_parm(iparm)
 614 #ifdef DEBUG
 615           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 616      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 617      &      wtor_d,wsccor,wbond
 618 #endif
 619           do ib=1,nT_h(iparm)
 620             if (rescale_mode.eq.1) then
 621               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 622               quotl=1.0d0
 623               kfacl=1.0d0
 624               do l=1,5
 625                 quotl1=quotl
 626                 quotl=quotl*quot
 627                 kfacl=kfacl*kfac
 628                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 629               enddo
 630 #if defined(FUNCTH)
 631               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 632               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 633 #elif defined(FUNCT)
 634               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 635 #else
 636               ft(6)=1.0d0
 637 #endif
 638             else if (rescale_mode.eq.2) then
 639               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 640               quotl=1.0d0
 641               do l=1,5
 642                 quotl=quotl*quot
 643                 fT(l)=1.12692801104297249644d0/
 644      &             dlog(dexp(quotl)+dexp(-quotl))
 645               enddo
 646 #if defined(FUNCTH)
 647               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 648               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 649 #elif defined(FUNCT)
 650               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 651 #else
 652               ft(6)=1.0d0
 653 #endif
 654 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 655             else if (rescale_mode.eq.0) then
 656               do l=1,6
 657                 fT(l)=1.0d0
 658               enddo
 659             else
 660               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 661      &          rescale_mode
 662               call flush(iout)
 663               return1
 664             endif
 665           evdw=enetb(21,i,iparm)
 666           evdw_t=enetb(1,i,iparm)
 667 #ifdef SCP14
 668           evdw2_14=enetb(17,i,iparm)
 669           evdw2=enetb(2,i,iparm)+evdw2_14
 670 #else
 671           evdw2=enetb(2,i,iparm)
 672           evdw2_14=0.0d0
 673 #endif
 674 #ifdef SPLITELE
 675           ees=enetb(3,i,iparm)
 676           evdw1=enetb(16,i,iparm)
 677 #else
 678           ees=enetb(3,i,iparm)
 679           evdw1=0.0d0
 680 #endif
 681           ecorr=enetb(4,i,iparm)
 682           ecorr5=enetb(5,i,iparm)
 683           ecorr6=enetb(6,i,iparm)
 684           eel_loc=enetb(7,i,iparm)
 685           eello_turn3=enetb(8,i,iparm)
 686           eello_turn4=enetb(9,i,iparm)
 687           eturn6=enetb(10,i,iparm)
 688           ebe=enetb(11,i,iparm)
 689           escloc=enetb(12,i,iparm)
 690           etors=enetb(13,i,iparm)
 691           etors_d=enetb(14,i,iparm)
 692           ehpb=enetb(15,i,iparm)
 693           estr=enetb(18,i,iparm)
 694           esccor=enetb(19,i,iparm)
 695           edihcnstr=enetb(20,i,iparm)
 696 C          edihcnstr=0.0d0
 697           eliptran=enetb(22,i,iparm)
 698             etube=enetb(25,i,iparm)
 699
 700 #ifdef SPLITELE
 701       if (shield_mode.gt.0) then
 702             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 703      &      +ft(1)*welec*ees
 704      &      +ft(1)*wvdwpp*evdw1
 705      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 706      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 707      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 708      &      +ft(2)*wturn3*eello_turn3
 709      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 710      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 711      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 712         else
 713             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 714      &      +wvdwpp*evdw1
 715      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 716      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 717      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 718      &      +ft(2)*wturn3*eello_turn3
 719      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 720      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 721      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 722       endif
 723 #else
 724       if (shield_mode.gt.0) then
 725             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
 726      &      +ft(1)*welec*(ees+evdw1)
 727      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 728      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 729      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 730      &      +ft(2)*wturn3*eello_turn3
 731      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 732      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 733      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 734        else
 735             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 736      &      +ft(1)*welec*(ees+evdw1)
 737      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 738      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 739      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 740      &      +ft(2)*wturn3*eello_turn3
 741      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 742      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 743      &      +wbond*estr+wliptran*eliptran+wtube*Etube
 744        endif
 745
 746 #endif
 747             etot=etot-entfac(i)/beta_h(ib,iparm)
 748             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 749
 750           enddo   ! ib
 751         enddo     ! iparm
 752       enddo       ! i
 753
 754
 755 #ifdef DEBUG
 756       write (iout,*) "The potEmin array before reduction"
 757       do i=1,nParmSet
 758         write (iout,*) "Parameter set",i
 759         do j=1,nT_h(i)
 760           write (iout,*) j,PotEmin_all(j,i)
 761         enddo
 762       enddo
 763       write (iout,*) "potEmin_min",potEmin_min
 764 #endif
 765 #ifdef MPI
 766 C Determine the minimum energes for all parameter sets and temperatures
 767       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
 768      &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 769       do i=1,nParmSet
 770         do j=1,nT_h(i)
 771           potEmin_all(j,i)=potEmin_t_all(j,i)
 772         enddo
 773       enddo
 774 #endif
 775       potEmin_min=potEmin_all(1,1)
 776       do i=1,nParmSet
 777         do j=1,nT_h(i)
 778           if (potEmin_all(j,i).lt.potEmin_min)
 779      &             potEmin_min=potEmin_all(j,i)
 780         enddo
 781       enddo
 782 #ifdef DEBUG
 783       write (iout,*) "The potEmin array"
 784       do i=1,nParmSet
 785         write (iout,*) "Parameter set",i
 786         do j=1,nT_h(i)
 787           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
 788      &      PotEmin_all(j,i)
 789         enddo
 790       enddo
 791       write (iout,*) "potEmin_min",potEmin_min
 792 #endif
 793
 794
 795 ! Now, put together the histograms from all simulations, in order to get the
 796 ! unbiased total histogram.
 797 #ifdef MPI
 798       do t=0,tmax
 799         hfin_ent_p(t)=0.0d0
 800       enddo
 801 #else
 802       do t=0,tmax
 803         hfin_ent(t)=0.0d0
 804       enddo
 805 #endif
 806       write (iout,*) "--------------hist"
 807 #ifdef MPI
 808       do iparm=1,nParmSet
 809         do i=0,nGridT
 810           sumW_p(i,iparm)=0.0d0
 811           sumE_p(i,iparm)=0.0d0
 812           sumEbis_p(i,iparm)=0.0d0
 813           sumEsq_p(i,iparm)=0.0d0
 814           do j=1,nQ+2
 815             sumQ_p(j,i,iparm)=0.0d0
 816             sumQsq_p(j,i,iparm)=0.0d0
 817             sumEQ_p(j,i,iparm)=0.0d0
 818           enddo
 819         enddo
 820       enddo
 821       upindE_p=0
 822 #else
 823       do iparm=1,nParmSet
 824         do i=0,nGridT
 825           sumW(i,iparm)=0.0d0
 826           sumE(i,iparm)=0.0d0
 827           sumEbis(i,iparm)=0.0d0
 828           sumEsq(i,iparm)=0.0d0
 829           do j=1,nQ+2
 830             sumQ(j,i,iparm)=0.0d0
 831             sumQsq(j,i,iparm)=0.0d0
 832             sumEQ(j,i,iparm)=0.0d0
 833           enddo
 834         enddo
 835       enddo
 836       upindE=0
 837 #endif
 838 c 8/26/05 entropy distribution
 839 #ifdef MPI
 840       entmin_p=1.0d10
 841       entmax_p=-1.0d10
 842       do t=1,scount(me1)
 843 c        ent=-dlog(entfac(t))
 844         ent=entfac(t)
 845         if (ent.lt.entmin_p) entmin_p=ent
 846         if (ent.gt.entmax_p) entmax_p=ent
 847       enddo
 848       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 849       call flush(iout)
 850       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 851      &  WHAM_COMM,IERROR)
 852       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 853      &  WHAM_COMM,IERROR)
 854 C      ientmax=entmax-entmin
 855 C      if (ientmax.gt.2000) ientmax=2000
 856       if ((-dlog(entmax)-entmin).lt.2000.0d0) then
 857       ientmax=-dlog(entmax)-entmin
 858       else
 859        ientmax=2000
 860       endif
 861       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 862       call flush(iout)
 863       do t=1,scount(me1)
 864 c        ient=-dlog(entfac(t))-entmin
 865         ient=entfac(t)-entmin
 866       write (iout,*) "ient",ient,entfac(t),entmin
 867 C        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 868       enddo
 869 C      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 870 C     &  MPI_SUM,WHAM_COMM,IERROR)
 871 C      if (me1.eq.Master) then
 872 C        write (iout,*) "Entropy histogram"
 873 C        do i=0,ientmax
 874 C          write(iout,'(f15.4,i10)') entmin+i,histent(i)
 875 C        enddo
 876 C      endif
 877 #else
 878       entmin=1.0d10
 879       entmax=-1.0d10
 880       do t=1,ntot(islice)
 881         ent=entfac(t)
 882         if (ent.lt.entmin) entmin=ent
 883         if (ent.gt.entmax) entmax=ent
 884       enddo
 885       if ((-dlog(entmax)-entmin).lt.2000.0d0) then
 886       ientmax=-dlog(entmax)-entmin
 887       else
 888        ientmax=2000
 889       endif
 890 C      do t=1,ntot(islice)
 891 C        ient=entfac(t)-entmin
 892 C        if (ient.le.2000) histent(ient)=histent(ient)+1
 893 C      enddo
 894 C      write (iout,*) "Entropy histogram"
 895 C      do i=0,ientmax
 896 C        write(iout,'(2f15.4)') entmin+i,histent(i)
 897 C      enddo
 898 #endif
 899
 900 #ifdef MPI
 901 c      write (iout,*) "me1",me1," scount",scount(me1)
 902
 903       do iparm=1,nParmSet
 904
 905 #ifdef MPI
 906         do ib=1,nT_h(iparm)
 907           do t=0,tmax
 908             hfin_p(t,ib)=0.0d0
 909           enddo
 910         enddo
 911         do i=1,maxindE
 912           histE_p(i)=0.0d0
 913         enddo
 914 #else
 915         do ib=1,nT_h(iparm)
 916           do t=0,tmax
 917             hfin(t,ib)=0.0d0
 918           enddo
 919         enddo
 920         do i=1,maxindE
 921           histE(i)=0.0d0
 922         enddo
 923 #endif
 924         do ib=1,nT_h(iparm)
 925           do i=0,MaxBinRms
 926             do j=0,MaxBinRgy
 927               hrmsrgy(j,i,ib)=0.0d0
 928 #ifdef MPI
 929               hrmsrgy_p(j,i,ib)=0.0d0
 930 #endif
 931             enddo
 932           enddo
 933         enddo
 934
 935         do t=1,scount(me1)
 936 #else
 937         do t=1,ntot(islice)
 938 #endif
 939           ind = ind_point(t)
 940 #ifdef MPI
 941           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 942 #else
 943           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 944 #endif
 945 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 946           call restore_parm(iparm)
 947           evdw=enetb(21,t,iparm)
 948           evdw_t=enetb(1,t,iparm)
 949 #ifdef SCP14
 950           evdw2_14=enetb(17,t,iparm)
 951           evdw2=enetb(2,t,iparm)+evdw2_14
 952 #else
 953           evdw2=enetb(2,t,iparm)
 954           evdw2_14=0.0d0
 955 #endif
 956 #ifdef SPLITELE
 957           ees=enetb(3,t,iparm)
 958           evdw1=enetb(16,t,iparm)
 959 #else
 960           ees=enetb(3,t,iparm)
 961           evdw1=0.0d0
 962 #endif
 963           ecorr=enetb(4,t,iparm)
 964           ecorr5=enetb(5,t,iparm)
 965           ecorr6=enetb(6,t,iparm)
 966           eel_loc=enetb(7,t,iparm)
 967           eello_turn3=enetb(8,t,iparm)
 968           eello_turn4=enetb(9,t,iparm)
 969           eturn6=enetb(10,t,iparm)
 970           ebe=enetb(11,t,iparm)
 971           escloc=enetb(12,t,iparm)
 972           etors=enetb(13,t,iparm)
 973           etors_d=enetb(14,t,iparm)
 974           ehpb=enetb(15,t,iparm)
 975           estr=enetb(18,t,iparm)
 976           esccor=enetb(19,t,iparm)
 977           edihcnstr=enetb(20,t,iparm)
 978 C          edihcnstr=0.0d0
 979           eliptran=enetb(22,i,iparm)
 980             etube=enetb(25,i,iparm)
 981
 982           do k=0,nGridT
 983             betaT=startGridT+k*delta_T
 984             temper=betaT
 985 c            fT=T0/betaT
 986 c            ft=2*T0/(T0+betaT)
 987             if (rescale_mode.eq.1) then
 988               quot=betaT/T0
 989               quotl=1.0d0
 990               kfacl=1.0d0
 991               do l=1,5
 992                 quotl1=quotl
 993                 quotl=quotl*quot
 994                 kfacl=kfacl*kfac
 995                 denom=kfacl-1.0d0+quotl
 996                 fT(l)=kfacl/denom
 997                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 998                 ftbis(l)=l*kfacl*quotl1*
 999      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1000               enddo
1001 #if defined(FUNCTH)
1002               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1003      &                  320.0d0
1004               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1005              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1006      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1007 #elif defined(FUNCT)
1008               fT(6)=betaT/T0
1009               ftprim(6)=1.0d0/T0
1010               ftbis(6)=0.0d0
1011 #else
1012               fT(6)=1.0d0
1013               ftprim(6)=0.0d0
1014               ftbis(6)=0.0d0
1015 #endif
1016             else if (rescale_mode.eq.2) then
1017               quot=betaT/T0
1018               quotl=1.0d0
1019               do l=1,5
1020                 quotl1=quotl
1021                 quotl=quotl*quot
1022                 eplus=dexp(quotl)
1023                 eminus=dexp(-quotl)
1024                 logfac=1.0d0/dlog(eplus+eminus)
1025                 tanhT=(eplus-eminus)/(eplus+eminus)
1026                 fT(l)=1.12692801104297249644d0*logfac
1027                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1028                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1029      &          2*l*quotl1/T0*logfac*
1030      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1031      &          +ftprim(l)*tanhT)
1032               enddo
1033 #if defined(FUNCTH)
1034               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1035      &                 320.0d0
1036               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1037              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1038      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1039 #elif defined(FUNCT)
1040               fT(6)=betaT/T0
1041               ftprim(6)=1.0d0/T0
1042               ftbis(6)=0.0d0
1043 #else
1044               fT(6)=1.0d0
1045               ftprim(6)=0.0d0
1046               ftbis(6)=0.0d0
1047 #endif
1048             else if (rescale_mode.eq.0) then
1049               do l=1,5
1050                 fT(l)=1.0d0
1051                 ftprim(l)=0.0d0
1052               enddo
1053             else
1054               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1055      &          rescale_mode
1056               call flush(iout)
1057               return1
1058             endif
1059 c            write (iout,*) "ftprim",ftprim
1060 c            write (iout,*) "ftbis",ftbis
1061             betaT=1.0d0/(1.987D-3*betaT)
1062             if (betaT.ge.beta_h(1,iparm)) then
1063               potEmin=potEmin_all(1,iparm)+
1064      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1065      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1066      &         (1.0/betaT-1.0/beta_h(1,iparm))
1067 #ifdef DEBUG
1068               write(iout,*) "first",temper,potEmin
1069 #endif
1070             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1071               potEmin=potEmin_all(nT_h(iparm),iparm)+
1072      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1073      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1074      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1075 #ifdef DEBUG
1076               write (iout,*) "last",temper,potEmin
1077 #endif
1078             else
1079               do l=1,nT_h(iparm)-1
1080                 if (betaT.le.beta_h(l,iparm) .and.
1081      &              betaT.gt.beta_h(l+1,iparm)) then
1082                   potEmin=potEmin_all(l,iparm)
1083 #ifdef DEBUG
1084                   write (iout,*) "l",l,
1085      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1086      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1087 #endif
1088                   exit
1089                 endif
1090               enddo
1091             endif
1092 #ifdef SPLITELE
1093       if (shield_mode.gt.0) then
1094             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1095      &      +ft(1)*welec*ees
1096      &      +ft(1)*wvdwpp*evdw1
1097      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1098      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1099      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1100      &      +ft(2)*wturn3*eello_turn3
1101      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1102      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1103      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1104             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1105 C     &            +ftprim(6)*evdw_t
1106      &            +ftprim(1)*wscp*evdw2
1107      &            +ftprim(1)*welec*ees
1108      &            +ftprim(1)*wvdwpp*evdw1
1109      &            +ftprim(1)*wtor*etors+
1110      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1111      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1112      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1113      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1114      &            ftprim(1)*wsccor*esccor
1115             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1116      &            +ftbis(1)*wscp*evdw2+
1117      &            ftbis(1)*welec*ees
1118      &            +ftbis(1)*wvdwpp*evdw
1119      &            +ftbis(1)*wtor*etors+
1120      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1121      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1122      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1123      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1124      &            ftbis(1)*wsccor*esccor
1125         else
1126             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1127      &      +wvdwpp*evdw1
1128      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1129      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1130      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1131      &      +ft(2)*wturn3*eello_turn3
1132      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1133      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1134      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1135             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1136      &            +ftprim(1)*wtor*etors+
1137      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1138      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1139      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1140      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1141      &            ftprim(1)*wsccor*esccor
1142             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1143      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1144      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1145      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1146      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1147      &            ftbis(1)*wsccor*esccor
1148       endif
1149 #else
1150       if (shield_mode.gt.0) then
1151             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1152      &      +ft(1)*welec*(ees+evdw1)
1153      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1154      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1155      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1156      &      +ft(2)*wturn3*eello_turn3
1157      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1158      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1159      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1160             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1161      &           +ftprim(1)*welec*(ees+evdw1)
1162      &           +ftprim(1)*wtor*etors+
1163      &            ftprim(1)*wscp*evdw2+
1164      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1165      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1166      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1167      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1168      &            ftprim(1)*wsccor*esccor
1169             ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1170      &           +ftbis(1)*wscp*evdw2
1171      &           +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1172      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1173      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1174      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1175      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1176      &            ftbis(1)*wsccor*esccor
1177        else
1178             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1179      &      +ft(1)*welec*(ees+evdw1)
1180      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1181      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1182      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1183      &      +ft(2)*wturn3*eello_turn3
1184      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1185      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1186      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1187             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1188      &           +ftprim(1)*wtor*etors+
1189      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1190      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1191      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1192      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1193      &            ftprim(1)*wsccor*esccor
1194             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1195      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1196      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1197      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1198      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1199      &            ftbis(1)*wsccor*esccor
1200
1201        endif
1202
1203 #endif
1204             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1205 #ifdef DEBUG
1206             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1207      &        " etot",etot," entfac",entfac(t),
1208      &        " weight",weight," ebis",ebis
1209 #endif
1210             etot=etot-temper*eprim
1211 #ifdef MPI
1212             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1213             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1214             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1215             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1216             do j=1,nQ+2
1217               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1218               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1219               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1220      &         +etot*q(j,t)*weight
1221             enddo
1222 #else
1223             sumW(k,iparm)=sumW(k,iparm)+weight
1224             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1225             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1226             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1227             do j=1,nQ+2
1228               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1229               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1230               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1231      &         +etot*q(j,t)*weight
1232             enddo
1233 #endif
1234           enddo
1235           indE = aint(potE(t,iparm)-aint(potEmin))
1236           if (indE.ge.0 .and. indE.le.maxinde) then
1237             if (indE.gt.upindE_p) upindE_p=indE
1238             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1239           endif
1240 #ifdef MPI
1241           do ib=1,nT_h(iparm)
1242             potEmin=potEmin_all(ib,iparm)
1243             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1244             hfin_p(ind,ib)=hfin_p(ind,ib)+
1245      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1246              if (rmsrgymap) then
1247                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1248                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1249                hrmsrgy_p(indrgy,indrms,ib)=
1250      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1251              endif
1252           enddo
1253 #else
1254           do ib=1,nT_h(iparm)
1255             potEmin=potEmin_all(ib,iparm)
1256             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1257             hfin(ind,ib)=hfin(ind,ib)+
1258      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1259              if (rmsrgymap) then
1260                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1261                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1262                hrmsrgy(indrgy,indrms,ib)=
1263      &           hrmsrgy(indrgy,indrms,ib)+expfac
1264              endif
1265           enddo
1266 #endif
1267         enddo ! t
1268         do ib=1,nT_h(iparm)
1269           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1270      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1271           if (rmsrgymap) then
1272           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1273      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1274      &       WHAM_COMM,IERROR)
1275           endif
1276         enddo
1277         call MPI_Reduce(upindE_p,upindE,1,
1278      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1279         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1280      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1281
1282         if (me1.eq.master) then
1283
1284         if (histout) then
1285
1286         write (iout,'(6x,$)')
1287         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1288      &   ib=1,nT_h(iparm))
1289         write (iout,*)
1290
1291         write (iout,'(/a)') 'Final histograms'
1292         if (histfile) then
1293           if (nslice.eq.1) then
1294             if (separate_parset) then
1295               write(licz3,"(bz,i3.3)") myparm
1296               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1297             else
1298               histname=prefix(:ilen(prefix))//'.hist'
1299             endif
1300           else
1301             if (separate_parset) then
1302               write(licz3,"(bz,i3.3)") myparm
1303               histname=prefix(:ilen(prefix))//'_par'//licz3//
1304      &         '_slice_'//licz2//'.hist'
1305             else
1306               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1307             endif
1308           endif
1309 #if defined(AIX) || defined(PGI)
1310           open (ihist,file=histname,position='append')
1311 #else
1312           open (ihist,file=histname,access='append')
1313 #endif
1314         endif
1315
1316         do t=0,tmax
1317           liczba=t
1318           sumH=0.0d0
1319           do ib=1,nT_h(iparm)
1320             sumH=sumH+hfin(t,ib)
1321           enddo
1322           if (sumH.gt.0.0d0) then
1323             do j=1,nQ
1324               jj = mod(liczba,nbin1)
1325               liczba=liczba/nbin1
1326               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1327               if (histfile)
1328      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1329             enddo
1330             do ib=1,nT_h(iparm)
1331               write (iout,'(e20.10,$)') hfin(t,ib)
1332               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1333             enddo
1334             write (iout,'(i5)') iparm
1335             if (histfile) write (ihist,'(i5)') iparm
1336           endif
1337         enddo
1338
1339         endif
1340
1341         if (entfile) then
1342           if (nslice.eq.1) then
1343             if (separate_parset) then
1344               write(licz3,"(bz,i3.3)") myparm
1345               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1346             else
1347               histname=prefix(:ilen(prefix))//'.ent'
1348             endif
1349           else
1350             if (separate_parset) then
1351               write(licz3,"(bz,i3.3)") myparm
1352               histname=prefix(:ilen(prefix))//'par_'//licz3//
1353      &           '_slice_'//licz2//'.ent'
1354             else
1355               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1356             endif
1357           endif
1358 #if defined(AIX) || defined(PGI)
1359           open (ihist,file=histname,position='append')
1360 #else
1361           open (ihist,file=histname,access='append')
1362 #endif
1363           write (ihist,'(a)') "# Microcanonical entropy"
1364           do i=0,upindE
1365             write (ihist,'(f8.0,$)') dint(potEmin)+i
1366             if (histE(i).gt.0.0e0) then
1367               write (ihist,'(f15.5,$)') dlog(histE(i))
1368             else
1369               write (ihist,'(f15.5,$)') 0.0d0
1370             endif
1371           enddo
1372           write (ihist,*)
1373           close(ihist)
1374         endif
1375         write (iout,*) "Microcanonical entropy"
1376         do i=0,upindE
1377           write (iout,'(f8.0,$)') dint(potEmin)+i
1378           if (histE(i).gt.0.0e0) then
1379             write (iout,'(f15.5,$)') dlog(histE(i))
1380           else
1381             write (iout,'(f15.5,$)') 0.0d0
1382           endif
1383           write (iout,*)
1384         enddo
1385         if (rmsrgymap) then
1386           if (nslice.eq.1) then
1387             if (separate_parset) then
1388               write(licz3,"(bz,i3.3)") myparm
1389               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1390             else
1391               histname=prefix(:ilen(prefix))//'.rmsrgy'
1392             endif
1393           else
1394             if (separate_parset) then
1395               write(licz3,"(bz,i3.3)") myparm
1396               histname=prefix(:ilen(prefix))//'_par'//licz3//
1397      &         '_slice_'//licz2//'.rmsrgy'
1398             else
1399              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1400             endif
1401           endif
1402 #if defined(AIX) || defined(PGI)
1403           open (ihist,file=histname,position='append')
1404 #else
1405           open (ihist,file=histname,access='append')
1406 #endif
1407           do i=0,nbin_rms
1408             do j=0,nbin_rgy
1409               write(ihist,'(2f8.2,$)')
1410      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1411               do ib=1,nT_h(iparm)
1412                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1413                   write(ihist,'(e14.5,$)')
1414      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1415      &              +potEmin
1416                 else
1417                   write(ihist,'(e14.5,$)') 1.0d6
1418                 endif
1419               enddo
1420               write (ihist,'(i2)') iparm
1421             enddo
1422           enddo
1423           close(ihist)
1424         endif
1425         endif
1426       enddo ! iparm
1427 #ifdef MPI
1428       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1429      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1430       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1431      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1432       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1433      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1434       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1435      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1436       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1437      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1438       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1439      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1440      &   WHAM_COMM,IERROR)
1441       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1442      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1443      &   WHAM_COMM,IERROR)
1444       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1445      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1446      &   WHAM_COMM,IERROR)
1447       if (me.eq.master) then
1448 #endif
1449       write (iout,'(/a)') 'Thermal characteristics of folding'
1450       if (nslice.eq.1) then
1451         nazwa=prefix
1452       else
1453         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1454       endif
1455       iln=ilen(nazwa)
1456       if (nparmset.eq.1 .and. .not.separate_parset) then
1457         nazwa=nazwa(:iln)//".thermal"
1458       else if (nparmset.eq.1 .and. separate_parset) then
1459         write(licz3,"(bz,i3.3)") myparm
1460         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1461       endif
1462       do iparm=1,nParmSet
1463       if (nparmset.gt.1) then
1464         write(licz3,"(bz,i3.3)") iparm
1465         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1466       endif
1467       open(34,file=nazwa)
1468       if (separate_parset) then
1469         write (iout,'(a,i3)') "Parameter set",myparm
1470       else
1471         write (iout,'(a,i3)') "Parameter set",iparm
1472       endif
1473       do i=0,NGridT
1474         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1475         if (betaT.ge.beta_h(1,iparm)) then
1476           potEmin=potEmin_all(1,iparm)
1477         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1478           potEmin=potEmin_all(nT_h(iparm),iparm)
1479         else
1480           do l=1,nT_h(iparm)-1
1481             if (betaT.le.beta_h(l,iparm) .and.
1482      &          betaT.gt.beta_h(l+1,iparm)) then
1483               potEmin=potEmin_all(l,iparm)
1484               exit
1485             endif
1486           enddo
1487         endif
1488
1489         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1490         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1491      &    sumW(i,iparm)
1492         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1493      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1494         do j=1,nQ+2
1495           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1496           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1497      &     -sumQ(j,i,iparm)**2
1498           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1499      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1500         enddo
1501         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1502      &   (startGridT+i*delta_T))+potEmin
1503         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1504      &   sumW(i,iparm),sumE(i,iparm)
1505         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1506         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1507      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1508         write (iout,*)
1509         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1510      &   sumW(i,iparm),sumE(i,iparm)
1511         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1512         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1513      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1514         write (34,*)
1515         call flush(34)
1516       enddo
1517       close(34)
1518       enddo
1519       if (histout) then
1520       do t=0,tmax
1521         if (hfin_ent(t).gt.0.0d0) then
1522           liczba=t
1523           jj = mod(liczba,nbin1)
1524           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1525      &     hfin_ent(t)
1526           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1527      &      dmin+(jj+0.5d0)*delta,
1528      &     hfin_ent(t)
1529         endif
1530       enddo
1531       if (histfile) close(ihist)
1532       endif
1533
1534 #ifdef ZSCORE
1535 ! Write data for zscore
1536       if (nslice.eq.1) then
1537         zscname=prefix(:ilen(prefix))//".zsc"
1538       else
1539         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1540       endif
1541 #if defined(AIX) || defined(PGI)
1542       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1543 #else
1544       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1545 #endif
1546       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1547       do iparm=1,nParmSet
1548         write (izsc,'("NT=",i1)') nT_h(iparm)
1549       do ib=1,nT_h(iparm)
1550         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1551      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1552         jj = min0(nR(ib,iparm),7)
1553         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1554         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1555         write (izsc,'("&")')
1556         if (nR(ib,iparm).gt.7) then
1557           do ii=8,nR(ib,iparm),9
1558             jj = min0(nR(ib,iparm),ii+8)
1559             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1560             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1561             write (izsc,'("&")')
1562           enddo
1563         endif
1564         write (izsc,'("FI=",$)')
1565         jj=min0(nR(ib,iparm),7)
1566         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1567         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1568         write (izsc,'("&")')
1569         if (nR(ib,iparm).gt.7) then
1570           do ii=8,nR(ib,iparm),9
1571             jj = min0(nR(ib,iparm),ii+8)
1572             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1573             if (jj.eq.nR(ib,iparm)) then
1574               write (izsc,*)
1575             else
1576               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1577               write (izsc,'(t80,"&")')
1578             endif
1579           enddo
1580         endif
1581         do i=1,nR(ib,iparm)
1582           write (izsc,'("KH=",$)')
1583           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1584           write (izsc,'(" Q0=",$)')
1585           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1586           write (izsc,*)
1587         enddo
1588       enddo
1589       enddo
1590       close(izsc)
1591 #endif
1592 #ifdef MPI
1593       endif
1594 #endif
1595
1596       return
1597 C#undef DEBUG
1598       end