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         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 867       enddo
 868       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 869      &  MPI_SUM,WHAM_COMM,IERROR)
 870       if (me1.eq.Master) then
 871         write (iout,*) "Entropy histogram"
 872         do i=0,ientmax
 873           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 874         enddo
 875       endif
 876 #else
 877       entmin=1.0d10
 878       entmax=-1.0d10
 879       do t=1,ntot(islice)
 880         ent=entfac(t)
 881         if (ent.lt.entmin) entmin=ent
 882         if (ent.gt.entmax) entmax=ent
 883       enddo
 884       if ((-dlog(entmax)-entmin).lt.2000.0d0) then
 885       ientmax=-dlog(entmax)-entmin
 886       else
 887        ientmax=2000
 888       endif
 889       do t=1,ntot(islice)
 890         ient=entfac(t)-entmin
 891         if (ient.le.2000) histent(ient)=histent(ient)+1
 892       enddo
 893       write (iout,*) "Entropy histogram"
 894       do i=0,ientmax
 895         write(iout,'(2f15.4)') entmin+i,histent(i)
 896       enddo
 897 #endif
 898
 899 #ifdef MPI
 900 c      write (iout,*) "me1",me1," scount",scount(me1)
 901
 902       do iparm=1,nParmSet
 903
 904 #ifdef MPI
 905         do ib=1,nT_h(iparm)
 906           do t=0,tmax
 907             hfin_p(t,ib)=0.0d0
 908           enddo
 909         enddo
 910         do i=1,maxindE
 911           histE_p(i)=0.0d0
 912         enddo
 913 #else
 914         do ib=1,nT_h(iparm)
 915           do t=0,tmax
 916             hfin(t,ib)=0.0d0
 917           enddo
 918         enddo
 919         do i=1,maxindE
 920           histE(i)=0.0d0
 921         enddo
 922 #endif
 923         do ib=1,nT_h(iparm)
 924           do i=0,MaxBinRms
 925             do j=0,MaxBinRgy
 926               hrmsrgy(j,i,ib)=0.0d0
 927 #ifdef MPI
 928               hrmsrgy_p(j,i,ib)=0.0d0
 929 #endif
 930             enddo
 931           enddo
 932         enddo
 933
 934         do t=1,scount(me1)
 935 #else
 936         do t=1,ntot(islice)
 937 #endif
 938           ind = ind_point(t)
 939 #ifdef MPI
 940           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 941 #else
 942           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 943 #endif
 944 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
 945           call restore_parm(iparm)
 946           evdw=enetb(21,t,iparm)
 947           evdw_t=enetb(1,t,iparm)
 948 #ifdef SCP14
 949           evdw2_14=enetb(17,t,iparm)
 950           evdw2=enetb(2,t,iparm)+evdw2_14
 951 #else
 952           evdw2=enetb(2,t,iparm)
 953           evdw2_14=0.0d0
 954 #endif
 955 #ifdef SPLITELE
 956           ees=enetb(3,t,iparm)
 957           evdw1=enetb(16,t,iparm)
 958 #else
 959           ees=enetb(3,t,iparm)
 960           evdw1=0.0d0
 961 #endif
 962           ecorr=enetb(4,t,iparm)
 963           ecorr5=enetb(5,t,iparm)
 964           ecorr6=enetb(6,t,iparm)
 965           eel_loc=enetb(7,t,iparm)
 966           eello_turn3=enetb(8,t,iparm)
 967           eello_turn4=enetb(9,t,iparm)
 968           eturn6=enetb(10,t,iparm)
 969           ebe=enetb(11,t,iparm)
 970           escloc=enetb(12,t,iparm)
 971           etors=enetb(13,t,iparm)
 972           etors_d=enetb(14,t,iparm)
 973           ehpb=enetb(15,t,iparm)
 974           estr=enetb(18,t,iparm)
 975           esccor=enetb(19,t,iparm)
 976           edihcnstr=enetb(20,t,iparm)
 977 C          edihcnstr=0.0d0
 978           eliptran=enetb(22,i,iparm)
 979             etube=enetb(25,i,iparm)
 980
 981           do k=0,nGridT
 982             betaT=startGridT+k*delta_T
 983             temper=betaT
 984 c            fT=T0/betaT
 985 c            ft=2*T0/(T0+betaT)
 986             if (rescale_mode.eq.1) then
 987               quot=betaT/T0
 988               quotl=1.0d0
 989               kfacl=1.0d0
 990               do l=1,5
 991                 quotl1=quotl
 992                 quotl=quotl*quot
 993                 kfacl=kfacl*kfac
 994                 denom=kfacl-1.0d0+quotl
 995                 fT(l)=kfacl/denom
 996                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 997                 ftbis(l)=l*kfacl*quotl1*
 998      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 999               enddo
