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,ecation_nucl
 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             ecation_nucl= enetb(50,i,iparm)
 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+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+wcatnucl*ecation_nucl
 421
 422
 423
 424 #else
 425             etot=wsc*evdw+wscp*evdw2 &
 426             +welec*(ees+evdw1) &
 427             +wang*ebe+wtor*etors+wscloc*escloc &
 428             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
 429             +wcorr6*ecorr6+wturn4*eello_turn4 &
 430             +wturn3*eello_turn3 &
 431             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
 432             +wtor_d*etors_d+wsccor*esccor &
 433             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 434             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 435             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 436             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 437             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 438             +wscbase*escbase&
 439             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl
 440
 441
 442
 443 #endif
 444
 445 #ifdef DEBUG
 446             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
 447               etot,potEmin
 448 #endif
 449 #ifdef DEBUG
 450             if (iparm.eq.1 .and. ib.eq.1) then
 451             write (iout,*)"Conformation",i
 452             energia(0)=etot
 453             do k=1,max_ene
 454               energia(k)=enetb(k,i,iparm)
 455             enddo
 456 !            call enerprint(energia(0),fT)
 457             call enerprint(energia(0))
 458             endif
 459 #endif
 460             do kk=1,nR(ib,iparm)
 461               Econstr=0.0d0
 462               do j=1,nQ
 463                 ddW = q(j,i)
 464                 Econstr=Econstr+Kh(j,kk,ib,iparm) &
 465                  *(ddW-q0(j,kk,ib,iparm))**2
 466               enddo
 467               v(i,kk,ib,iparm)= &
 468                 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 469 #ifdef DEBUG
 470               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
 471                etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 472 #endif
 473             enddo ! kk
 474           enddo   ! ib
 475         enddo     ! iparm
 476       enddo       ! i
 477 ! Simple iteration to calculate free energies corresponding to all simulation
 478 ! runs.
 479       do iter=1,maxit
 480
 481 ! Compute new free-energy values corresponding to the righ-hand side of the
 482 ! equation and their derivatives.
 483         write (iout,*) "------------------------fi"
 484 #ifdef MPI
 485         do t=1,scount(me1)
 486 #else
 487         do t=1,ntot(islice)
 488 #endif
 489           vmax=-1.0d+20
 490           do i=1,nParmSet
 491             do k=1,nT_h(i)
 492               do l=1,nR(k,i)
 493                 vf=v(t,l,k,i)+f(l,k,i)
 494                 if (vf.gt.vmax) vmax=vf
 495               enddo
 496             enddo
 497           enddo
 498           denom=0.0d0
 499           do i=1,nParmSet
 500             do k=1,nT_h(i)
 501               do l=1,nR(k,i)
 502                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 503                 if (aux.gt.-200.0d0) &
 504                 denom=denom+snk(l,k,i,islice)*dexp(aux)
 505               enddo
 506             enddo
 507           enddo
 508           entfac(t)=-dlog(denom)-vmax
 509           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 510 #ifdef DEBUG
 511           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 512 #endif
 513         enddo
 514         do iparm=1,nParmSet
 515           do iib=1,nT_h(iparm)
 516             do ii=1,nR(iib,iparm)
 517 #ifdef MPI
 518               fi_p_min(ii,iib,iparm)=-1.0d5
 519               do t=1,scount(me)
 520 !                fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 521 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 522                 aux=v(t,ii,iib,iparm)+entfac(t)
 523                 if (aux.gt.fi_p_min(ii,iib,iparm))   fi_p_min(ii,iib,iparm)=aux
 524
 525 !#define DEBUG
 526 #ifdef DEBUG
 527               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
 528                v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
 529 #endif
 530 !#undef DEBUG
 531               enddo
 532 #else
 533 !              fi(ii,iib,iparm)=0.0d0
 534               do t=1,ntot(islice)
 535 !                fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 536 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 537                 aux=v(t,ii,iib,iparm)+entfac(t)
 538                 if (aux.gt.fi_min(ii,iib,iparm))
 539      &                                     fi_min(ii,iib,iparm)=aux
 540
 541               enddo
 542 #endif
 543             enddo ! ii
 544           enddo ! iib
 545         enddo ! iparm
 546
 547 #ifdef MPI
 548 #ifdef DEBUG
 549         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 550         do iparm=1,nParmSet
 551           do ib=1,nT_h(nparmset)
 552             write (iout,*) "iparm",iparm," ib",ib
 553             write (iout,*) "beta=",beta_h(ib,iparm)
 554             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 555             write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
 556           enddo
 557         enddo
 558 #endif
 559         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
 560               MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 561
 562 #ifdef DEBUG
 563         write (iout,*) "fi_min after AllReduce"
 564         do i=1,nParmSet
 565           do j=1,nT_h(i)
 566             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 567           enddo
 568         enddo
 569 #endif
 570 !#undef DEBUG
 571 #endif
 572         do iparm=1,nParmSet
 573           do iib=1,nT_h(iparm)
 574             do ii=1,nR(iib,iparm)
 575 #ifdef MPI
 576               fi_p(ii,iib,iparm)=0.0d0
 577               do t=1,scount(me)
 578                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 579                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 580 #ifdef DEBUG
 581               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
 582               v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
 583               fi_p(ii,iib,iparm)
 584 #endif
 585               enddo
 586 #else
 587               fi(ii,iib,iparm)=0.0d0
 588               do t=1,ntot(islice)
 589                 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 590                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 591               enddo
 592 #endif
 593             enddo ! ii
 594           enddo ! iib
 595         enddo ! iparm
 596
 597 #ifdef MPI
 598 #ifdef DEBUG
 599         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 600         do iparm=1,nParmSet
 601           do ib=1,nT_h(nparmset)
 602             write (iout,*) "iparm",iparm," ib",ib
 603             write (iout,*) "beta=",beta_h(ib,iparm)
 604             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 605           enddo
 606         enddo
 607 #endif
 608 #ifdef DEBUG
 609         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
 610         maxR*MaxT_h*nParmSet
 611         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
 612         " WHAM_COMM",WHAM_COMM
 613 #endif
 614
 615         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
 616          maxR*MaxT_h*nParmSet
 617         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
 618          " WHAM_COMM",WHAM_COMM
 619         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
 620          MPI_DOUBLE_PRECISION,&
 621          MPI_SUM,Master,WHAM_COMM,IERROR)
 622 #ifdef DEBUG
 623         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 624         do iparm=1,nParmSet
 625           write (iout,*) "iparm",iparm
 626           do ib=1,nT_h(iparm)
 627             write (iout,*) "beta=",beta_h(ib,iparm)
 628             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 629           enddo
 630         enddo
 631 #endif
 632         if (me1.eq.Master) then
 633 #endif
 634         avefi=0.0d0
 635         do iparm=1,nParmSet
 636           do ib=1,nT_h(iparm)
 637             do i=1,nR(ib,iparm)
 638               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 639               avefi=avefi+fi(i,ib,iparm)
 640             enddo
 641           enddo
 642         enddo
 643         avefi=avefi/nFi
 644         do iparm=1,nParmSet
 645           write (iout,*) "Parameter set",iparm
 646           do ib =1,nT_h(iparm)
 647             write (iout,*) "beta=",beta_h(ib,iparm)
 648             do i=1,nR(ib,iparm)
 649               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 650             enddo
 651             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 652             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 653           enddo
 654         enddo
 655
 656 ! Compute the norm of free-energy increments.
 657         finorm=0.0d0
 658         do iparm=1,nParmSet
 659           do ib=1,nT_h(iparm)
 660             do i=1,nR(ib,iparm)
 661               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 662               f(i,ib,iparm)=fi(i,ib,iparm)
 663             enddo
 664           enddo
 665         enddo
 666
 667         write (iout,*) 'Iteration',iter,' finorm',finorm
 668
 669 #ifdef MPI
 670         endif
 671         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
 672          MPI_DOUBLE_PRECISION,Master,&
 673          WHAM_COMM,IERROR)
 674         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
 675          WHAM_COMM,IERROR)
 676 #endif
 677 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 678         if (finorm.lt.fimin) then
 679           write (iout,*) 'Iteration converged'
 680           goto 20
 681         endif
 682
 683       enddo ! iter
 684
 685    20 continue
 686 ! Now, put together the histograms from all simulations, in order to get the
 687 ! unbiased total histogram.
 688 !C Determine the minimum free energies
 689 #ifdef MPI
 690       do i=1,scount(me1)
 691 #else
 692       do i=1,ntot(islice)
 693 #endif
 694 !c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 695         do iparm=1,nParmSet
 696 #ifdef DEBUG
 697           write (iout,'(2i5,21f8.2)') i,iparm, &
 698           (enetb(k,i,iparm),k=1,26)
 699 #endif
 700           call restore_parm(iparm)
 701 #ifdef DEBUG
 702           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
 703            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
 704            wtor_d,wsccor,wbond
 705 #endif
 706           do ib=1,nT_h(iparm)
 707             if (rescale_modeW.eq.1) then
 708               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 709               quotl=1.0d0
 710               kfacl=1.0d0
 711               do l=1,5
 712                 quotl1=quotl
 713                 quotl=quotl*quot
 714                 kfacl=kfacl*kfac
 715                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 716               enddo
 717 #if defined(FUNCTH)
 718               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 719               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 720 #elif defined(FUNCT)
 721               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 722 #else
 723               ft(6)=1.0d0
 724 #endif
 725             else if (rescale_modeW.eq.2) then
 726               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 727               quotl=1.0d0
 728               do l=1,5
 729                 quotl=quotl*quot
 730                 fT(l)=1.12692801104297249644d0/ &
 731                   dlog(dexp(quotl)+dexp(-quotl))
 732                enddo
 733 #if defined(FUNCTH)
 734               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 735               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 736 #elif defined(FUNCT)
 737               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 738 #else
 739               ft(6)=1.0d0
 740 #endif
 741 !c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 742             else if (rescale_modeW.eq.0) then
 743               do l=1,6
 744                 fT(l)=1.0d0
 745               enddo
 746             else
 747               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
 748                rescale_mode
 749               call flush(iout)
 750               return 1
 751             endif
 752             evdw=enetb(1,i,iparm)
 753 !            evdw_t=enetb(21,i,iparm)
 754             evdw_t=enetb(20,i,iparm)
 755 #ifdef SCP14
 756 !            evdw2_14=enetb(17,i,iparm)
 757             evdw2_14=enetb(18,i,iparm)
 758             evdw2=enetb(2,i,iparm)+evdw2_14
 759 #else
 760             evdw2=enetb(2,i,iparm)
 761             evdw2_14=0.0d0
 762 #endif
 763 #ifdef SPLITELE
 764             ees=enetb(3,i,iparm)
 765             evdw1=enetb(16,i,iparm)
 766 #else
 767             ees=enetb(3,i,iparm)
 768             evdw1=0.0d0
 769 #endif
 770             ecorr=enetb(4,i,iparm)
 771             ecorr5=enetb(5,i,iparm)
 772             ecorr6=enetb(6,i,iparm)
 773             eel_loc=enetb(7,i,iparm)
 774             eello_turn3=enetb(8,i,iparm)
 775             eello_turn4=enetb(9,i,iparm)
 776             eello_turn6=enetb(10,i,iparm)
 777             ebe=enetb(11,i,iparm)
 778             escloc=enetb(12,i,iparm)
 779             etors=enetb(13,i,iparm)
 780             etors_d=enetb(14,i,iparm)
 781             ehpb=enetb(15,i,iparm)
 782 !            estr=enetb(18,i,iparm)
 783             estr=enetb(17,i,iparm)
 784 !            esccor=enetb(19,i,iparm)
 785             esccor=enetb(21,i,iparm)
 786 !            edihcnstr=enetb(20,i,iparm)
 787             edihcnstr=enetb(19,i,iparm)
 788 !            ehomology_constr=enetb(22,i,iparm)
 789 !            esaxs=enetb(26,i,iparm)
 790           ecationcation=enetb(41,i,iparm)
 791           ecation_prot=enetb(42,i,iparm)
 792             evdwpp =      enetb(26,i,iparm)
 793             eespp =      enetb(27,i,iparm)
 794             evdwpsb =      enetb(28,i,iparm)
 795             eelpsb =      enetb(29,i,iparm)
 796             evdwsb =      enetb(30,i,iparm)
 797             eelsb =      enetb(31,i,iparm)
 798             estr_nucl =      enetb(32,i,iparm)
 799             ebe_nucl =      enetb(33,i,iparm)
 800             esbloc =      enetb(34,i,iparm)
 801             etors_nucl =      enetb(35,i,iparm)
 802             etors_d_nucl =      enetb(36,i,iparm)
 803             ecorr_nucl =      enetb(37,i,iparm)
 804             ecorr3_nucl =      enetb(38,i,iparm)
 805             epeppho=   enetb(49,i,iparm)
 806             escpho=    enetb(48,i,iparm)
 807             epepbase=  enetb(47,i,iparm)
 808             escbase=   enetb(46,i,iparm)
 809             ecation_nucl=enetb(50,i,iparm)
 810 #ifdef SPLITELE
 811             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 812             +wvdwpp*evdw1 &
 813             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 814             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 815             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 816             +ft(2)*wturn3*eello_turn3 &
 817             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
 818             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 819             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
 820             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 821             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 822             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 823             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 824             +wcorr3_nucl*ecorr3_nucl&
 825             +wscbase*escbase&
 826             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl
 827
 828 #else
 829             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
 830             +ft(1)*welec*(ees+evdw1) &
 831             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 832             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 833             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 834             +ft(2)*wturn3*eello_turn3 &
 835             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 836             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 837             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 838             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 839             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 840             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 841             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 842             +wcorr3_nucl*ecorr3_nucl&
 843             +wscbase*escbase&
 844             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl
 845
 846
 847 #endif
 848             write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
 849             etot=etot-entfac(i)/beta_h(ib,iparm)
 850             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 851
 852           enddo   ! ib
 853         enddo     ! iparm
 854       enddo       ! i
 855 #ifdef DEBUG
 856       write (iout,*) "The potEmin array before reduction"
 857       do i=1,nParmSet
 858         write (iout,*) "Parameter set",i
 859         do j=1,nT_h(i)
 860           write (iout,*) j,PotEmin_all(j,i)
 861         enddo
 862       enddo
 863       write (iout,*) "potEmin_min",potEmin_min
 864 #endif
 865 #ifdef MPI
 866 !C Determine the minimum energes for all parameter sets and temperatures
 867       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
 868        maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 869       do i=1,nParmSet
 870         do j=1,nT_h(i)
 871           potEmin_all(j,i)=potEmin_t_all(j,i)
 872         enddo
 873       enddo
 874 #endif
 875       potEmin_min=potEmin_all(1,1)
 876       do i=1,nParmSet
 877         do j=1,nT_h(i)
 878           if (potEmin_all(j,i).lt.potEmin_min) &
 879                   potEmin_min=potEmin_all(j,i)
 880         enddo
 881       enddo
 882 #ifdef DEBUG
 883       write (iout,*) "The potEmin array"
 884       do i=1,nParmSet
 885         write (iout,*) "Parameter set",i
 886         do j=1,nT_h(i)
 887           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
 888            PotEmin_all(j,i)
 889         enddo
 890       enddo
 891       write (iout,*) "potEmin_min",potEmin_min
 892 #endif
 893
 894 #ifdef MPI
 895       do t=0,tmax
 896         hfin_ent_p(t)=0.0d0
 897       enddo
 898 #else
 899       do t=0,tmax
 900         hfin_ent(t)=0.0d0
 901       enddo
 902 #endif
 903       write (iout,*) "--------------hist"
 904 #ifdef MPI
 905       do iparm=1,nParmSet
 906         do i=0,nGridT
 907           sumW_p(i,iparm)=0.0d0
 908           sumE_p(i,iparm)=0.0d0
 909           sumEbis_p(i,iparm)=0.0d0
 910           sumEsq_p(i,iparm)=0.0d0
 911           do j=1,nQ+2
 912             sumQ_p(j,i,iparm)=0.0d0
 913             sumQsq_p(j,i,iparm)=0.0d0
 914             sumEQ_p(j,i,iparm)=0.0d0
 915           enddo
 916         enddo
 917       enddo
 918       upindE_p=0
 919 #else
 920       do iparm=1,nParmSet
 921         do i=0,nGridT
 922           sumW(i,iparm)=0.0d0
 923           sumE(i,iparm)=0.0d0
 924           sumEbis(i,iparm)=0.0d0
 925           sumEsq(i,iparm)=0.0d0
 926           do j=1,nQ+2
 927             sumQ(j,i,iparm)=0.0d0
 928             sumQsq(j,i,iparm)=0.0d0
 929             sumEQ(j,i,iparm)=0.0d0
 930           enddo
 931         enddo
 932       enddo
 933       upindE=0
 934 #endif
 935 ! 8/26/05 entropy distribution
 936 !#ifdef MPI
 937 !      entmin_p=1.0d10
 938 !      entmax_p=-1.0d10
 939 !      do t=1,scount(me1)
 940 !!        ent=-dlog(entfac(t))
 941 !        ent=entfac(t)
 942 !        if (ent.lt.entmin_p) entmin_p=ent
 943 !        if (ent.gt.entmax_p) entmax_p=ent
 944 !      enddo
 945 !      write (iout,*) "entmin",entmin_p," entmax",entmax_p
 946 !!      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
 947 !      call flush(iout)
 948 !      call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
 949 !        WHAM_COMM,IERROR)
 950 !      call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
 951 !        WHAM_COMM,IERROR)
 952 !      write (iout,*) "entmin",entmin," entmax",entmax
 953 !      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
 954 !      ientmax=entmax-entmin
 955 !iientmax=entmax-entmin !el
 956 !write (iout,*) "ientmax",ientmax,entmax,entmin
 957 !write (iout,*) "iientmax",iientmax
 958 !      if (ientmax.gt.2000) ientmax=2000
 959 !      write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 960 !      call flush(iout)
 961 !      do t=1,scount(me1)
 962 !!        ient=-dlog(entfac(t))-entmin
 963 !        ient=entfac(t)-entmin
 964 !        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 965 !      enddo
 966 !      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
 967 !        MPI_SUM,WHAM_COMM,IERROR)
 968 !      if (me1.eq.Master) then
 969 !        write (iout,*) "Entropy histogram"
 970 !        do i=0,ientmax
 971 !          write(iout,'(f15.4,i10)') entmin+i,histent(i)
 972 !        enddo
 973 !      endif
 974 !#else
 975 !      entmin=1.0d10
 976 !      entmax=-1.0d10
 977 !      do t=1,ntot(islice)
 978 !        ent=entfac(t)
 979 !        if (ent.lt.entmin) entmin=ent
 980 !        if (ent.gt.entmax) entmax=ent
 981 !      enddo
 982 !      ientmax=-dlog(entmax)-entmin
 983 !      if (ientmax.gt.2000) ientmax=2000
 984 !      do t=1,ntot(islice)
 985 !        ient=entfac(t)-entmin
 986 !        if (ient.le.2000) histent(ient)=histent(ient)+1
 987 !      enddo
 988 !      write (iout,*) "Entropy histogram"
 989 !      do i=0,ientmax
 990 !        write(iout,'(2f15.4)') entmin+i,histent(i)
 991 !      enddo
 992 !#endif
 993
 994 #ifdef MPI
 995       write (iout,*) "me1",me1," scount",scount(me1) !d
 996
 997       do iparm=1,nParmSet
 998
 999 #ifdef MPI
