source/wham/wham_calc.F90

   1       module wham_calc
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use wham_data
   5 !
   6       use ene_calc
   7 #ifdef MPI
   8       use MPI_data
   9 !      include "COMMON.MPI"
  10 #endif
  11       implicit none
  12 !-----------------------------------------------------------------------------
  13 !
  14 !
  15 !-----------------------------------------------------------------------------
  16       contains
  17 !-----------------------------------------------------------------------------
  18
  19       subroutine WHAMCALC(islice,*)
  20 ! Weighed Histogram Analysis Method (WHAM) code
  21 ! Written by A. Liwo based on the work of Kumar et al.,
  22 ! J.Comput.Chem., 13, 1011 (1992)
  23 !
  24 ! 2/1/05 Multiple temperatures allowed.
  25 ! 2/2/05 Free energies calculated directly from data points
  26 !  acc. to Eq. (21) of Kumar et al.; final histograms also
  27 !  constructed based on this equation.
  28 ! 2/12/05 Multiple parameter sets included
  29 !
  30 ! 2/2/05 Parallel version
  31 !      use wham_data
  32 !      use io_units
  33       use names
  34       use io_base, only:ilen
  35       use energy_data
  36 #ifdef MPI
  37       include "mpif.h"
  38 #endif
  39 !      include "DIMENSIONS"
  40 !      include "DIMENSIONS.ZSCOPT"
  41 !      include "DIMENSIONS.FREE"
  42       integer,parameter :: NGridT=400
  43       integer,parameter :: MaxBinRms=100,MaxBinRgy=100
  44       integer,parameter :: MaxHdim=200
  45 !      parameter (MaxHdim=200000)
  46       integer,parameter :: maxinde=200
  47 #ifdef MPI
  48       integer :: ierror,errcode,status(MPI_STATUS_SIZE)
  49 #endif
  50 !      include "COMMON.CONTROL"
  51 !      include "COMMON.IOUNITS"
  52 !      include "COMMON.FREE"
  53 !      include "COMMON.ENERGIES"
  54 !      include "COMMON.FFIELD"
  55 !      include "COMMON.SBRIDGE"
  56 !      include "COMMON.PROT"
  57 !      include "COMMON.ENEPS"
  58       integer,parameter :: MaxPoint=MaxStr,&
  59               MaxPointProc=MaxStr_Proc
  60       real(kind=8),parameter :: finorm_max=1.0d0
  61       real(kind=8) :: potfac,expfac,vf
  62 !      real(kind=8) :: potfac,entmin,entmax,expfac,vf
  63       integer :: islice
  64       integer :: i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  65       integer :: start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,&
  66               nbin_rmsrgy,liczbaW,iparm,nFi,indrgy,indrms
  67 ! 4/17/17 AKS & AL: histent is obsolete
  68       integer :: htot(0:MaxHdim)!,histent(0:2000)
  69       real(kind=8) :: v(MaxPointProc,MaxR,MaxT_h,nParmSet)  !(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  70       real(kind=8) :: energia(0:n_ene)
  71 !el      real(kind=8) :: energia(0:max_ene)
  72 #ifdef MPI
  73       integer :: tmax_t,upindE_p
  74       real(kind=8) :: fi_p(MaxR,MaxT_h,nParmSet), &
  75                      fi_p_min(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
  76       real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
  77               sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
  78       real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
  79               sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
  80       real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
  81               hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
  82               hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  83       real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  84       real(kind=8) :: potEmin_t,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p
  85 !      integer :: histent_p(0:2000)
  86       logical :: lprint=.true.
  87 #endif
  88       real(kind=8) :: delta_T=1.0d0,iientmax
  89       real(kind=8) :: rgymin,rmsmin,rgymax,rmsmax
  90       real(kind=8),dimension(0:nGridT,nParmSet) :: sumW,sumE,&
  91               sumEsq,sumEprim,sumEbis !(0:NGridT,Max_Parm)
  92       real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ,&
  93               sumQsq,sumEQ !(MaxQ1,0:NGridT,Max_Parm)
  94       real(kind=8) :: betaT,weight,econstr
  95       real(kind=8) :: fi(MaxR,MaxT_h,nParmSet),& !(MaxR,maxT_h,Max_Parm)
  96               ddW,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,&
  97               hfin(0:MaxHdim,maxT_h),histE(0:maxindE),&
  98               hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),&
  99               potEmin,ent,&
 100               hfin_ent(0:MaxHdim),vmax,aux,fi_min(MaxR,maxT_h,nParmSet), &
 101               potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min
 102       real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
 103         eprim,ebis,temper,kfac=2.4d0,T0=300.0d0,startGridT=200.0d0,&
 104         eplus,eminus,logfac,tanhT,tt
 105       real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
 106         escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
 107         eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
 108         ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
 109         evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
 110         ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
 111         elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran
 112
 113
 114
 115       integer :: ind_point(maxpoint),upindE,indE
 116       character(len=16) :: plik
 117       character(len=1) :: licz1
 118       character(len=2) :: licz2
 119       character(len=3) :: licz3
 120       character(len=128) :: nazwa
 121 !      integer ilen
 122 !      external ilen
 123 !el      ientmax=0
 124 !el      ent=0.0d0
 125       write(licz2,'(bz,i2.2)') islice
 126       nbin1 = 1.0d0/delta
 127       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
 128         i2/80(1h-)//)') islice
 129       write (iout,*) "delta",delta," nbin1",nbin1
 130       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 131       call flush(iout)
 132       dmin=0.0d0
 133       tmax=0
 134       potEmin=1.0d10
 135       do i=1,nParmset
 136         do j=1,nT_h(i)
 137           potEmin_all(j,i)=1.0d15
 138         enddo
 139       enddo
 140       entfac_min=1.0d10
 141       rgymin=1.0d10
 142       rmsmin=1.0d10
 143       rgymax=0.0d0
 144       rmsmax=0.0d0
 145       do t=0,MaxN
 146         htot(t)=0
 147       enddo
 148 #ifdef MPI
 149       do i=1,scount(me1)
 150 #else
 151       do i=1,ntot(islice)
 152 #endif
 153         do j=1,nParmSet
 154           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 155         enddo
 156         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 157         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 158         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 159         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 160         ind_point(i)=0
 161         do j=nQ,1,-1
 162           ind=(q(j,i)-dmin+1.0d-8)/delta
 163           if (j.eq.1) then
 164             ind_point(i)=ind_point(i)+ind
 165           else
 166             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 167           endif
 168 !          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 169 !     &      ind_point(i)
 170           call flush(iout)
 171           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 172             write (iout,*) "Error - index exceeds range for point",i,&
 173             " q=",q(j,i)," ind",ind_point(i)
 174 #ifdef MPI
 175             write (iout,*) "Processor",me1
 176             call flush(iout)
 177             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 178 #endif
 179             stop
 180           endif
 181         enddo ! j
 182         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 183         htot(ind_point(i))=htot(ind_point(i))+1
 184 #ifdef DEBUG
 185         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
 186          " htot",htot(ind_point(i))
 187         call flush(iout)
 188 #endif
 189       enddo ! i
 190       call flush(iout)
 191
 192       nbin=nbin1**nQ-1
 193       write (iout,'(a)') "Numbers of counts in Q bins"
 194       do t=0,tmax
 195         if (htot(t).gt.0) then
 196         write (iout,'(i15,$)') t
 197         liczbaW=t
 198         do j=1,nQ
 199           jj = mod(liczbaW,nbin1)
 200           liczbaW=liczbaW/nbin1
 201           write (iout,'(i5,$)') jj
 202         enddo
 203         write (iout,'(i8)') htot(t)
 204         endif
 205       enddo
 206       do iparm=1,nParmSet
 207       write (iout,'(a,i3)') "Number of data points for parameter set",&
 208        iparm
 209       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
 210        ib=1,nT_h(iparm))
 211       write (iout,'(i8)') stot(islice)
 212       write (iout,'(a)')
 213       enddo
 214       call flush(iout)
 215
 216 #ifdef MPI
 217       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
 218         WHAM_COMM,IERROR)
 219       tmax=tmax_t
 220 !      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
 221 !        MPI_MIN,WHAM_COMM,IERROR) !????
 222       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
 223         MPI_MIN,WHAM_COMM,IERROR)
 224       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
 225         MPI_MAX,WHAM_COMM,IERROR)
 226       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,&
 227         MPI_MIN,WHAM_COMM,IERROR)
 228       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
 229         MPI_MAX,WHAM_COMM,IERROR)
 230 !      potEmin=potEmin_t !/2 try now??
 231       rgymin=rgymin_t
 232       rgymax=rgymax_t
 233       rmsmin=rmsmin_t
 234       rmsmax=rmsmax_t
 235       write (iout,*) "potEmin",potEmin
 236 #endif
 237       rmsmin=deltrms*dint(rmsmin/deltrms)
 238       rmsmax=deltrms*dint(rmsmax/deltrms)
 239       rgymin=deltrms*dint(rgymin/deltrgy)
 240       rgymax=deltrms*dint(rgymax/deltrgy)
 241       nbin_rms=(rmsmax-rmsmin)/deltrms
 242       nbin_rgy=(rgymax-rgymin)/deltrgy
 243       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
 244        " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 245       nFi=0
 246       do i=1,nParmSet
 247         do j=1,nT_h(i)
 248           nFi=nFi+nR(j,i)
 249         enddo
 250       enddo
 251       write (iout,*) "nFi",nFi
 252 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 253 ! simulations.
 254 #ifdef MPI
 255       do i=1,scount(me1)
 256 #else
 257       do i=1,ntot(islice)
 258 #endif
 259 !        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 260         do iparm=1,nParmSet
 261 #ifdef DEBUG
 262           write (iout,'(2i5,21f8.2)') i,iparm,&
 263            (enetb(k,i,iparm),k=1,21)
 264 #endif
 265           call restore_parm(iparm)
 266 #ifdef DEBUG
 267           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 268             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 269             wtor_d,wsccor,wbond,wcatcat
 270 #endif
 271           do ib=1,nT_h(iparm)
 272 !el old rascale weights
 273 !
 274 !            if (rescale_modeW.eq.1) then
 275 !              quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 276 !              quotl=1.0d0
 277 !              kfacl=1.0d0
 278 !              do l=1,5
 279 !                quotl1=quotl
 280 !                quotl=quotl*quot
 281 !                kfacl=kfacl*kfac
 282 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
 283 !              enddo
 284 !#if defined(FUNCTH)
 285 !              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 286 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 287 !#elif defined(FUNCT)
 288 !              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 289 !#else
 290 !              ft(6)=1.0d0
 291 !#endif
 292 !            else if (rescale_modeW.eq.2) then
 293 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 294 !              quotl=1.0d0
 295 !              do l=1,5
 296 !                quotl=quotl*quot
 297 !                fT(l)=1.12692801104297249644d0/ &
 298 !                   dlog(dexp(quotl)+dexp(-quotl))
 299 !              enddo
 300 !#if defined(FUNCTH)
 301 !              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 302 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 303 !#elif defined(FUNCT)
 304 !              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 305 !#else
 306 !              ft(6)=1.0d0
 307 !#endif
 308 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 309 !            else if (rescale_modeW.eq.0) then
 310 !              do l=1,6
 311 !                fT(l)=1.0d0
 312 !              enddo
 313 !            else
 314 !              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
 315 !                rescale_modeW
 316 !              call flush(iout)
 317 !              return 1
 318 !            endif
 319 ! el end old rescale weights
 320             call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
 321
 322 !            call etot(enetb(0,i,iparm))
 323             evdw=enetb(1,i,iparm)
 324 !            evdw_t=enetb(21,i,iparm)
 325             evdw_t=enetb(20,i,iparm)
 326 #ifdef SCP14
 327 !            evdw2_14=enetb(17,i,iparm)
 328             evdw2_14=enetb(18,i,iparm)
 329             evdw2=enetb(2,i,iparm)+evdw2_14
 330 #else
 331             evdw2=enetb(2,i,iparm)
 332             evdw2_14=0.0d0
 333 #endif
 334 #ifdef SPLITELE
 335             ees=enetb(3,i,iparm)
 336             evdw1=enetb(16,i,iparm)
 337 #else
 338             ees=enetb(3,i,iparm)
 339             evdw1=0.0d0
 340 #endif
 341             ecorr=enetb(4,i,iparm)
 342             ecorr5=enetb(5,i,iparm)
 343             ecorr6=enetb(6,i,iparm)
 344             eel_loc=enetb(7,i,iparm)
 345             eello_turn3=enetb(8,i,iparm)
 346             eello_turn4=enetb(9,i,iparm)
 347             eello_turn6=enetb(10,i,iparm)
 348             ebe=enetb(11,i,iparm)
 349             escloc=enetb(12,i,iparm)
 350             etors=enetb(13,i,iparm)
 351             etors_d=enetb(14,i,iparm)
 352             ehpb=enetb(15,i,iparm)
 353 !            estr=enetb(18,i,iparm)
 354             estr=enetb(17,i,iparm)
 355 !            esccor=enetb(19,i,iparm)
 356             esccor=enetb(21,i,iparm)
 357 !            edihcnstr=enetb(20,i,iparm)
 358             edihcnstr=enetb(19,i,iparm)
 359           ecationcation=enetb(41,i,iparm)
 360           ecation_prot=enetb(42,i,iparm)
 361             evdwpp =      enetb(26,i,iparm)
 362             eespp =      enetb(27,i,iparm)
 363             evdwpsb =      enetb(28,i,iparm)
 364             eelpsb =      enetb(29,i,iparm)
 365             evdwsb =      enetb(30,i,iparm)
 366             eelsb =      enetb(31,i,iparm)
 367             estr_nucl =      enetb(32,i,iparm)
 368             ebe_nucl =      enetb(33,i,iparm)
 369             esbloc =      enetb(34,i,iparm)
 370             etors_nucl =      enetb(35,i,iparm)
 371             etors_d_nucl =      enetb(36,i,iparm)
 372             ecorr_nucl =      enetb(37,i,iparm)
 373             ecorr3_nucl =      enetb(38,i,iparm)
 374             epeppho=   enetb(49,i,iparm)
 375             escpho=    enetb(48,i,iparm)
 376             epepbase=  enetb(47,i,iparm)
 377             escbase=   enetb(46,i,iparm)
 378             ecation_nucl= enetb(50,i,iparm)
 379             elipbond=enetb(52,i,iparm)
 380             elipang=enetb(53,i,iparm)
 381             eliplj=enetb(54,i,iparm)
 382             elipelec=enetb(55,i,iparm)
 383             ecat_prottran=enetb(56,i,iparm)
 384             ecation_protang=enetb(57,i,iparm)
 385             eliptran=enetb(22,i,iparm)
 386 #ifdef DEBUG
 387             write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
 388              evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
 389              etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation
 390 #endif
 391
 392 !#ifdef SPLITELE
 393 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 394 !            +wvdwpp*evdw1 &
 395 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 396 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 397 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 398 !            +ft(2)*wturn3*eello_turn3 &
 399 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc &
 400 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 401 !            +wbond*estr
 402 !#else
 403 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
 404 !            +ft(1)*welec*(ees+evdw1) &
 405 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 406 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 407 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 408 !            +ft(2)*wturn3*eello_turn3 &
 409 !            +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 410 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 411 !            +wbond*estr
 412 !#endif
 413
 414 #ifdef SPLITELE
 415             etot=wsc*evdw+wscp*evdw2+welec*ees &
 416             +wvdwpp*evdw1 &
 417             +wang*ebe+wtor*etors+wscloc*escloc &
 418             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
 419             +wcorr6*ecorr6+wturn4*eello_turn4 &
 420             +wturn3*eello_turn3 &
 421             +wturn6*eello_turn6+wel_loc*eel_loc &
 422             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
 423             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 424             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 425             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 426             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 427             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 428             +wscbase*escbase&
 429             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
 430             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 431             +wliptran*eliptran
 432
 433
 434 #else
 435             etot=wsc*evdw+wscp*evdw2 &
 436             +welec*(ees+evdw1) &
 437             +wang*ebe+wtor*etors+wscloc*escloc &
 438             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
 439             +wcorr6*ecorr6+wturn4*eello_turn4 &
 440             +wturn3*eello_turn3 &
 441             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
 442             +wtor_d*etors_d+wsccor*esccor &
 443             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 444             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 445             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 446             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 447             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 448             +wscbase*escbase&
 449             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
 450             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 451             +wliptran*eliptran
 452
 453
 454
 455
 456 #endif
 457
 458 #ifdef DEBUG
 459             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
 460               etot,potEmin
 461 #endif
 462 #ifdef DEBUG
 463             if (iparm.eq.1 .and. ib.eq.1) then
 464             write (iout,*)"Conformation",i
 465             energia(0)=etot
 466             do k=1,max_ene
 467               energia(k)=enetb(k,i,iparm)
 468             enddo
 469 !            call enerprint(energia(0),fT)
 470             call enerprint(energia(0))
 471             endif
 472 #endif
 473             do kk=1,nR(ib,iparm)
 474               Econstr=0.0d0
 475               do j=1,nQ
 476                 ddW = q(j,i)
 477                 Econstr=Econstr+Kh(j,kk,ib,iparm) &
 478                  *(ddW-q0(j,kk,ib,iparm))**2
 479               enddo
 480               v(i,kk,ib,iparm)= &
 481                 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 482 #ifdef DEBUG
 483               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
 484                etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 485 #endif
 486             enddo ! kk
 487           enddo   ! ib
 488         enddo     ! iparm
 489       enddo       ! i
 490 ! Simple iteration to calculate free energies corresponding to all simulation
 491 ! runs.
 492       do iter=1,maxit
 493
 494 ! Compute new free-energy values corresponding to the righ-hand side of the
 495 ! equation and their derivatives.
 496         write (iout,*) "------------------------fi"
 497 #ifdef MPI
 498         do t=1,scount(me1)
 499 #else
 500         do t=1,ntot(islice)
 501 #endif
 502           vmax=-1.0d+20
 503           do i=1,nParmSet
 504             do k=1,nT_h(i)
 505               do l=1,nR(k,i)
 506                 vf=v(t,l,k,i)+f(l,k,i)
 507                 if (vf.gt.vmax) vmax=vf
 508               enddo
 509             enddo
 510           enddo
 511           denom=0.0d0
 512           do i=1,nParmSet
 513             do k=1,nT_h(i)
 514               do l=1,nR(k,i)
 515                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 516                 if (aux.gt.-200.0d0) &
 517                 denom=denom+snk(l,k,i,islice)*dexp(aux)
 518               enddo
 519             enddo
 520           enddo
 521           entfac(t)=-dlog(denom)-vmax
 522           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 523 #ifdef DEBUG
 524           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 525 #endif
 526         enddo
 527         do iparm=1,nParmSet
 528           do iib=1,nT_h(iparm)
 529             do ii=1,nR(iib,iparm)
 530 #ifdef MPI
 531               fi_p_min(ii,iib,iparm)=-1.0d5
 532               do t=1,scount(me)
 533 !                fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 534 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 535                 aux=v(t,ii,iib,iparm)+entfac(t)
 536                 if (aux.gt.fi_p_min(ii,iib,iparm))   fi_p_min(ii,iib,iparm)=aux
 537
 538 !#define DEBUG
 539 #ifdef DEBUG
 540               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
 541                v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
 542 #endif
 543 !#undef DEBUG
 544               enddo
 545 #else
 546 !              fi(ii,iib,iparm)=0.0d0
 547               do t=1,ntot(islice)
 548 !                fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 549 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 550                 aux=v(t,ii,iib,iparm)+entfac(t)
 551                 if (aux.gt.fi_min(ii,iib,iparm))
 552      &                                     fi_min(ii,iib,iparm)=aux
 553
 554               enddo
 555 #endif
 556             enddo ! ii
 557           enddo ! iib
 558         enddo ! iparm
 559
 560 #ifdef MPI
 561 #ifdef DEBUG
 562         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 563         do iparm=1,nParmSet
 564           do ib=1,nT_h(nparmset)
 565             write (iout,*) "iparm",iparm," ib",ib
 566             write (iout,*) "beta=",beta_h(ib,iparm)
 567             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 568             write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
 569           enddo
 570         enddo
 571 #endif
 572         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
 573               MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 574
 575 #ifdef DEBUG
 576         write (iout,*) "fi_min after AllReduce"
 577         do i=1,nParmSet
 578           do j=1,nT_h(i)
 579             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 580           enddo
 581         enddo
 582 #endif
 583 !#undef DEBUG
 584 #endif
 585         do iparm=1,nParmSet
 586           do iib=1,nT_h(iparm)
 587             do ii=1,nR(iib,iparm)
 588 #ifdef MPI
 589               fi_p(ii,iib,iparm)=0.0d0
 590               do t=1,scount(me)
 591                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 592                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 593 #ifdef DEBUG
 594               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
 595               v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
 596               fi_p(ii,iib,iparm)
 597 #endif
 598               enddo
 599 #else
 600               fi(ii,iib,iparm)=0.0d0
 601               do t=1,ntot(islice)
 602                 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 603                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 604               enddo
 605 #endif
 606             enddo ! ii
 607           enddo ! iib
 608         enddo ! iparm
 609
 610 #ifdef MPI
 611 #ifdef DEBUG
 612         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 613         do iparm=1,nParmSet
 614           do ib=1,nT_h(nparmset)
 615             write (iout,*) "iparm",iparm," ib",ib
 616             write (iout,*) "beta=",beta_h(ib,iparm)
 617             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 618           enddo
 619         enddo
 620 #endif
 621 #ifdef DEBUG
 622         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
 623         maxR*MaxT_h*nParmSet
 624         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
 625         " WHAM_COMM",WHAM_COMM
 626 #endif
 627
 628         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
 629          maxR*MaxT_h*nParmSet
 630         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
 631          " WHAM_COMM",WHAM_COMM
 632         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
 633          MPI_DOUBLE_PRECISION,&
 634          MPI_SUM,Master,WHAM_COMM,IERROR)
 635 #ifdef DEBUG
 636         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 637         do iparm=1,nParmSet
 638           write (iout,*) "iparm",iparm
 639           do ib=1,nT_h(iparm)
 640             write (iout,*) "beta=",beta_h(ib,iparm)
 641             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 642           enddo
 643         enddo
 644 #endif
 645         if (me1.eq.Master) then
 646 #endif
 647         avefi=0.0d0
 648         do iparm=1,nParmSet
 649           do ib=1,nT_h(iparm)
 650             do i=1,nR(ib,iparm)
 651               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 652               avefi=avefi+fi(i,ib,iparm)
 653             enddo
 654           enddo
 655         enddo
 656         avefi=avefi/nFi
 657         do iparm=1,nParmSet
 658           write (iout,*) "Parameter set",iparm
 659           do ib =1,nT_h(iparm)
 660             write (iout,*) "beta=",beta_h(ib,iparm)
 661             do i=1,nR(ib,iparm)
 662               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 663             enddo
 664             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 665             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 666           enddo
 667         enddo
 668
 669 ! Compute the norm of free-energy increments.
 670         finorm=0.0d0
 671         do iparm=1,nParmSet
 672           do ib=1,nT_h(iparm)
 673             do i=1,nR(ib,iparm)
 674               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 675               f(i,ib,iparm)=fi(i,ib,iparm)
 676             enddo
 677           enddo
 678         enddo
 679
 680         write (iout,*) 'Iteration',iter,' finorm',finorm
 681
 682 #ifdef MPI
 683         endif
 684         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
 685          MPI_DOUBLE_PRECISION,Master,&
 686          WHAM_COMM,IERROR)
 687         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
 688          WHAM_COMM,IERROR)
 689 #endif
 690 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 691         if (finorm.lt.fimin) then
 692           write (iout,*) 'Iteration converged'
 693           goto 20
 694         endif
 695
 696       enddo ! iter
 697
 698    20 continue
 699 ! Now, put together the histograms from all simulations, in order to get the
 700 ! unbiased total histogram.
 701 !C Determine the minimum free energies
 702 #ifdef MPI
 703       do i=1,scount(me1)
 704 #else
 705       do i=1,ntot(islice)
 706 #endif
 707 !c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 708         do iparm=1,nParmSet
 709 #ifdef DEBUG
 710           write (iout,'(2i5,21f8.2)') i,iparm, &
 711           (enetb(k,i,iparm),k=1,26)
 712 #endif
 713           call restore_parm(iparm)
 714 #ifdef DEBUG
 715           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
 716            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
 717            wtor_d,wsccor,wbond
 718 #endif
 719           do ib=1,nT_h(iparm)
 720             if (rescale_modeW.eq.1) then
 721               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 722               quotl=1.0d0
 723               kfacl=1.0d0
 724               do l=1,5
 725                 quotl1=quotl
 726                 quotl=quotl*quot
 727                 kfacl=kfacl*kfac
 728                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 729               enddo
 730 #if defined(FUNCTH)
 731               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 732               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 733 #elif defined(FUNCT)
 734               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 735 #else
 736               ft(6)=1.0d0
 737 #endif
 738             else if (rescale_modeW.eq.2) then
 739               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 740               quotl=1.0d0
 741               do l=1,5
 742                 quotl=quotl*quot
 743                 fT(l)=1.12692801104297249644d0/ &
 744                   dlog(dexp(quotl)+dexp(-quotl))
 745                enddo
 746 #if defined(FUNCTH)
 747               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 748               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 749 #elif defined(FUNCT)
 750               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 751 #else
 752               ft(6)=1.0d0
 753 #endif
 754 !c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 755             else if (rescale_modeW.eq.0) then
 756               do l=1,6
 757                 fT(l)=1.0d0
 758               enddo
 759             else
 760               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
 761                rescale_mode
 762               call flush(iout)
 763               return 1
 764             endif
 765             evdw=enetb(1,i,iparm)
 766 !            evdw_t=enetb(21,i,iparm)
 767             evdw_t=enetb(20,i,iparm)
 768 #ifdef SCP14
 769 !            evdw2_14=enetb(17,i,iparm)
 770             evdw2_14=enetb(18,i,iparm)
 771             evdw2=enetb(2,i,iparm)+evdw2_14
 772 #else
 773             evdw2=enetb(2,i,iparm)
 774             evdw2_14=0.0d0
 775 #endif
 776 #ifdef SPLITELE
 777             ees=enetb(3,i,iparm)
 778             evdw1=enetb(16,i,iparm)
 779 #else
 780             ees=enetb(3,i,iparm)
 781             evdw1=0.0d0
 782 #endif
 783             ecorr=enetb(4,i,iparm)
 784             ecorr5=enetb(5,i,iparm)
 785             ecorr6=enetb(6,i,iparm)
 786             eel_loc=enetb(7,i,iparm)
 787             eello_turn3=enetb(8,i,iparm)
 788             eello_turn4=enetb(9,i,iparm)
 789             eello_turn6=enetb(10,i,iparm)
 790             ebe=enetb(11,i,iparm)
 791             escloc=enetb(12,i,iparm)
 792             etors=enetb(13,i,iparm)
 793             etors_d=enetb(14,i,iparm)
 794             ehpb=enetb(15,i,iparm)
 795 !            estr=enetb(18,i,iparm)
 796             estr=enetb(17,i,iparm)
 797 !            esccor=enetb(19,i,iparm)
 798             esccor=enetb(21,i,iparm)
 799 !            edihcnstr=enetb(20,i,iparm)
 800             edihcnstr=enetb(19,i,iparm)
 801 !            ehomology_constr=enetb(22,i,iparm)
 802 !            esaxs=enetb(26,i,iparm)
 803           ecationcation=enetb(41,i,iparm)
 804           ecation_prot=enetb(42,i,iparm)
 805             evdwpp =      enetb(26,i,iparm)
 806             eespp =      enetb(27,i,iparm)
 807             evdwpsb =      enetb(28,i,iparm)
 808             eelpsb =      enetb(29,i,iparm)
 809             evdwsb =      enetb(30,i,iparm)
 810             eelsb =      enetb(31,i,iparm)
 811             estr_nucl =      enetb(32,i,iparm)
 812             ebe_nucl =      enetb(33,i,iparm)
 813             esbloc =      enetb(34,i,iparm)
 814             etors_nucl =      enetb(35,i,iparm)
 815             etors_d_nucl =      enetb(36,i,iparm)
 816             ecorr_nucl =      enetb(37,i,iparm)
 817             ecorr3_nucl =      enetb(38,i,iparm)
 818             epeppho=   enetb(49,i,iparm)
 819             escpho=    enetb(48,i,iparm)
 820             epepbase=  enetb(47,i,iparm)
 821             escbase=   enetb(46,i,iparm)
 822             ecation_nucl=enetb(50,i,iparm)
 823             elipbond=enetb(52,i,iparm)
 824             elipang=enetb(53,i,iparm)
 825             eliplj=enetb(54,i,iparm)
 826             elipelec=enetb(55,i,iparm)
 827             ecat_prottran=enetb(56,i,iparm)
 828             ecation_protang=enetb(57,i,iparm)
 829             eliptran=enetb(22,i,iparm)
 830 #ifdef SPLITELE
 831             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 832             +wvdwpp*evdw1 &
 833             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 834             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 835             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 836             +ft(2)*wturn3*eello_turn3 &
 837             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
 838             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 839             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
 840             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 841             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 842             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 843             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 844             +wcorr3_nucl*ecorr3_nucl&
 845             +wscbase*escbase&
 846             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
 847             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 848             +wliptran*eliptran
 849
 850
 851 #else
 852             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
 853             +ft(1)*welec*(ees+evdw1) &
 854             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 855             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 856             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 857             +ft(2)*wturn3*eello_turn3 &
 858             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 859             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 860             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 861             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 862             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 863             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 864             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 865             +wcorr3_nucl*ecorr3_nucl&
 866             +wscbase*escbase&
 867             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
 868             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 869             +wliptran*eliptran
 870
 871
 872
 873 #endif
 874             write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
 875             etot=etot-entfac(i)/beta_h(ib,iparm)
 876             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 877
 878           enddo   ! ib
 879         enddo     ! iparm
 880       enddo       ! i
 881 #ifdef DEBUG
 882       write (iout,*) "The potEmin array before reduction"
 883       do i=1,nParmSet
 884         write (iout,*) "Parameter set",i
 885         do j=1,nT_h(i)
 886           write (iout,*) j,PotEmin_all(j,i)
 887         enddo
 888       enddo
 889       write (iout,*) "potEmin_min",potEmin_min
 890 #endif
 891 #ifdef MPI
 892 !C Determine the minimum energes for all parameter sets and temperatures
 893       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
 894        maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 895       do i=1,nParmSet
 896         do j=1,nT_h(i)
 897           potEmin_all(j,i)=potEmin_t_all(j,i)
 898         enddo
 899       enddo
 900 #endif
 901       potEmin_min=potEmin_all(1,1)
 902       do i=1,nParmSet
 903         do j=1,nT_h(i)
 904           if (potEmin_all(j,i).lt.potEmin_min) &
 905                   potEmin_min=potEmin_all(j,i)
 906         enddo
 907       enddo
 908 #ifdef DEBUG
 909       write (iout,*) "The potEmin array"
 910       do i=1,nParmSet
 911         write (iout,*) "Parameter set",i
 912         do j=1,nT_h(i)
 913           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
 914            PotEmin_all(j,i)
 915         enddo
 916       enddo
 917       write (iout,*) "potEmin_min",potEmin_min
 918 #endif
 919
 920 #ifdef MPI
 921       do t=0,tmax
 922         hfin_ent_p(t)=0.0d0
 923       enddo
 924 #else
 925       do t=0,tmax
 926         hfin_ent(t)=0.0d0
 927       enddo
 928 #endif
 929       write (iout,*) "--------------hist"
 930 #ifdef MPI
 931       do iparm=1,nParmSet
 932         do i=0,nGridT
 933           sumW_p(i,iparm)=0.0d0
 934           sumE_p(i,iparm)=0.0d0
 935           sumEbis_p(i,iparm)=0.0d0
 936           sumEsq_p(i,iparm)=0.0d0
 937           do j=1,nQ+2
 938             sumQ_p(j,i,iparm)=0.0d0
 939             sumQsq_p(j,i,iparm)=0.0d0
 940             sumEQ_p(j,i,iparm)=0.0d0
 941           enddo
 942         enddo
 943       enddo
 944       upindE_p=0
 945 #else
 946       do iparm=1,nParmSet
 947         do i=0,nGridT
 948           sumW(i,iparm)=0.0d0
 949           sumE(i,iparm)=0.0d0
 950           sumEbis(i,iparm)=0.0d0
 951           sumEsq(i,iparm)=0.0d0
 952           do j=1,nQ+2
 953             sumQ(j,i,iparm)=0.0d0
 954             sumQsq(j,i,iparm)=0.0d0
 955             sumEQ(j,i,iparm)=0.0d0
 956           enddo
 957         enddo
 958       enddo
 959       upindE=0
 960 #endif
 961 ! 8/26/05 entropy distribution
 962 !#ifdef MPI
 963 !      entmin_p=1.0d10
 964 !      entmax_p=-1.0d10
 965 !      do t=1,scount(me1)
 966 !!        ent=-dlog(entfac(t))
 967 !        ent=entfac(t)
 968 !        if (ent.lt.entmin_p) entmin_p=ent
 969 !        if (ent.gt.entmax_p) entmax_p=ent
 970 !      enddo
 971 !      write (iout,*) "entmin",entmin_p," entmax",entmax_p
 972 !!      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
 973 !      call flush(iout)
 974 !      call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
 975 !        WHAM_COMM,IERROR)
 976 !      call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
 977 !        WHAM_COMM,IERROR)
 978 !      write (iout,*) "entmin",entmin," entmax",entmax
 979 !      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
 980 !      ientmax=entmax-entmin
 981 !iientmax=entmax-entmin !el
 982 !write (iout,*) "ientmax",ientmax,entmax,entmin
 983 !write (iout,*) "iientmax",iientmax
 984 !      if (ientmax.gt.2000) ientmax=2000
 985 !      write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 986 !      call flush(iout)
 987 !      do t=1,scount(me1)
 988 !!        ient=-dlog(entfac(t))-entmin
 989 !        ient=entfac(t)-entmin
 990 !        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 991 !      enddo
 992 !      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
 993 !        MPI_SUM,WHAM_COMM,IERROR)
 994 !      if (me1.eq.Master) then
 995 !        write (iout,*) "Entropy histogram"
 996 !        do i=0,ientmax
 997 !          write(iout,'(f15.4,i10)') entmin+i,histent(i)
 998 !        enddo
 999 !      endif
