source/wham/src-M/wham_calc1.F

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