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