1000 #if defined(FUNCTH)
1001               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1002      &                  320.0d0
1003               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1004              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1005      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1006 #elif defined(FUNCT)
1007               fT(6)=betaT/T0
1008               ftprim(6)=1.0d0/T0
1009               ftbis(6)=0.0d0
1010 #else
1011               fT(6)=1.0d0
1012               ftprim(6)=0.0d0
1013               ftbis(6)=0.0d0
1014 #endif
1015             else if (rescale_mode.eq.2) then
1016               quot=betaT/T0
1017               quotl=1.0d0
1018               do l=1,5
1019                 quotl1=quotl
1020                 quotl=quotl*quot
1021                 eplus=dexp(quotl)
1022                 eminus=dexp(-quotl)
1023                 logfac=1.0d0/dlog(eplus+eminus)
1024                 tanhT=(eplus-eminus)/(eplus+eminus)
1025                 fT(l)=1.12692801104297249644d0*logfac
1026                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1027                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1028      &          2*l*quotl1/T0*logfac*
1029      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1030      &          +ftprim(l)*tanhT)
1031               enddo
1032 #if defined(FUNCTH)
1033               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1034      &                 320.0d0
1035               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1036              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1037      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1038 #elif defined(FUNCT)
1039               fT(6)=betaT/T0
1040               ftprim(6)=1.0d0/T0
1041               ftbis(6)=0.0d0
1042 #else
1043               fT(6)=1.0d0
1044               ftprim(6)=0.0d0
1045               ftbis(6)=0.0d0
1046 #endif
1047             else if (rescale_mode.eq.0) then
1048               do l=1,5
1049                 fT(l)=1.0d0
1050                 ftprim(l)=0.0d0
1051               enddo
1052             else
1053               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1054      &          rescale_mode
1055               call flush(iout)
1056               return1
1057             endif
1058 c            write (iout,*) "ftprim",ftprim
1059 c            write (iout,*) "ftbis",ftbis
1060             betaT=1.0d0/(1.987D-3*betaT)
1061             if (betaT.ge.beta_h(1,iparm)) then
1062               potEmin=potEmin_all(1,iparm)+
1063      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1064      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1065      &         (1.0/betaT-1.0/beta_h(1,iparm))
1066 #ifdef DEBUG
1067               write(iout,*) "first",temper,potEmin
1068 #endif
1069             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1070               potEmin=potEmin_all(nT_h(iparm),iparm)+
1071      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1072      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1073      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1074 #ifdef DEBUG
1075               write (iout,*) "last",temper,potEmin
1076 #endif
1077             else
1078               do l=1,nT_h(iparm)-1
1079                 if (betaT.le.beta_h(l,iparm) .and.
1080      &              betaT.gt.beta_h(l+1,iparm)) then
1081                   potEmin=potEmin_all(l,iparm)
1082 #ifdef DEBUG
1083                   write (iout,*) "l",l,
1084      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1085      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1086 #endif
1087                   exit
1088                 endif
1089               enddo
1090             endif
1091 #ifdef SPLITELE
1092       if (shield_mode.gt.0) then
1093             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1094      &      +ft(1)*welec*ees
1095      &      +ft(1)*wvdwpp*evdw1
1096      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1097      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1098      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1099      &      +ft(2)*wturn3*eello_turn3
1100      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1101      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1102      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1103             eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1104 C     &            +ftprim(6)*evdw_t
1105      &            +ftprim(1)*wscp*evdw2
1106      &            +ftprim(1)*welec*ees
1107      &            +ftprim(1)*wvdwpp*evdw1
1108      &            +ftprim(1)*wtor*etors+
1109      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1110      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1111      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1112      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1113      &            ftprim(1)*wsccor*esccor
1114             ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1115      &            +ftbis(1)*wscp*evdw2+
1116      &            ftbis(1)*welec*ees
1117      &            +ftbis(1)*wvdwpp*evdw
1118      &            +ftbis(1)*wtor*etors+
1119      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1120      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1121      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1122      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1123      &            ftbis(1)*wsccor*esccor
1124         else
1125             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1126      &      +wvdwpp*evdw1
1127      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1128      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1129      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1130      &      +ft(2)*wturn3*eello_turn3
1131      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1132      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1133      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1134             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1135      &            +ftprim(1)*wtor*etors+
1136      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1137      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1138      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1139      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1140      &            ftprim(1)*wsccor*esccor
1141             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1142      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1143      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1144      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1145      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1146      &            ftbis(1)*wsccor*esccor
1147       endif
1148 #else
1149       if (shield_mode.gt.0) then
1150             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
1151      &      +ft(1)*welec*(ees+evdw1)
1152      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1153      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1154      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1155      &      +ft(2)*wturn3*eello_turn3
1156      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1157      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1158      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1159             eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
1160      &           +ftprim(1)*welec*(ees+evdw1)
1161      &           +ftprim(1)*wtor*etors+
1162      &            ftprim(1)*wscp*evdw2+
1163      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1164      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1165      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1166      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1167      &            ftprim(1)*wsccor*esccor
1168             ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
1169      &           +ftbis(1)*wscp*evdw2
1170      &           +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1171      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1172      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1173      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1174      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1175      &            ftbis(1)*wsccor*esccor
1176        else
1177             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1178      &      +ft(1)*welec*(ees+evdw1)
1179      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1180      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1181      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1182      &      +ft(2)*wturn3*eello_turn3
1183      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1184      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1185      &      +wbond*estr+wliptran*eliptran+wtube*Etube
1186             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1187      &           +ftprim(1)*wtor*etors+
1188      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1189      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1190      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1191      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1192      &            ftprim(1)*wsccor*esccor
1193             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1194      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1195      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1196      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1197      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1198      &            ftbis(1)*wsccor*esccor
1199
1200        endif
1201
1202 #endif
1203             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1204 #ifdef DEBUG
1205             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
1206      &        " etot",etot," entfac",entfac(t),
1207      &        " weight",weight," ebis",ebis
1208 #endif
1209             etot=etot-temper*eprim
1210 #ifdef MPI
1211             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1212             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1213             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1214             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1215             do j=1,nQ+2
1216               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1217               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1218               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1219      &         +etot*q(j,t)*weight
1220             enddo
1221 #else
1222             sumW(k,iparm)=sumW(k,iparm)+weight
1223             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1224             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1225             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1226             do j=1,nQ+2
1227               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1228               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1229               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1230      &         +etot*q(j,t)*weight
1231             enddo
1232 #endif
1233           enddo
1234           indE = aint(potE(t,iparm)-aint(potEmin))
1235           if (indE.ge.0 .and. indE.le.maxinde) then
1236             if (indE.gt.upindE_p) upindE_p=indE
1237             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1238           endif
1239 #ifdef MPI
1240           do ib=1,nT_h(iparm)
1241             potEmin=potEmin_all(ib,iparm)
1242             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1243             hfin_p(ind,ib)=hfin_p(ind,ib)+
1244      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1245              if (rmsrgymap) then
1246                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1247                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1248                hrmsrgy_p(indrgy,indrms,ib)=
1249      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1250              endif
1251           enddo
1252 #else
1253           do ib=1,nT_h(iparm)
1254             potEmin=potEmin_all(ib,iparm)
1255             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1256             hfin(ind,ib)=hfin(ind,ib)+
1257      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1258              if (rmsrgymap) then
1259                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1260                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1261                hrmsrgy(indrgy,indrms,ib)=
1262      &           hrmsrgy(indrgy,indrms,ib)+expfac
1263              endif
1264           enddo
1265 #endif
1266         enddo ! t
1267         do ib=1,nT_h(iparm)
1268           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1269      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1270           if (rmsrgymap) then
1271           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1272      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1273      &       WHAM_COMM,IERROR)
1274           endif
1275         enddo
1276         call MPI_Reduce(upindE_p,upindE,1,
1277      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1278         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1279      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1280
1281         if (me1.eq.master) then
1282
1283         if (histout) then
1284
1285         write (iout,'(6x,$)')
1286         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1287      &   ib=1,nT_h(iparm))
1288         write (iout,*)
1289
1290         write (iout,'(/a)') 'Final histograms'
1291         if (histfile) then
1292           if (nslice.eq.1) then
1293             if (separate_parset) then
1294               write(licz3,"(bz,i3.3)") myparm
1295               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1296             else
1297               histname=prefix(:ilen(prefix))//'.hist'
1298             endif
1299           else
1300             if (separate_parset) then
1301               write(licz3,"(bz,i3.3)") myparm
1302               histname=prefix(:ilen(prefix))//'_par'//licz3//
1303      &         '_slice_'//licz2//'.hist'
1304             else
1305               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1306             endif
1307           endif
1308 #if defined(AIX) || defined(PGI)
1309           open (ihist,file=histname,position='append')
1310 #else
1311           open (ihist,file=histname,access='append')
1312 #endif
1313         endif
1314
1315         do t=0,tmax
1316           liczba=t
1317           sumH=0.0d0
1318           do ib=1,nT_h(iparm)
1319             sumH=sumH+hfin(t,ib)
1320           enddo
1321           if (sumH.gt.0.0d0) then
1322             do j=1,nQ
1323               jj = mod(liczba,nbin1)
1324               liczba=liczba/nbin1
1325               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1326               if (histfile)
1327      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1328             enddo
1329             do ib=1,nT_h(iparm)
1330               write (iout,'(e20.10,$)') hfin(t,ib)
1331               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1332             enddo
1333             write (iout,'(i5)') iparm
1334             if (histfile) write (ihist,'(i5)') iparm
1335           endif
1336         enddo
1337
1338         endif
1339
1340         if (entfile) then
1341           if (nslice.eq.1) then
1342             if (separate_parset) then
1343               write(licz3,"(bz,i3.3)") myparm
1344               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1345             else
1346               histname=prefix(:ilen(prefix))//'.ent'
1347             endif
1348           else
1349             if (separate_parset) then
1350               write(licz3,"(bz,i3.3)") myparm
1351               histname=prefix(:ilen(prefix))//'par_'//licz3//
1352      &           '_slice_'//licz2//'.ent'
1353             else
1354               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1355             endif
1356           endif
1357 #if defined(AIX) || defined(PGI)
1358           open (ihist,file=histname,position='append')
1359 #else
1360           open (ihist,file=histname,access='append')
1361 #endif
1362           write (ihist,'(a)') "# Microcanonical entropy"
1363           do i=0,upindE
1364             write (ihist,'(f8.0,$)') dint(potEmin)+i
1365             if (histE(i).gt.0.0e0) then
1366               write (ihist,'(f15.5,$)') dlog(histE(i))
1367             else
1368               write (ihist,'(f15.5,$)') 0.0d0
1369             endif
1370           enddo
1371           write (ihist,*)
1372           close(ihist)
1373         endif
1374         write (iout,*) "Microcanonical entropy"
1375         do i=0,upindE
1376           write (iout,'(f8.0,$)') dint(potEmin)+i
1377           if (histE(i).gt.0.0e0) then
1378             write (iout,'(f15.5,$)') dlog(histE(i))
1379           else
1380             write (iout,'(f15.5,$)') 0.0d0
1381           endif
1382           write (iout,*)
1383         enddo
1384         if (rmsrgymap) then
1385           if (nslice.eq.1) then
1386             if (separate_parset) then
1387               write(licz3,"(bz,i3.3)") myparm
1388               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1389             else
1390               histname=prefix(:ilen(prefix))//'.rmsrgy'
1391             endif
1392           else
1393             if (separate_parset) then
1394               write(licz3,"(bz,i3.3)") myparm
1395               histname=prefix(:ilen(prefix))//'_par'//licz3//
1396      &         '_slice_'//licz2//'.rmsrgy'
1397             else
1398              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1399             endif
1400           endif
1401 #if defined(AIX) || defined(PGI)
1402           open (ihist,file=histname,position='append')
1403 #else
1404           open (ihist,file=histname,access='append')
1405 #endif
1406           do i=0,nbin_rms
1407             do j=0,nbin_rgy
1408               write(ihist,'(2f8.2,$)')
1409      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1410               do ib=1,nT_h(iparm)
1411                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1412                   write(ihist,'(e14.5,$)')
1413      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1414      &              +potEmin
1415                 else
1416                   write(ihist,'(e14.5,$)') 1.0d6
1417                 endif
1418               enddo
1419               write (ihist,'(i2)') iparm
1420             enddo
1421           enddo
1422           close(ihist)
1423         endif
1424         endif
1425       enddo ! iparm
1426 #ifdef MPI
1427       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1428      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1429       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1430      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1431       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1432      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1433       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1434      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1435       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1436      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1437       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1438      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1439      &   WHAM_COMM,IERROR)
1440       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1441      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1442      &   WHAM_COMM,IERROR)
1443       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1444      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1445      &   WHAM_COMM,IERROR)
1446       if (me.eq.master) then
1447 #endif
1448       write (iout,'(/a)') 'Thermal characteristics of folding'
1449       if (nslice.eq.1) then
1450         nazwa=prefix
1451       else
1452         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1453       endif
1454       iln=ilen(nazwa)
1455       if (nparmset.eq.1 .and. .not.separate_parset) then
1456         nazwa=nazwa(:iln)//".thermal"
1457       else if (nparmset.eq.1 .and. separate_parset) then
1458         write(licz3,"(bz,i3.3)") myparm
1459         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1460       endif
1461       do iparm=1,nParmSet
1462       if (nparmset.gt.1) then
1463         write(licz3,"(bz,i3.3)") iparm
1464         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1465       endif
1466       open(34,file=nazwa)
1467       if (separate_parset) then
1468         write (iout,'(a,i3)') "Parameter set",myparm
1469       else
1470         write (iout,'(a,i3)') "Parameter set",iparm
1471       endif
1472       do i=0,NGridT
1473         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1474         if (betaT.ge.beta_h(1,iparm)) then
1475           potEmin=potEmin_all(1,iparm)
1476         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1477           potEmin=potEmin_all(nT_h(iparm),iparm)
1478         else
1479           do l=1,nT_h(iparm)-1
1480             if (betaT.le.beta_h(l,iparm) .and.
1481      &          betaT.gt.beta_h(l+1,iparm)) then
1482               potEmin=potEmin_all(l,iparm)
1483               exit
1484             endif
1485           enddo
1486         endif
1487
1488         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1489         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1490      &    sumW(i,iparm)
1491         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1492      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1493         do j=1,nQ+2
1494           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1495           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1496      &     -sumQ(j,i,iparm)**2
1497           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1498      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1499         enddo
1500         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1501      &   (startGridT+i*delta_T))+potEmin
1502         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1503      &   sumW(i,iparm),sumE(i,iparm)
1504         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1505         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1506      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1507         write (iout,*)
1508         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1509      &   sumW(i,iparm),sumE(i,iparm)
1510         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1511         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1512      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1513         write (34,*)
1514         call flush(34)
1515       enddo
1516       close(34)
1517       enddo
1518       if (histout) then
1519       do t=0,tmax
1520         if (hfin_ent(t).gt.0.0d0) then
1521           liczba=t
1522           jj = mod(liczba,nbin1)
1523           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1524      &     hfin_ent(t)
1525           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1526      &      dmin+(jj+0.5d0)*delta,
1527      &     hfin_ent(t)
1528         endif
1529       enddo
1530       if (histfile) close(ihist)
1531       endif
1532
1533 #ifdef ZSCORE
1534 ! Write data for zscore
1535       if (nslice.eq.1) then
1536         zscname=prefix(:ilen(prefix))//".zsc"
1537       else
1538         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1539       endif
1540 #if defined(AIX) || defined(PGI)
1541       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1542 #else
1543       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1544 #endif
1545       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1546       do iparm=1,nParmSet
1547         write (izsc,'("NT=",i1)') nT_h(iparm)
1548       do ib=1,nT_h(iparm)
1549         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1550      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1551         jj = min0(nR(ib,iparm),7)
1552         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1553         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1554         write (izsc,'("&")')
1555         if (nR(ib,iparm).gt.7) then
1556           do ii=8,nR(ib,iparm),9
1557             jj = min0(nR(ib,iparm),ii+8)
1558             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1559             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1560             write (izsc,'("&")')
1561           enddo
1562         endif
1563         write (izsc,'("FI=",$)')
1564         jj=min0(nR(ib,iparm),7)
1565         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1566         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1567         write (izsc,'("&")')
1568         if (nR(ib,iparm).gt.7) then
1569           do ii=8,nR(ib,iparm),9
1570             jj = min0(nR(ib,iparm),ii+8)
1571             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1572             if (jj.eq.nR(ib,iparm)) then
1573               write (izsc,*)
1574             else
1575               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1576               write (izsc,'(t80,"&")')
1577             endif
1578           enddo
1579         endif
1580         do i=1,nR(ib,iparm)
1581           write (izsc,'("KH=",$)')
1582           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1583           write (izsc,'(" Q0=",$)')
1584           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1585           write (izsc,*)
1586         enddo
1587       enddo
1588       enddo
1589       close(izsc)
1590 #endif
1591 #ifdef MPI
1592       endif
1593 #endif
1594
1595       return
1596 C#undef DEBUG
1597       end