1000         do ib=1,nT_h(iparm)
1001           do t=0,tmax
1002             hfin_p(t,ib)=0.0d0
1003           enddo
1004         enddo
1005         do i=1,maxindE
1006           histE_p(i)=0.0d0
1007         enddo
1008 #else
1009         do ib=1,nT_h(iparm)
1010           do t=0,tmax
1011             hfin(t,ib)=0.0d0
1012           enddo
1013         enddo
1014         do i=1,maxindE
1015           histE(i)=0.0d0
1016         enddo
1017 #endif
1018         do ib=1,nT_h(iparm)
1019           do i=0,MaxBinRms
1020             do j=0,MaxBinRgy
1021               hrmsrgy(j,i,ib)=0.0d0
1022 #ifdef MPI
1023               hrmsrgy_p(j,i,ib)=0.0d0
1024 #endif
1025             enddo
1026           enddo
1027         enddo
1028
1029         do t=1,scount(me1)
1030 #else
1031         do t=1,ntot(islice)
1032 #endif
1033           ind = ind_point(t)
1034 #ifdef MPI
1035           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1036 #else
1037           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1038 #endif
1039 !          write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
1040           call restore_parm(iparm)
1041 !          evdw=enetb(21,t,iparm)
1042           evdw=enetb(20,t,iparm)
1043           evdw_t=enetb(1,t,iparm)
1044 #ifdef SCP14
1045 !          evdw2_14=enetb(17,t,iparm)
1046           evdw2_14=enetb(18,t,iparm)
1047           evdw2=enetb(2,t,iparm)+evdw2_14
1048 #else
1049           evdw2=enetb(2,t,iparm)
1050           evdw2_14=0.0d0
1051 #endif
1052 #ifdef SPLITELE
1053           ees=enetb(3,t,iparm)
1054           evdw1=enetb(16,t,iparm)
1055 #else
1056           ees=enetb(3,t,iparm)
1057           evdw1=0.0d0
1058 #endif
1059           ecorr=enetb(4,t,iparm)
1060           ecorr5=enetb(5,t,iparm)
1061           ecorr6=enetb(6,t,iparm)
1062           eel_loc=enetb(7,t,iparm)
1063           eello_turn3=enetb(8,t,iparm)
1064           eello_turn4=enetb(9,t,iparm)
1065           eello_turn6=enetb(10,t,iparm)
1066           ebe=enetb(11,t,iparm)
1067           escloc=enetb(12,t,iparm)
1068           etors=enetb(13,t,iparm)
1069           etors_d=enetb(14,t,iparm)
1070           ehpb=enetb(15,t,iparm)
1071 !          estr=enetb(18,t,iparm)
1072           estr=enetb(17,t,iparm)
1073 !          esccor=enetb(19,t,iparm)
1074           esccor=enetb(21,t,iparm)
1075 !          edihcnstr=enetb(20,t,iparm)
1076           edihcnstr=enetb(19,t,iparm)
1077           edihcnstr=0.0d0
1078           ecationcation=enetb(41,t,iparm)
1079           ecation_prot=enetb(42,t,iparm)
1080             evdwpp =      enetb(26,t,iparm)
1081             eespp =      enetb(27,t,iparm)
1082             evdwpsb =      enetb(28,t,iparm)
1083             eelpsb =      enetb(29,t,iparm)
1084             evdwsb =      enetb(30,t,iparm)
1085             eelsb =      enetb(31,t,iparm)
1086             estr_nucl =      enetb(32,t,iparm)
1087             ebe_nucl =      enetb(33,t,iparm)
1088             esbloc =      enetb(34,t,iparm)
1089             etors_nucl =      enetb(35,t,iparm)
1090             etors_d_nucl =      enetb(36,t,iparm)
1091             ecorr_nucl =      enetb(37,t,iparm)
1092             ecorr3_nucl =      enetb(38,t,iparm)
1093             epeppho=   enetb(49,t,iparm)
1094             escpho=    enetb(48,t,iparm)
1095             epepbase=  enetb(47,t,iparm)
1096             escbase=   enetb(46,t,iparm)
1097             ecation_nucl=enetb(50,t,iparm)
1098
1099           do k=0,nGridT
1100             betaT=startGridT+k*delta_T
1101             temper=betaT
1102 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
1103 !d            fT=T0/betaT
1104 !d            ft=2*T0/(T0+betaT)
1105             if (rescale_modeW.eq.1) then
1106               quot=betaT/T0
1107               quotl=1.0d0
1108               kfacl=1.0d0
1109               do l=1,5
1110                 quotl1=quotl
1111                 quotl=quotl*quot
1112                 kfacl=kfacl*kfac
1113                 denom=kfacl-1.0d0+quotl
1114                 fT(l)=kfacl/denom
1115                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1116                 ftbis(l)=l*kfacl*quotl1* &
1117                  (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1118               enddo
1119 #if defined(FUNCTH)
1120               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1121                         320.0d0
1122               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1123              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1124                     /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1125 #elif defined(FUNCT)
1126               fT(6)=betaT/T0
1127               ftprim(6)=1.0d0/T0
1128               ftbis(6)=0.0d0
1129 #else
1130               fT(6)=1.0d0
1131               ftprim(6)=0.0d0
1132               ftbis(6)=0.0d0
1133 #endif
1134             else if (rescale_modeW.eq.2) then
1135               quot=betaT/T0
1136               quotl=1.0d0
1137               do l=1,5
1138                 quotl1=quotl
1139                 quotl=quotl*quot
1140                 eplus=dexp(quotl)
1141                 eminus=dexp(-quotl)
1142                 logfac=1.0d0/dlog(eplus+eminus)
1143                 tanhT=(eplus-eminus)/(eplus+eminus)
1144                 fT(l)=1.12692801104297249644d0*logfac
1145                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1146                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1147                 2*l*quotl1/T0*logfac* &
1148                 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1149                 +ftprim(l)*tanhT)
1150               enddo
1151 #if defined(FUNCTH)
1152               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1153                        320.0d0
1154               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1155              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1156                      /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1157 #elif defined(FUNCT)
1158               fT(6)=betaT/T0
1159               ftprim(6)=1.0d0/T0
1160               ftbis(6)=0.0d0
1161 #else
1162               fT(6)=1.0d0
1163               ftprim(6)=0.0d0
1164               ftbis(6)=0.0d0
1165 #endif
1166             else if (rescale_modeW.eq.0) then
1167               do l=1,5
1168                 fT(l)=1.0d0
1169                 ftprim(l)=0.0d0
1170               enddo
1171             else
1172               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
1173                 rescale_modeW
1174               call flush(iout)
1175               return 1
1176             endif
1177 !            write (iout,*) "ftprim",ftprim
1178 !            write (iout,*) "ftbis",ftbis
1179             betaT=1.0d0/(1.987D-3*betaT)
1180 #ifdef SPLITELE
1181             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1182             +wvdwpp*evdw1 &
1183             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1184             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1185             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1186             +ft(2)*wturn3*eello_turn3 &
1187             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1188             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1189             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1190             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1191             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1192             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1193             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1194             +wcorr3_nucl*ecorr3_nucl&
1195             +wscbase*escbase&
1196             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1197             + wcatnucl*ecation_nucl
1198
1199
1200             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
1201                  +ftprim(1)*wtor*etors+ &
1202                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1203                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1204                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1205                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1206                  ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
1207                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1208             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
1209                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1210                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1211                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1212                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1213                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1214                  +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
1215 #else
1216             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1217             +ft(1)*welec*(ees+evdw1) &
1218             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1219             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1220             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1221             +ft(2)*wturn3*eello_turn3 &
1222             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1223             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1224             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1225             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1226             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1227             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1228             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1229             +wcorr3_nucl*ecorr3_nucl&
1230             +wscbase*escbase&
1231             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1232             + wcatnucl*ecation_nucl
1233
1234
1235             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1236                 +ftprim(1)*wtor*etors+ &
1237                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1238                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1239                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1240                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1241                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1242                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1243             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1244                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1245                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1246                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1247                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1248                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1249                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1250
1251 #endif
1252             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1253 #ifdef DEBUG
1254             write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
1255               " etot",etot," entfac",entfac(t),&
1256               " weight",weight," ebis",ebis
1257 #endif
1258             etot=etot-temper*eprim
1259 #ifdef MPI
1260             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1261             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1262             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1263             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1264             do j=1,nQ+2
1265               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1266               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1267               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
1268                +etot*q(j,t)*weight
1269             enddo
1270 #else
1271             sumW(k,iparm)=sumW(k,iparm)+weight
1272             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1273             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1274             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1275             do j=1,nQ+2
1276               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1277               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1278               sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
1279                +etot*q(j,t)*weight
1280             enddo
1281 #endif
1282           enddo
1283           indE = aint(potE(t,iparm)-aint(potEmin))
1284           if (indE.ge.0 .and. indE.le.maxinde) then
1285             if (indE.gt.upindE_p) upindE_p=indE
1286             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1287           endif
1288 #ifdef MPI
1289           do ib=1,nT_h(iparm)
1290             potEmin=potEmin_all(ib,iparm)
1291             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1292             hfin_p(ind,ib)=hfin_p(ind,ib)+ &
1293              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1294              if (rmsrgymap) then
1295                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1296                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1297                hrmsrgy_p(indrgy,indrms,ib)= &
1298                  hrmsrgy_p(indrgy,indrms,ib)+expfac
1299              endif
1300           enddo
1301 #else
1302           do ib=1,nT_h(iparm)
1303             potEmin=potEmin_all(ib,iparm)
1304             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1305             hfin(ind,ib)=hfin(ind,ib)+ &
1306              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1307              if (rmsrgymap) then
1308                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1309                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1310                hrmsrgy(indrgy,indrms,ib)= &
1311                  hrmsrgy(indrgy,indrms,ib)+expfac
1312              endif
1313           enddo
1314 #endif
1315         enddo ! t
1316         do ib=1,nT_h(iparm)
1317           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
1318           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1319           if (rmsrgymap) then
1320           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
1321          (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1322              WHAM_COMM,IERROR)
1323           endif
1324         enddo
1325         call MPI_Reduce(upindE_p,upindE,1,&
1326           MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1327         call MPI_Reduce(histE_p(0),histE(0),maxindE,&
1328           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1329
1330         if (me1.eq.master) then
1331
1332         if (histout) then
1333
1334         write (iout,'(6x,$)')
1335         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
1336          ib=1,nT_h(iparm))
1337         write (iout,*)
1338
1339         write (iout,'(/a)') 'Final histograms'
1340         if (histfile) then
1341           if (nslice.eq.1) then
1342             if (separate_parset) then
1343               write(licz3,"(bz,i3.3)") myparm
1344               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1345             else
1346               histname=prefix(:ilen(prefix))//'.hist'
1347             endif
1348           else
1349             if (separate_parset) then
1350               write(licz3,"(bz,i3.3)") myparm
1351               histname=prefix(:ilen(prefix))//'_par'//licz3// &
1352                '_slice_'//licz2//'.hist'
1353             else
1354               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1355             endif
1356           endif
1357 #if defined(AIX) || defined(PGI)
1358           open (ihist,file=histname,position='append')
1359 #else
1360           open (ihist,file=histname,access='append')
1361 #endif
1362         endif
1363
1364         do t=0,tmax
1365           liczbaW=t
1366           sumH=0.0d0
1367           do ib=1,nT_h(iparm)
1368             sumH=sumH+hfin(t,ib)
1369           enddo
1370           if (sumH.gt.0.0d0) then
1371             do j=1,nQ
1372               jj = mod(liczbaW,nbin1)
1373               liczbaW=liczbaW/nbin1
1374               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1375               if (histfile) &
1376                  write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1377             enddo
1378             do ib=1,nT_h(iparm)
1379               write (iout,'(e20.10,$)') hfin(t,ib)
1380               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1381             enddo
1382             write (iout,'(i5)') iparm
1383             if (histfile) write (ihist,'(i5)') iparm
1384           endif
1385         enddo
1386
1387         endif
1388
1389         if (entfile) then
1390           if (nslice.eq.1) then
1391             if (separate_parset) then
1392               write(licz3,"(bz,i3.3)") myparm
1393               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1394             else
1395               histname=prefix(:ilen(prefix))//'.ent'
1396             endif
1397           else
1398             if (separate_parset) then
1399               write(licz3,"(bz,i3.3)") myparm
1400               histname=prefix(:ilen(prefix))//'par_'//licz3// &
1401                  '_slice_'//licz2//'.ent'
1402             else
1403               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1404             endif
1405           endif
1406 #if defined(AIX) || defined(PGI)
1407           open (ihist,file=histname,position='append')
1408 #else
1409           open (ihist,file=histname,access='append')
1410 #endif
1411           write (ihist,'(a)') "# Microcanonical entropy"
1412           do i=0,upindE
1413             write (ihist,'(f8.0,$)') dint(potEmin)+i
1414             if (histE(i).gt.0.0e0) then
1415               write (ihist,'(f15.5,$)') dlog(histE(i))
1416             else
1417               write (ihist,'(f15.5,$)') 0.0d0
1418             endif
1419           enddo
1420           write (ihist,*)
1421           close(ihist)
1422         endif
1423         write (iout,*) "Microcanonical entropy"
1424         do i=0,upindE
1425           write (iout,'(f8.0,$)') dint(potEmin)+i
1426           if (histE(i).gt.0.0e0) then
1427             write (iout,'(f15.5,$)') dlog(histE(i))
1428           else
1429             write (iout,'(f15.5,$)') 0.0d0
1430           endif
1431           write (iout,*)
1432         enddo
1433         if (rmsrgymap) then
1434           if (nslice.eq.1) then
1435             if (separate_parset) then
1436               write(licz3,"(bz,i3.3)") myparm
1437               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1438             else
1439               histname=prefix(:ilen(prefix))//'.rmsrgy'
1440             endif
1441           else
1442             if (separate_parset) then
1443               write(licz3,"(bz,i3.3)") myparm
1444               histname=prefix(:ilen(prefix))//'_par'//licz3// &
1445                '_slice_'//licz2//'.rmsrgy'
1446             else
1447              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1448             endif
1449           endif
1450 #if defined(AIX) || defined(PGI)
1451           open (ihist,file=histname,position='append')
1452 #else
1453           open (ihist,file=histname,access='append')
1454 #endif
1455           do i=0,nbin_rms
1456             do j=0,nbin_rgy
1457               write(ihist,'(2f8.2,$)') &
1458                 rgymin+deltrgy*j,rmsmin+deltrms*i
1459               do ib=1,nT_h(iparm)
1460                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1461                   write(ihist,'(e14.5,$)') &
1462                     -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
1463                     +potEmin
1464                 else
1465                   write(ihist,'(e14.5,$)') 1.0d6
1466                 endif
1467               enddo
1468               write (ihist,'(i2)') iparm
1469             enddo
1470           enddo
1471           close(ihist)
1472         endif
1473         endif
1474       enddo ! iparm
1475 #ifdef MPI
1476       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
1477          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1478       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
1479          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1480       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
1481          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1482       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
1483          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1484       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
1485          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1486       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
1487          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1488          WHAM_COMM,IERROR)
1489       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
1490          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1491          WHAM_COMM,IERROR)
1492       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
1493          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1494          WHAM_COMM,IERROR)
1495       if (me.eq.master) then
1496 #endif
1497       write (iout,'(/a)') 'Thermal characteristics of folding'
1498       if (nslice.eq.1) then
1499         nazwa=prefix
1500       else
1501         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1502       endif
1503       iln=ilen(nazwa)
1504       if (nparmset.eq.1 .and. .not.separate_parset) then
1505         nazwa=nazwa(:iln)//".thermal"
1506       else if (nparmset.eq.1 .and. separate_parset) then
1507         write(licz3,"(bz,i3.3)") myparm
1508         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1509       endif
1510       do iparm=1,nParmSet
1511       if (nparmset.gt.1) then
1512         write(licz3,"(bz,i3.3)") iparm
1513         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1514       endif
1515       open(34,file=nazwa)
1516       if (separate_parset) then
1517         write (iout,'(a,i3)') "Parameter set",myparm
1518       else
1519         write (iout,'(a,i3)') "Parameter set",iparm
1520       endif
1521       do i=0,NGridT
1522         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1523         if (betaT.ge.beta_h(1,iparm)) then
1524           potEmin=potEmin_all(1,iparm)
1525         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1526           potEmin=potEmin_all(nT_h(iparm),iparm)
1527         else
1528           do l=1,nT_h(iparm)-1
1529             if (betaT.le.beta_h(l,iparm) .and. &
1530                betaT.gt.beta_h(l+1,iparm)) then
1531               potEmin=potEmin_all(l,iparm)
1532               exit
1533             endif
1534           enddo
1535         endif
1536
1537         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1538         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
1539           sumW(i,iparm)
1540         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
1541           -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1542         do j=1,nQ+2
1543           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1544           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
1545            -sumQ(j,i,iparm)**2
1546           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
1547            -sumQ(j,i,iparm)*sumE(i,iparm)
1548         enddo
1549         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
1550          (startGridT+i*delta_T))+potEmin
1551         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1552          sumW(i,iparm),sumE(i,iparm)
1553         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1554         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1555          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1556         write (iout,*)
1557         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1558          sumW(i,iparm),sumE(i,iparm)
1559         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1560         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1561          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1562         write (34,*)
1563         call flush(34)
1564       enddo
1565       close(34)
1566       enddo
1567       if (histout) then
1568       do t=0,tmax
1569         if (hfin_ent(t).gt.0.0d0) then
1570           liczbaW=t
1571           jj = mod(liczbaW,nbin1)
1572           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
1573            hfin_ent(t)
1574           if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
1575             dmin+(jj+0.5d0)*delta,&
1576            hfin_ent(t)
1577         endif
1578       enddo
1579       if (histfile) close(ihist)
1580       endif
1581
1582 #ifdef ZSCORE
1583 ! Write data for zscore
1584       if (nslice.eq.1) then
1585         zscname=prefix(:ilen(prefix))//".zsc"
1586       else
1587         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1588       endif
1589 #if defined(AIX) || defined(PGI)
1590       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1591 #else
1592       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1593 #endif
1594       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1595       do iparm=1,nParmSet
1596         write (izsc,'("NT=",i1)') nT_h(iparm)
1597       do ib=1,nT_h(iparm)
1598         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
1599           1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1600         jj = min0(nR(ib,iparm),7)
1601         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1602         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1603         write (izsc,'("&")')
1604         if (nR(ib,iparm).gt.7) then
1605           do ii=8,nR(ib,iparm),9
1606             jj = min0(nR(ib,iparm),ii+8)
1607             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1608             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1609             write (izsc,'("&")')
1610           enddo
1611         endif
1612         write (izsc,'("FI=",$)')
1613         jj=min0(nR(ib,iparm),7)
1614         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1615         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1616         write (izsc,'("&")')
1617         if (nR(ib,iparm).gt.7) then
1618           do ii=8,nR(ib,iparm),9
1619             jj = min0(nR(ib,iparm),ii+8)
1620             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1621             if (jj.eq.nR(ib,iparm)) then
1622               write (izsc,*)
1623             else
1624               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1625               write (izsc,'(t80,"&")')
1626             endif
1627           enddo
1628         endif
1629         do i=1,nR(ib,iparm)
1630           write (izsc,'("KH=",$)')
1631           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1632           write (izsc,'(" Q0=",$)')
1633           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1634           write (izsc,*)
1635         enddo
1636       enddo
1637       enddo
1638       close(izsc)
1639 #endif
1640 #ifdef MPI
1641       endif
1642 #endif
1643       return
1644       end subroutine WHAMCALC
1645 !-----------------------------------------------------------------------------
1646       end module wham_calc
1647