source/wham/src-HCD-5D/wham_calc1.F.safe

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