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