source/wham/wham_calc.F90

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