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