1000 !#else
1001 !      entmin=1.0d10
1002 !      entmax=-1.0d10
1003 !      do t=1,ntot(islice)
1004 !        ent=entfac(t)
1005 !        if (ent.lt.entmin) entmin=ent
1006 !        if (ent.gt.entmax) entmax=ent
1007 !      enddo
1008 !      ientmax=-dlog(entmax)-entmin
1009 !      if (ientmax.gt.2000) ientmax=2000
1010 !      do t=1,ntot(islice)
1011 !        ient=entfac(t)-entmin
1012 !        if (ient.le.2000) histent(ient)=histent(ient)+1
1013 !      enddo
1014 !      write (iout,*) "Entropy histogram"
1015 !      do i=0,ientmax
1016 !        write(iout,'(2f15.4)') entmin+i,histent(i)
1017 !      enddo
1018 !#endif
1019
1020 #ifdef MPI
1021       write (iout,*) "me1",me1," scount",scount(me1) !d
1022
1023       do iparm=1,nParmSet
1024
1025 #ifdef MPI
1026         do ib=1,nT_h(iparm)
1027           do t=0,tmax
1028             hfin_p(t,ib)=0.0d0
1029           enddo
1030         enddo
1031         do i=1,maxindE
1032           histE_p(i)=0.0d0
1033         enddo
1034 #else
1035         do ib=1,nT_h(iparm)
1036           do t=0,tmax
1037             hfin(t,ib)=0.0d0
1038           enddo
1039         enddo
1040         do i=1,maxindE
1041           histE(i)=0.0d0
1042         enddo
1043 #endif
1044         do ib=1,nT_h(iparm)
1045           do i=0,MaxBinRms
1046             do j=0,MaxBinRgy
1047               hrmsrgy(j,i,ib)=0.0d0
1048 #ifdef MPI
1049               hrmsrgy_p(j,i,ib)=0.0d0
1050 #endif
1051             enddo
1052           enddo
1053         enddo
1054
1055         do t=1,scount(me1)
1056 #else
1057         do t=1,ntot(islice)
1058 #endif
1059           ind = ind_point(t)
1060 #ifdef MPI
1061           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1062 #else
1063           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1064 #endif
1065 !          write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
1066           call restore_parm(iparm)
1067 !          evdw=enetb(21,t,iparm)
1068           evdw=enetb(20,t,iparm)
1069           evdw_t=enetb(1,t,iparm)
1070 #ifdef SCP14
1071 !          evdw2_14=enetb(17,t,iparm)
1072           evdw2_14=enetb(18,t,iparm)
1073           evdw2=enetb(2,t,iparm)+evdw2_14
1074 #else
1075           evdw2=enetb(2,t,iparm)
1076           evdw2_14=0.0d0
1077 #endif
1078 #ifdef SPLITELE
1079           ees=enetb(3,t,iparm)
1080           evdw1=enetb(16,t,iparm)
1081 #else
1082           ees=enetb(3,t,iparm)
1083           evdw1=0.0d0
1084 #endif
1085           ecorr=enetb(4,t,iparm)
1086           ecorr5=enetb(5,t,iparm)
1087           ecorr6=enetb(6,t,iparm)
1088           eel_loc=enetb(7,t,iparm)
1089           eello_turn3=enetb(8,t,iparm)
1090           eello_turn4=enetb(9,t,iparm)
1091           eello_turn6=enetb(10,t,iparm)
1092           ebe=enetb(11,t,iparm)
1093           escloc=enetb(12,t,iparm)
1094           etors=enetb(13,t,iparm)
1095           etors_d=enetb(14,t,iparm)
1096           ehpb=enetb(15,t,iparm)
1097 !          estr=enetb(18,t,iparm)
1098           estr=enetb(17,t,iparm)
1099 !          esccor=enetb(19,t,iparm)
1100           esccor=enetb(21,t,iparm)
1101 !          edihcnstr=enetb(20,t,iparm)
1102           edihcnstr=enetb(19,t,iparm)
1103           edihcnstr=0.0d0
1104           ecationcation=enetb(41,t,iparm)
1105           ecation_prot=enetb(42,t,iparm)
1106             evdwpp =      enetb(26,t,iparm)
1107             eespp =      enetb(27,t,iparm)
1108             evdwpsb =      enetb(28,t,iparm)
1109             eelpsb =      enetb(29,t,iparm)
1110             evdwsb =      enetb(30,t,iparm)
1111             eelsb =      enetb(31,t,iparm)
1112             estr_nucl =      enetb(32,t,iparm)
1113             ebe_nucl =      enetb(33,t,iparm)
1114             esbloc =      enetb(34,t,iparm)
1115             etors_nucl =      enetb(35,t,iparm)
1116             etors_d_nucl =      enetb(36,t,iparm)
1117             ecorr_nucl =      enetb(37,t,iparm)
1118             ecorr3_nucl =      enetb(38,t,iparm)
1119             epeppho=   enetb(49,t,iparm)
1120             escpho=    enetb(48,t,iparm)
1121             epepbase=  enetb(47,t,iparm)
1122             escbase=   enetb(46,t,iparm)
1123             ecation_nucl=enetb(50,t,iparm)
1124             elipbond=enetb(52,t,iparm)
1125             elipang=enetb(53,t,iparm)
1126             eliplj=enetb(54,t,iparm)
1127             elipelec=enetb(55,t,iparm)
1128             ecat_prottran=enetb(56,t,iparm)
1129             ecation_protang=enetb(57,t,iparm)
1130             eliptran=enetb(22,t,iparm)
1131           do k=0,nGridT
1132             betaT=startGridT+k*delta_T
1133             temper=betaT
1134 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
1135 !d            fT=T0/betaT
1136 !d            ft=2*T0/(T0+betaT)
1137             if (rescale_modeW.eq.1) then
1138               quot=betaT/T0
1139               quotl=1.0d0
1140               kfacl=1.0d0
1141               do l=1,5
1142                 quotl1=quotl
1143                 quotl=quotl*quot
1144                 kfacl=kfacl*kfac
1145                 denom=kfacl-1.0d0+quotl
1146                 fT(l)=kfacl/denom
1147                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1148                 ftbis(l)=l*kfacl*quotl1* &
1149                  (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1150               enddo
1151 #if defined(FUNCTH)
1152               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1153                         320.0d0
1154               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1155              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1156                     /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1157 #elif defined(FUNCT)
1158               fT(6)=betaT/T0
1159               ftprim(6)=1.0d0/T0
1160               ftbis(6)=0.0d0
1161 #else
1162               fT(6)=1.0d0
1163               ftprim(6)=0.0d0
1164               ftbis(6)=0.0d0
1165 #endif
1166             else if (rescale_modeW.eq.2) then
1167               quot=betaT/T0
1168               quotl=1.0d0
1169               do l=1,5
1170                 quotl1=quotl
1171                 quotl=quotl*quot
1172                 eplus=dexp(quotl)
1173                 eminus=dexp(-quotl)
1174                 logfac=1.0d0/dlog(eplus+eminus)
1175                 tanhT=(eplus-eminus)/(eplus+eminus)
1176                 fT(l)=1.12692801104297249644d0*logfac
1177                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1178                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1179                 2*l*quotl1/T0*logfac* &
1180                 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1181                 +ftprim(l)*tanhT)
1182               enddo
1183 #if defined(FUNCTH)
1184               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1185                        320.0d0
1186               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1187              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1188                      /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1189 #elif defined(FUNCT)
1190               fT(6)=betaT/T0
1191               ftprim(6)=1.0d0/T0
1192               ftbis(6)=0.0d0
1193 #else
1194               fT(6)=1.0d0
1195               ftprim(6)=0.0d0
1196               ftbis(6)=0.0d0
1197 #endif
1198             else if (rescale_modeW.eq.0) then
1199               do l=1,5
1200                 fT(l)=1.0d0
1201                 ftprim(l)=0.0d0
1202               enddo
1203             else
1204               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
1205                 rescale_modeW
1206               call flush(iout)
1207               return 1
1208             endif
1209 !            write (iout,*) "ftprim",ftprim
1210 !            write (iout,*) "ftbis",ftbis
1211             betaT=1.0d0/(1.987D-3*betaT)
1212 #ifdef SPLITELE
1213             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1214             +wvdwpp*evdw1 &
1215             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1216             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1217             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1218             +ft(2)*wturn3*eello_turn3 &
1219             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1220             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1221             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1222             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1223             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1224             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1225             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1226             +wcorr3_nucl*ecorr3_nucl&
1227             +wscbase*escbase&
1228             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1229             + wcatnucl*ecation_nucl&
1230             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1231             +wliptran*eliptran
1232
1233
1234
1235             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
1236                  +ftprim(1)*wtor*etors+ &
1237                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1238                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1239                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1240                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1241                  ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
1242                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1243             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
1244                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1245                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1246                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1247                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1248                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1249                  +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
1250 #else
1251             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1252             +ft(1)*welec*(ees+evdw1) &
1253             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1254             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1255             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1256             +ft(2)*wturn3*eello_turn3 &
1257             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1258             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1259             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1260             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1261             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1262             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1263             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1264             +wcorr3_nucl*ecorr3_nucl&
1265             +wscbase*escbase&
1266             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1267             + wcatnucl*ecation_nucl&
1268             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1269             +wliptran*eliptran
1270
1271
1272
1273             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1274                 +ftprim(1)*wtor*etors+ &
1275                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1276                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1277                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1278                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1279                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1280                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1281             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1282                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1283                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1284                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1285                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1286                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1287                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1288
1289 #endif
1290             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1291 #ifdef DEBUG
1292             write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
1293               " etot",etot," entfac",entfac(t),&
1294               " weight",weight," ebis",ebis
1295 #endif
1296             etot=etot-temper*eprim
1297 #ifdef MPI
1298             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1299             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1300             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1301             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1302             do j=1,nQ+2
1303               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1304               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1305               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
1306                +etot*q(j,t)*weight
1307             enddo
1308 #else
1309             sumW(k,iparm)=sumW(k,iparm)+weight
1310             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1311             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1312             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1313             do j=1,nQ+2
1314               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1315               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1316               sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
1317                +etot*q(j,t)*weight
1318             enddo
1319 #endif
1320           enddo
1321           indE = aint(potE(t,iparm)-aint(potEmin))
1322           if (indE.ge.0 .and. indE.le.maxinde) then
1323             if (indE.gt.upindE_p) upindE_p=indE
1324             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1325           endif
1326 #ifdef MPI
1327           do ib=1,nT_h(iparm)
1328             potEmin=potEmin_all(ib,iparm)
1329             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1330             hfin_p(ind,ib)=hfin_p(ind,ib)+ &
1331              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1332              if (rmsrgymap) then
1333                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1334                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1335                hrmsrgy_p(indrgy,indrms,ib)= &
1336                  hrmsrgy_p(indrgy,indrms,ib)+expfac
1337              endif
1338           enddo
1339 #else
1340           do ib=1,nT_h(iparm)
1341             potEmin=potEmin_all(ib,iparm)
1342             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1343             hfin(ind,ib)=hfin(ind,ib)+ &
1344              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1345              if (rmsrgymap) then
1346                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1347                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1348                hrmsrgy(indrgy,indrms,ib)= &
1349                  hrmsrgy(indrgy,indrms,ib)+expfac
1350              endif
1351           enddo
1352 #endif
1353         enddo ! t
1354         do ib=1,nT_h(iparm)
1355           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
1356           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1357           if (rmsrgymap) then
1358           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
1359          (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1360              WHAM_COMM,IERROR)
1361           endif
1362         enddo
1363         call MPI_Reduce(upindE_p,upindE,1,&
1364           MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1365         call MPI_Reduce(histE_p(0),histE(0),maxindE,&
1366           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1367
1368         if (me1.eq.master) then
1369
1370         if (histout) then
1371
1372         write (iout,'(6x,$)')
1373         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
1374          ib=1,nT_h(iparm))
1375         write (iout,*)
1376
1377         write (iout,'(/a)') 'Final histograms'
1378         if (histfile) then
1379           if (nslice.eq.1) then
1380             if (separate_parset) then
1381               write(licz3,"(bz,i3.3)") myparm
1382               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1383             else
1384               histname=prefix(:ilen(prefix))//'.hist'
1385             endif
1386           else
1387             if (separate_parset) then
1388               write(licz3,"(bz,i3.3)") myparm
1389               histname=prefix(:ilen(prefix))//'_par'//licz3// &
1390                '_slice_'//licz2//'.hist'
1391             else
1392               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1393             endif
1394           endif
1395 #if defined(AIX) || defined(PGI)
1396           open (ihist,file=histname,position='append')
1397 #else
1398           open (ihist,file=histname,access='append')
1399 #endif
1400         endif
1401
1402         do t=0,tmax
1403           liczbaW=t
1404           sumH=0.0d0
1405           do ib=1,nT_h(iparm)
1406             sumH=sumH+hfin(t,ib)
1407           enddo
1408           if (sumH.gt.0.0d0) then
1409             do j=1,nQ
1410               jj = mod(liczbaW,nbin1)
1411               liczbaW=liczbaW/nbin1
1412               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1413               if (histfile) &
1414                  write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1415             enddo
1416             do ib=1,nT_h(iparm)
1417               write (iout,'(e20.10,$)') hfin(t,ib)
1418               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1419             enddo
1420             write (iout,'(i5)') iparm
1421             if (histfile) write (ihist,'(i5)') iparm
1422           endif
1423         enddo
1424
1425         endif
1426
1427         if (entfile) then
1428           if (nslice.eq.1) then
1429             if (separate_parset) then
1430               write(licz3,"(bz,i3.3)") myparm
1431               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1432             else
1433               histname=prefix(:ilen(prefix))//'.ent'
1434             endif
1435           else
1436             if (separate_parset) then
1437               write(licz3,"(bz,i3.3)") myparm
1438               histname=prefix(:ilen(prefix))//'par_'//licz3// &
1439                  '_slice_'//licz2//'.ent'
1440             else
1441               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1442             endif
1443           endif
1444 #if defined(AIX) || defined(PGI)
1445           open (ihist,file=histname,position='append')
1446 #else
1447           open (ihist,file=histname,access='append')
1448 #endif
1449           write (ihist,'(a)') "# Microcanonical entropy"
1450           do i=0,upindE
1451             write (ihist,'(f8.0,$)') dint(potEmin)+i
1452             if (histE(i).gt.0.0e0) then
1453               write (ihist,'(f15.5,$)') dlog(histE(i))
1454             else
1455               write (ihist,'(f15.5,$)') 0.0d0
1456             endif
1457           enddo
1458           write (ihist,*)
1459           close(ihist)
1460         endif
1461         write (iout,*) "Microcanonical entropy"
1462         do i=0,upindE
1463           write (iout,'(f8.0,$)') dint(potEmin)+i
1464           if (histE(i).gt.0.0e0) then
1465             write (iout,'(f15.5,$)') dlog(histE(i))
1466           else
1467             write (iout,'(f15.5,$)') 0.0d0
1468           endif
1469           write (iout,*)
1470         enddo
1471         if (rmsrgymap) then
1472           if (nslice.eq.1) then
1473             if (separate_parset) then
1474               write(licz3,"(bz,i3.3)") myparm
1475               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1476             else
1477               histname=prefix(:ilen(prefix))//'.rmsrgy'
1478             endif
1479           else
1480             if (separate_parset) then
1481               write(licz3,"(bz,i3.3)") myparm
1482               histname=prefix(:ilen(prefix))//'_par'//licz3// &
1483                '_slice_'//licz2//'.rmsrgy'
1484             else
1485              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1486             endif
1487           endif
1488 #if defined(AIX) || defined(PGI)
1489           open (ihist,file=histname,position='append')
1490 #else
1491           open (ihist,file=histname,access='append')
1492 #endif
1493           do i=0,nbin_rms
1494             do j=0,nbin_rgy
1495               write(ihist,'(2f8.2,$)') &
1496                 rgymin+deltrgy*j,rmsmin+deltrms*i
1497               do ib=1,nT_h(iparm)
1498                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1499                   write(ihist,'(e14.5,$)') &
1500                     -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
1501                     +potEmin
1502                 else
1503                   write(ihist,'(e14.5,$)') 1.0d6
1504                 endif
1505               enddo
1506               write (ihist,'(i2)') iparm
1507             enddo
1508           enddo
1509           close(ihist)
1510         endif
1511         endif
1512       enddo ! iparm
1513 #ifdef MPI
1514       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
1515          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1516       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
1517          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1518       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
1519          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1520       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
1521          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1522       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
1523          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1524       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
1525          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1526          WHAM_COMM,IERROR)
1527       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
1528          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1529          WHAM_COMM,IERROR)
1530       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
1531          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
1532          WHAM_COMM,IERROR)
1533       if (me.eq.master) then
1534 #endif
1535       write (iout,'(/a)') 'Thermal characteristics of folding'
1536       if (nslice.eq.1) then
1537         nazwa=prefix
1538       else
1539         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1540       endif
1541       iln=ilen(nazwa)
1542       if (nparmset.eq.1 .and. .not.separate_parset) then
1543         nazwa=nazwa(:iln)//".thermal"
1544       else if (nparmset.eq.1 .and. separate_parset) then
1545         write(licz3,"(bz,i3.3)") myparm
1546         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1547       endif
1548       do iparm=1,nParmSet
1549       if (nparmset.gt.1) then
1550         write(licz3,"(bz,i3.3)") iparm
1551         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1552       endif
1553       open(34,file=nazwa)
1554       if (separate_parset) then
1555         write (iout,'(a,i3)') "Parameter set",myparm
1556       else
1557         write (iout,'(a,i3)') "Parameter set",iparm
1558       endif
1559       do i=0,NGridT
1560         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1561         if (betaT.ge.beta_h(1,iparm)) then
1562           potEmin=potEmin_all(1,iparm)
1563         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1564           potEmin=potEmin_all(nT_h(iparm),iparm)
1565         else
1566           do l=1,nT_h(iparm)-1
1567             if (betaT.le.beta_h(l,iparm) .and. &
1568                betaT.gt.beta_h(l+1,iparm)) then
1569               potEmin=potEmin_all(l,iparm)
1570               exit
1571             endif
1572           enddo
1573         endif
1574
1575         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1576         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
1577           sumW(i,iparm)
1578         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
1579           -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1580         do j=1,nQ+2
1581           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1582           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
1583            -sumQ(j,i,iparm)**2
1584           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
1585            -sumQ(j,i,iparm)*sumE(i,iparm)
1586         enddo
1587         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
1588          (startGridT+i*delta_T))+potEmin
1589         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1590          sumW(i,iparm),sumE(i,iparm)
1591         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1592         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1593          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1594         write (iout,*)
1595         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
1596          sumW(i,iparm),sumE(i,iparm)
1597         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1598         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
1599          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1600         write (34,*)
1601         call flush(34)
1602       enddo
1603       close(34)
1604       enddo
1605       if (histout) then
1606       do t=0,tmax
1607         if (hfin_ent(t).gt.0.0d0) then
1608           liczbaW=t
1609           jj = mod(liczbaW,nbin1)
1610           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
1611            hfin_ent(t)
1612           if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
1613             dmin+(jj+0.5d0)*delta,&
1614            hfin_ent(t)
1615         endif
1616       enddo
1617       if (histfile) close(ihist)
1618       endif
1619
1620 #ifdef ZSCORE
1621 ! Write data for zscore
1622       if (nslice.eq.1) then
1623         zscname=prefix(:ilen(prefix))//".zsc"
1624       else
1625         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1626       endif
1627 #if defined(AIX) || defined(PGI)
1628       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1629 #else
1630       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1631 #endif
1632       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1633       do iparm=1,nParmSet
1634         write (izsc,'("NT=",i1)') nT_h(iparm)
1635       do ib=1,nT_h(iparm)
1636         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
1637           1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1638         jj = min0(nR(ib,iparm),7)
1639         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1640         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1641         write (izsc,'("&")')
1642         if (nR(ib,iparm).gt.7) then
1643           do ii=8,nR(ib,iparm),9
1644             jj = min0(nR(ib,iparm),ii+8)
1645             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1646             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1647             write (izsc,'("&")')
1648           enddo
1649         endif
1650         write (izsc,'("FI=",$)')
1651         jj=min0(nR(ib,iparm),7)
1652         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1653         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1654         write (izsc,'("&")')
1655         if (nR(ib,iparm).gt.7) then
1656           do ii=8,nR(ib,iparm),9
1657             jj = min0(nR(ib,iparm),ii+8)
1658             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1659             if (jj.eq.nR(ib,iparm)) then
1660               write (izsc,*)
1661             else
1662               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1663               write (izsc,'(t80,"&")')
1664             endif
1665           enddo
1666         endif
1667         do i=1,nR(ib,iparm)
1668           write (izsc,'("KH=",$)')
1669           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1670           write (izsc,'(" Q0=",$)')
1671           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1672           write (izsc,*)
1673         enddo
1674       enddo
1675       enddo
1676       close(izsc)
1677 #endif
1678 #ifdef MPI
1679       endif
1680 #endif
1681       return
1682       end subroutine WHAMCALC
1683 !-----------------------------------------------------------------------------
1684       end module wham_calc
1685