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