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