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