source/wham/src/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,16f12.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+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+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 #ifdef DEBUG
 360             write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
 361 #endif
 362             if (homol_nset.gt.1) then
 363
 364             do kk=1,nR(ib,iparm)
 365               Econstr=waga_homology(kk)*ehomology_constr
 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,Econstr,v(i,kk,ib,iparm)
 371 #endif
 372             enddo ! kk
 373
 374             else
 375
 376             etot=etot+ehomology_constr
 377             do kk=1,nR(ib,iparm)
 378               Econstr=0.0d0
 379               do j=1,nQ
 380                 dd = q(j,i)
 381                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 382      &           *(dd-q0(j,kk,ib,iparm))**2
 383               enddo
 384               v(i,kk,ib,iparm)=
 385      &          -beta_h(ib,iparm)*(etot+Econstr)
 386 #ifdef DEBUG
 387               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 388      &         etot,v(i,kk,ib,iparm)
 389 #endif
 390             enddo ! kk
 391
 392             endif
 393           enddo   ! ib
 394         enddo     ! iparm
 395       enddo       ! i
 396 ! Simple iteration to calculate free energies corresponding to all simulation
 397 ! runs.
 398       do iter=1,maxit
 399
 400 ! Compute new free-energy values corresponding to the righ-hand side of the
 401 ! equation and their derivatives.
 402         write (iout,*) "------------------------fi"
 403         entfac_min=1.0d10
 404 #ifdef MPI
 405         do t=1,scount(me1)
 406 #else
 407         do t=1,ntot(islice)
 408 #endif
 409           vmax=-1.0d+20
 410           do i=1,nParmSet
 411             do k=1,nT_h(i)
 412               do l=1,nR(k,i)
 413                 vf=v(t,l,k,i)+f(l,k,i)
 414                 if (vf.gt.vmax) vmax=vf
 415               enddo
 416             enddo
 417           enddo
 418           denom=0.0d0
 419           do i=1,nParmSet
 420             do k=1,nT_h(i)
 421               do l=1,nR(k,i)
 422                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 423                 if (aux.gt.-200.0d0)
 424      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 425               enddo
 426             enddo
 427           enddo
 428           entfac(t)=-dlog(denom)-vmax
 429           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 430 #ifdef DEBUG
 431           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 432 #endif
 433         enddo
 434 c#ifdef MPI
 435 c        write (iout,*) "entfac_min before AllReduce",entfac_min
 436 c        call MPI_AllReduce(entfac_min,entfac_min_t,1,
 437 c     &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 438 c        entfac_min=entfac_min_t
 439 c        write (iout,*) "entfac_min after AllReduce",entfac_min
 440 c#endif
 441 c#ifdef MPI
 442 c        do t=1,scount(me)
 443 c          entfac(t)=entfac(t)-entfac_min
 444 c        enddo
 445 c#else
 446 c        do t=1,ntot(islice)
 447 c          entfac(t)=entfac(t)-entfac_min
 448 c        enddo
 449 c#endif
 450         do iparm=1,nParmSet
 451           do iib=1,nT_h(iparm)
 452             do ii=1,nR(iib,iparm)
 453 #ifdef MPI
 454               fi_p_min(ii,iib,iparm)=-1.0d10
 455               do t=1,scount(me)
 456                 aux=v(t,ii,iib,iparm)+entfac(t)
 457                 if (aux.gt.fi_p_min(ii,iib,iparm))
 458      &                                   fi_p_min(ii,iib,iparm)=aux
 459               enddo
 460 #else
 461               do t=1,ntot(islice)
 462                 aux=v(t,ii,iib,iparm)+entfac(t)
 463                 if (aux.gt.fi_min(ii,iib,iparm))
 464      &                                     fi_min(ii,iib,iparm)=aux
 465               enddo
 466 #endif
 467             enddo ! ii
 468           enddo ! iib
 469         enddo ! iparm
 470 #ifdef MPI
 471 #ifdef DEBUG
 472         write (iout,*) "fi_min before AllReduce"
 473         do i=1,nParmSet
 474           do j=1,nT_h(i)
 475             write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
 476           enddo
 477         enddo
 478 #endif
 479         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
 480      &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 481 #ifdef DEBUG
 482         write (iout,*) "fi_min after AllReduce"
 483         do i=1,nParmSet
 484           do j=1,nT_h(i)
 485             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 486           enddo
 487         enddo
 488 #endif
 489 #endif
 490         do iparm=1,nParmSet
 491           do iib=1,nT_h(iparm)
 492             do ii=1,nR(iib,iparm)
 493 #ifdef MPI
 494               fi_p(ii,iib,iparm)=0.0d0
 495               do t=1,scount(me)
 496                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 497      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 498 #ifdef DEBUG
 499               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
 500      &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
 501      &         fi_p(ii,iib,iparm)
 502 #endif
 503               enddo
 504 #else
 505               fi(ii,iib,iparm)=0.0d0
 506               do t=1,ntot(islice)
 507                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 508      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 509               enddo
 510 #endif
 511             enddo ! ii
 512           enddo ! iib
 513         enddo ! iparm
 514
 515 #ifdef MPI
 516 #ifdef DEBUG
 517         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 518         do iparm=1,nParmSet
 519           do ib=1,nT_h(iparm)
 520             write (iout,*) "iparm",iparm," ib",ib
 521             write (iout,*) "beta=",beta_h(ib,iparm)
 522             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 523           enddo
 524         enddo
 525 #endif
 526 #ifdef DEBUG
 527         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 528      &   maxR*MaxT_h*nParmSet
 529         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 530      &   " WHAM_COMM",WHAM_COMM
 531 #endif
 532         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 533      &   MPI_DOUBLE_PRECISION,
 534      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 535 #ifdef DEBUG
 536         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 537         do iparm=1,nParmSet
 538           write (iout,*) "iparm",iparm
 539           do ib=1,nT_h(iparm)
 540             write (iout,*) "beta=",beta_h(ib,iparm)
 541             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 542           enddo
 543         enddo
 544 #endif
 545         if (me1.eq.Master) then
 546 #endif
 547         avefi=0.0d0
 548         do iparm=1,nParmSet
 549           do ib=1,nT_h(iparm)
 550             do i=1,nR(ib,iparm)
 551               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 552               avefi=avefi+fi(i,ib,iparm)
 553             enddo
 554           enddo
 555         enddo
 556         avefi=avefi/nFi
 557         do iparm=1,nParmSet
 558           write (iout,*) "Parameter set",iparm
 559           do ib =1,nT_h(iparm)
 560             write (iout,*) "beta=",beta_h(ib,iparm)
 561             do i=1,nR(ib,iparm)
 562               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 563             enddo
 564             write (iout,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 565             write (iout,'(6f15.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 566           enddo
 567         enddo
 568
 569 ! Compute the norm of free-energy increments.
 570         finorm=0.0d0
 571         do iparm=1,nParmSet
 572           do ib=1,nT_h(iparm)
 573             do i=1,nR(ib,iparm)
 574               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 575               f(i,ib,iparm)=fi(i,ib,iparm)
 576             enddo
 577           enddo
 578         enddo
 579
 580         write (iout,*) 'Iteration',iter,' finorm',finorm
 581
 582 #ifdef MPI
 583         endif
 584         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 585      &   MPI_DOUBLE_PRECISION,Master,
 586      &   WHAM_COMM,IERROR)
 587         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 588      &   WHAM_COMM,IERROR)
 589 #endif
 590 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 591         if (finorm.lt.fimin) then
 592           write (iout,*) 'Iteration converged'
 593           goto 20
 594         endif
 595
 596       enddo ! iter
 597
 598    20 continue
 599 ! Now, put together the histograms from all simulations, in order to get the
 600 ! unbiased total histogram.
 601
 602 C Determine the minimum free energies
 603 #ifdef MPI
 604       do i=1,scount(me1)
 605 #else
 606       do i=1,ntot(islice)
 607 #endif
 608 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 609         do iparm=1,nParmSet
 610 #ifdef DEBUG
 611           write (iout,'(2i5,21f8.2)') i,iparm,
 612      &     (enetb(k,i,iparm),k=1,22)
 613 #endif
 614           call restore_parm(iparm)
 615 #ifdef DEBUG
 616           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 617      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 618      &      wtor_d,wsccor,wbond
 619 #endif
 620           do ib=1,nT_h(iparm)
 621             if (rescale_mode.eq.1) then
 622               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 623               quotl=1.0d0
 624               kfacl=1.0d0
 625               do l=1,5
 626                 quotl1=quotl
 627                 quotl=quotl*quot
 628                 kfacl=kfacl*kfac
 629                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 630               enddo
 631 #if defined(FUNCTH)
 632               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 633               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 634 #elif defined(FUNCT)
 635               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 636 #else
 637               ft(6)=1.0d0
 638 #endif
 639             else if (rescale_mode.eq.2) then
 640               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 641               quotl=1.0d0
 642               do l=1,5
 643                 quotl=quotl*quot
 644                 fT(l)=1.12692801104297249644d0/
 645      &             dlog(dexp(quotl)+dexp(-quotl))
 646               enddo
 647 #if defined(FUNCTH)
 648               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 649               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 650 #elif defined(FUNCT)
 651               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 652 #else
 653               ft(6)=1.0d0
 654 #endif
 655 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 656             else if (rescale_mode.eq.0) then
 657               do l=1,6
 658                 fT(l)=1.0d0
 659               enddo
 660             else
 661               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 662      &          rescale_mode
 663               call flush(iout)
 664               return1
 665             endif
 666             evdw=enetb(1,i,iparm)
 667             evdw_t=enetb(21,i,iparm)
 668 #ifdef SCP14
 669             evdw2_14=enetb(17,i,iparm)
 670             evdw2=enetb(2,i,iparm)+evdw2_14
 671 #else
 672             evdw2=enetb(2,i,iparm)
 673             evdw2_14=0.0d0
 674 #endif
 675 #ifdef SPLITELE
 676             ees=enetb(3,i,iparm)
 677             evdw1=enetb(16,i,iparm)
 678 #else
 679             ees=enetb(3,i,iparm)
 680             evdw1=0.0d0
 681 #endif
 682             ecorr=enetb(4,i,iparm)
 683             ecorr5=enetb(5,i,iparm)
 684             ecorr6=enetb(6,i,iparm)
 685             eel_loc=enetb(7,i,iparm)
 686             eello_turn3=enetb(8,i,iparm)
 687             eello_turn4=enetb(9,i,iparm)
 688             eturn6=enetb(10,i,iparm)
 689             ebe=enetb(11,i,iparm)
 690             escloc=enetb(12,i,iparm)
 691             etors=enetb(13,i,iparm)
 692             etors_d=enetb(14,i,iparm)
 693             ehpb=enetb(15,i,iparm)
 694             estr=enetb(18,i,iparm)
 695             esccor=enetb(19,i,iparm)
 696             edihcnstr=enetb(20,i,iparm)
 697             ehomology_constr=enetb(22,i,iparm)
 698             if (homol_nset.gt.1)
 699      &       ehomology_constr=waga_homology(ihset)*ehomology_constr
 700             edfadis=enetb(23,i,iparm)
 701             edfator=enetb(24,i,iparm)
 702             edfanei=enetb(25,i,iparm)
 703             edfabet=enetb(26,i,iparm)
 704 #ifdef DEBUG
 705             write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
 706      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 707      &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
 708      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
 709      &       wdfa_nei*edfanei+wdfa_beta*edfabet
 710 #endif
 711
 712 #ifdef SPLITELE
 713             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 714      &      +wvdwpp*evdw1
 715      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 716      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 717      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 718      &      +ft(2)*wturn3*eello_turn3
 719      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 720      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 721      &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
 722      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 723 #else
 724             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 725      &      +ft(1)*welec*(ees+evdw1)
 726      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 727      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 728      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 729      &      +ft(2)*wturn3*eello_turn3
 730      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 731      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 732      &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
 733      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 734 #endif
 735 #ifdef DEBUG
 736             write (iout,*) "i",i," ib",ib,
 737      &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
 738      &      " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm)
 739 #endif
 740             etot=etot-entfac(i)/beta_h(ib,iparm)
 741             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 742 #ifdef DEBUG
 743             write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
 744 #endif
 745           enddo   ! ib
 746         enddo     ! iparm
 747       enddo       ! i
 748 #ifdef DEBUG
 749       write (iout,*) "The potEmin array before reduction"
 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 #ifdef MPI
 759 C Determine the minimum energes for all parameter sets and temperatures
 760       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
 761      &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 762       do i=1,nParmSet
 763         do j=1,nT_h(i)
 764           potEmin_all(j,i)=potEmin_t_all(j,i)
 765         enddo
 766       enddo
 767 #endif
 768       potEmin_min=potEmin_all(1,1)
 769       do i=1,nParmSet
 770         do j=1,nT_h(i)
 771           if (potEmin_all(j,i).lt.potEmin_min)
 772      &             potEmin_min=potEmin_all(j,i)
 773         enddo
 774       enddo
 775 #ifdef DEBUG
 776       write (iout,*) "The potEmin array"
 777       do i=1,nParmSet
 778         write (iout,*) "Parameter set",i
 779         do j=1,nT_h(i)
 780           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
 781      &      PotEmin_all(j,i)
 782         enddo
 783       enddo
 784       write (iout,*) "potEmin_min",potEmin_min
 785 #endif
 786
 787 #ifdef MPI
 788       do t=0,tmax
 789         hfin_ent_p(t)=0.0d0
 790       enddo
 791 #else
 792       do t=0,tmax
 793         hfin_ent(t)=0.0d0
 794       enddo
 795 #endif
 796       write (iout,*) "--------------hist"
 797 #ifdef MPI
 798       do iparm=1,nParmSet
 799         do i=0,nGridT
 800           sumW_p(i,iparm)=0.0d0
 801           sumE_p(i,iparm)=0.0d0
 802           sumEbis_p(i,iparm)=0.0d0
 803           sumEsq_p(i,iparm)=0.0d0
 804           do j=1,nQ+2
 805             sumQ_p(j,i,iparm)=0.0d0
 806             sumQsq_p(j,i,iparm)=0.0d0
 807             sumEQ_p(j,i,iparm)=0.0d0
 808           enddo
 809         enddo
 810       enddo
 811       upindE_p=0
 812 #else
 813       do iparm=1,nParmSet
 814         do i=0,nGridT
 815           sumW(i,iparm)=0.0d0
 816           sumE(i,iparm)=0.0d0
 817           sumEbis(i,iparm)=0.0d0
 818           sumEsq(i,iparm)=0.0d0
 819           do j=1,nQ+2
 820             sumQ(j,i,iparm)=0.0d0
 821             sumQsq(j,i,iparm)=0.0d0
 822             sumEQ(j,i,iparm)=0.0d0
 823           enddo
 824         enddo
 825       enddo
 826       upindE=0
 827 #endif
 828 c 8/26/05 entropy distribution
 829 #ifdef MPI
 830       entmin_p=1.0d10
 831       entmax_p=-1.0d10
 832       do t=1,scount(me1)
 833 c        ent=-dlog(entfac(t))
 834         ent=entfac(t)
 835         if (ent.lt.entmin_p) entmin_p=ent
 836         if (ent.gt.entmax_p) entmax_p=ent
 837       enddo
 838       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 839       call flush(iout)
 840       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 841      &  WHAM_COMM,IERROR)
 842       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 843      &  WHAM_COMM,IERROR)
 844       ientmax=entmax-entmin
 845       if (ientmax.gt.2000) ientmax=2000
 846       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 847       call flush(iout)
 848       do t=1,scount(me1)
 849 c        ient=-dlog(entfac(t))-entmin
 850         ient=entfac(t)-entmin
 851         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 852       enddo
 853       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 854      &  MPI_SUM,WHAM_COMM,IERROR)
 855       if (me1.eq.Master) then
 856         write (iout,*) "Entropy histogram"
 857         do i=0,ientmax
 858           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 859         enddo
 860       endif
 861 #else
 862       entmin=1.0d10
 863       entmax=-1.0d10
 864       do t=1,ntot(islice)
 865         ent=entfac(t)
 866         if (ent.lt.entmin) entmin=ent
 867         if (ent.gt.entmax) entmax=ent
 868       enddo
 869       ientmax=-dlog(entmax)-entmin
 870       if (ientmax.gt.2000) ientmax=2000
 871       do t=1,ntot(islice)
 872         ient=entfac(t)-entmin
 873         if (ient.le.2000) histent(ient)=histent(ient)+1
 874       enddo
 875       write (iout,*) "Entropy histogram"
 876       do i=0,ientmax
 877         write(iout,'(2f15.4)') entmin+i,histent(i)
 878       enddo
 879 #endif
 880
 881 #ifdef MPI
 882 c      write (iout,*) "me1",me1," scount",scount(me1)
 883
 884       do iparm=1,nParmSet
 885
 886 #ifdef MPI
 887         do ib=1,nT_h(iparm)
 888           do t=0,tmax
 889             hfin_p(t,ib)=0.0d0
 890           enddo
 891         enddo
 892         do i=1,maxindE
 893           histE_p(i)=0.0d0
 894         enddo
 895 #else
 896         do ib=1,nT_h(iparm)
 897           do t=0,tmax
 898             hfin(t,ib)=0.0d0
 899           enddo
 900         enddo
 901         do i=1,maxindE
 902           histE(i)=0.0d0
 903         enddo
 904 #endif
 905         do ib=1,nT_h(iparm)
 906           do i=0,MaxBinRms
 907             do j=0,MaxBinRgy
 908               hrmsrgy(j,i,ib)=0.0d0
 909 #ifdef MPI
 910               hrmsrgy_p(j,i,ib)=0.0d0
 911 #endif
 912             enddo
 913           enddo
 914         enddo
 915
 916         do t=1,scount(me1)
 917 #else
 918         do t=1,ntot(islice)
 919 #endif
 920           ind = ind_point(t)
 921 #ifdef MPI
 922           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 923 #else
 924           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 925 #endif
 926           call restore_parm(iparm)
 927           evdw=enetb(21,t,iparm)
 928           evdw_t=enetb(1,t,iparm)
 929 #ifdef SCP14
 930           evdw2_14=enetb(17,t,iparm)
 931           evdw2=enetb(2,t,iparm)+evdw2_14
 932 #else
 933           evdw2=enetb(2,t,iparm)
 934           evdw2_14=0.0d0
 935 #endif
 936 #ifdef SPLITELE
 937           ees=enetb(3,t,iparm)
 938           evdw1=enetb(16,t,iparm)
 939 #else
 940           ees=enetb(3,t,iparm)
 941           evdw1=0.0d0
 942 #endif
 943           ecorr=enetb(4,t,iparm)
 944           ecorr5=enetb(5,t,iparm)
 945           ecorr6=enetb(6,t,iparm)
 946           eel_loc=enetb(7,t,iparm)
 947           eello_turn3=enetb(8,t,iparm)
 948           eello_turn4=enetb(9,t,iparm)
 949           eturn6=enetb(10,t,iparm)
 950           ebe=enetb(11,t,iparm)
 951           escloc=enetb(12,t,iparm)
 952           etors=enetb(13,t,iparm)
 953           etors_d=enetb(14,t,iparm)
 954           ehpb=enetb(15,t,iparm)
 955           estr=enetb(18,t,iparm)
 956           esccor=enetb(19,t,iparm)
 957           edihcnstr=enetb(20,t,iparm)
 958           ehomology_constr=enetb(22,t,iparm)
 959           if (homol_nset.gt.1)
 960      &       ehomology_constr=waga_homology(ihset)*ehomology_constr
 961           edfadis=enetb(23,t,iparm)
 962           edfator=enetb(24,t,iparm)
 963           edfanei=enetb(25,t,iparm)
 964           edfabet=enetb(26,t,iparm)
 965 c          do k=0,nGridT
 966           do k=0,nGridT
 967             betaT=startGridT+k*delta_T
 968             temper=betaT
 969 c            fT=T0/betaT
 970 c            ft=2*T0/(T0+betaT)
 971             if (rescale_mode.eq.1) then
 972               quot=betaT/T0
 973               quotl=1.0d0
 974               kfacl=1.0d0
 975               do l=1,5
 976                 quotl1=quotl
 977                 quotl=quotl*quot
 978                 kfacl=kfacl*kfac
 979                 denom=kfacl-1.0d0+quotl
 980                 fT(l)=kfacl/denom
 981                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 982                 ftbis(l)=l*kfacl*quotl1*
 983      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 984               enddo
 985 #if defined(FUNCTH)
 986               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 987      &                  320.0d0
 988               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 989              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 990      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 991 #elif defined(FUNCT)
 992               fT(6)=betaT/T0
 993               ftprim(6)=1.0d0/T0
 994               ftbis(6)=0.0d0
 995 #else
 996               fT(6)=1.0d0
 997               ftprim(6)=0.0d0
 998               ftbis(6)=0.0d0
 999 #endif
1000             else if (rescale_mode.eq.2) then
1001               quot=betaT/T0
1002               quotl=1.0d0
1003               do l=1,5
1004                 quotl1=quotl
1005                 quotl=quotl*quot
1006                 eplus=dexp(quotl)
1007                 eminus=dexp(-quotl)
1008                 logfac=1.0d0/dlog(eplus+eminus)
1009                 tanhT=(eplus-eminus)/(eplus+eminus)
1010                 fT(l)=1.12692801104297249644d0*logfac
1011                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1012                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1013      &          2*l*quotl1/T0*logfac*
1014      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1015      &          +ftprim(l)*tanhT)
1016               enddo
1017 #if defined(FUNCTH)
1018               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1019      &                 320.0d0
1020               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1021              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1022      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1023 #elif defined(FUNCT)
1024               fT(6)=betaT/T0
1025               ftprim(6)=1.0d0/T0
1026               ftbis(6)=0.0d0
1027 #else
1028               fT(6)=1.0d0
1029               ftprim(6)=0.0d0
1030               ftbis(6)=0.0d0
1031 #endif
1032             else if (rescale_mode.eq.0) then
1033               do l=1,5
1034                 fT(l)=1.0d0
1035                 ftprim(l)=0.0d0
1036               enddo
1037             else
1038               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1039      &          rescale_mode
1040               call flush(iout)
1041               return1
1042             endif
1043 c            write (iout,*) "ftprim",ftprim
1044 c            write (iout,*) "ftbis",ftbis
1045             betaT=1.0d0/(1.987D-3*betaT)
1046             if (betaT.ge.beta_h(1,iparm)) then
1047               potEmin=potEmin_all(1,iparm)+
1048      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
1049      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
1050      &         (1.0/betaT-1.0/beta_h(1,iparm))
1051 #ifdef DEBUG
1052               write(iout,*) "first",temper,potEmin
1053 #endif
1054             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1055               potEmin=potEmin_all(nT_h(iparm),iparm)+
1056      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1057      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1058      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1059 #ifdef DEBUG
1060               write (iout,*) "last",temper,potEmin
1061 #endif
1062             else
1063               do l=1,nT_h(iparm)-1
1064                 if (betaT.le.beta_h(l,iparm) .and.
1065      &              betaT.gt.beta_h(l+1,iparm)) then
1066                   potEmin=potEmin_all(l,iparm)
1067 #ifdef DEBUG
1068                   write (iout,*) "l",l,
1069      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1070      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1071 #endif
1072                   exit
1073                 endif
1074               enddo
1075             endif
1076 #ifdef DEBUG
1077             write (iout,*) "k",k," potEmin",potEmin
1078 #endif
1079 #ifdef SPLITELE
1080             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1081      &      +wvdwpp*evdw1
1082      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1083      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1084      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1085      &      +ft(2)*wturn3*eello_turn3
1086      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1087      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1088      &      +wbond*estr+ehomology_constr
1089             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1090      &            +ftprim(1)*wtor*etors+
1091      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1092      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1093      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1094      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1095      &            ftprim(1)*wsccor*esccor
1096             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1097      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1098      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1099      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1100      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1101      &            ftbis(1)*wsccor*esccor
1102 #else
1103             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1104      &      +ft(1)*welec*(ees+evdw1)
1105      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1106      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1107      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1108      &      +ft(2)*wturn3*eello_turn3
1109      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1110      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1111      &      +wbond*estr+ehomology_constr
1112             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1113      &           +ftprim(1)*wtor*etors+
1114      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1115      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1116      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1117      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1118      &            ftprim(1)*wsccor*esccor
1119             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1120      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1121      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1122      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1123      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1124      &            ftprim(1)*wsccor*esccor
1125 #endif
1126             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1127 #ifdef DEBUG
1128             write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
1129      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
1130      &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
1131      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
1132      &       wdfa_nei*edfanei+wdfa_beta*edfabet
1133             write (iout,*) "iparm",iparm," t",t," temper",temper,
1134      &        " etot",etot," entfac",entfac(t),
1135      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1136      &        " boltz",-betaT*(etot-potEmin)+entfac(t),
1137      &        " weight",weight," ebis",ebis
1138 #endif
1139             etot=etot-temper*eprim
1140 #ifdef MPI
1141             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1142             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1143             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1144             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1145             do j=1,nQ+2
1146               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1147               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1148               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1149      &         +etot*q(j,t)*weight
1150             enddo
1151 #else
1152             sumW(k,iparm)=sumW(k,iparm)+weight
1153             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1154             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1155             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1156             do j=1,nQ+2
1157               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1158               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1159               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1160      &         +etot*q(j,t)*weight
1161             enddo
1162 #endif
1163           enddo
1164           indE = aint(potE(t,iparm)-aint(potEmin))
1165           if (indE.ge.0 .and. indE.le.maxinde) then
1166             if (indE.gt.upindE_p) upindE_p=indE
1167             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1168           endif
1169 #ifdef MPI
1170           do ib=1,nT_h(iparm)
1171             potEmin=potEmin_all(ib,iparm)
1172             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1173             hfin_p(ind,ib)=hfin_p(ind,ib)+
1174      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1175              if (rmsrgymap) then
1176                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1177                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1178                hrmsrgy_p(indrgy,indrms,ib)=
1179      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1180              endif
1181           enddo
1182 #else
1183           do ib=1,nT_h(iparm)
1184             potEmin=potEmin_all(ib,iparm)
1185             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1186             hfin(ind,ib)=hfin(ind,ib)+
1187      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1188              if (rmsrgymap) then
1189                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1190                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1191                hrmsrgy(indrgy,indrms,ib)=
1192      &           hrmsrgy(indrgy,indrms,ib)+expfac
1193              endif
1194           enddo
1195 #endif
1196         enddo ! t
1197         do ib=1,nT_h(iparm)
1198           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1199      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1200           if (rmsrgymap) then
1201           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1202      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1203      &       WHAM_COMM,IERROR)
1204           endif
1205         enddo
1206         call MPI_Reduce(upindE_p,upindE,1,
1207      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1208         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1209      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1210
1211         if (me1.eq.master) then
1212
1213         if (histout) then
1214
1215         write (iout,'(6x,$)')
1216         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1217      &   ib=1,nT_h(iparm))
1218         write (iout,*)
1219
1220         write (iout,'(/a)') 'Final histograms'
1221         if (histfile) then
1222           if (nslice.eq.1) then
1223             if (separate_parset) then
1224               write(licz3,"(bz,i3.3)") myparm
1225               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1226             else
1227               histname=prefix(:ilen(prefix))//'.hist'
1228             endif
1229           else
1230             if (separate_parset) then
1231               write(licz3,"(bz,i3.3)") myparm
1232               histname=prefix(:ilen(prefix))//'_par'//licz3//
1233      &         '_slice_'//licz2//'.hist'
1234             else
1235               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1236             endif
1237           endif
1238 #if defined(AIX) || defined(PGI)
1239           open (ihist,file=histname,position='append')
1240 #else
1241           open (ihist,file=histname,access='append')
1242 #endif
1243         endif
1244
1245         do t=0,tmax
1246           liczba=t
1247           sumH=0.0d0
1248           do ib=1,nT_h(iparm)
1249             sumH=sumH+hfin(t,ib)
1250           enddo
1251           if (sumH.gt.0.0d0) then
1252             do j=1,nQ
1253               jj = mod(liczba,nbin1)
1254               liczba=liczba/nbin1
1255               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1256               if (histfile)
1257      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1258             enddo
1259             do ib=1,nT_h(iparm)
1260               write (iout,'(e20.10,$)') hfin(t,ib)
1261               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1262             enddo
1263             write (iout,'(i5)') iparm
1264             if (histfile) write (ihist,'(i5)') iparm
1265           endif
1266         enddo
1267
1268         endif
1269
1270         if (entfile) then
1271           if (nslice.eq.1) then
1272             if (separate_parset) then
1273               write(licz3,"(bz,i3.3)") myparm
1274               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1275             else
1276               histname=prefix(:ilen(prefix))//'.ent'
1277             endif
1278           else
1279             if (separate_parset) then
1280               write(licz3,"(bz,i3.3)") myparm
1281               histname=prefix(:ilen(prefix))//'par_'//licz3//
1282      &           '_slice_'//licz2//'.ent'
1283             else
1284               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1285             endif
1286           endif
1287 #if defined(AIX) || defined(PGI)
1288           open (ihist,file=histname,position='append')
1289 #else
1290           open (ihist,file=histname,access='append')
1291 #endif
1292           write (ihist,'(a)') "# Microcanonical entropy"
1293           do i=0,upindE
1294             write (ihist,'(f8.0,$)') dint(potEmin)+i
1295             if (histE(i).gt.0.0e0) then
1296               write (ihist,'(f15.5,$)') dlog(histE(i))
1297             else
1298               write (ihist,'(f15.5,$)') 0.0d0
1299             endif
1300           enddo
1301           write (ihist,*)
1302           close(ihist)
1303         endif
1304         write (iout,*) "Microcanonical entropy"
1305         do i=0,upindE
1306           write (iout,'(f8.0,$)') dint(potEmin)+i
1307           if (histE(i).gt.0.0e0) then
1308             write (iout,'(f15.5,$)') dlog(histE(i))
1309           else
1310             write (iout,'(f15.5,$)') 0.0d0
1311           endif
1312           write (iout,*)
1313         enddo
1314         if (rmsrgymap) then
1315           if (nslice.eq.1) then
1316             if (separate_parset) then
1317               write(licz3,"(bz,i3.3)") myparm
1318               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1319             else
1320               histname=prefix(:ilen(prefix))//'.rmsrgy'
1321             endif
1322           else
1323             if (separate_parset) then
1324               write(licz3,"(bz,i3.3)") myparm
1325               histname=prefix(:ilen(prefix))//'_par'//licz3//
1326      &         '_slice_'//licz2//'.rmsrgy'
1327             else
1328              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1329             endif
1330           endif
1331 #if defined(AIX) || defined(PGI)
1332           open (ihist,file=histname,position='append')
1333 #else
1334           open (ihist,file=histname,access='append')
1335 #endif
1336           do i=0,nbin_rms
1337             do j=0,nbin_rgy
1338               write(ihist,'(2f8.2,$)')
1339      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1340               do ib=1,nT_h(iparm)
1341                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1342                   write(ihist,'(e14.5,$)')
1343      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1344      &              +potEmin
1345                 else
1346                   write(ihist,'(e14.5,$)') 1.0d6
1347                 endif
1348               enddo
1349               write (ihist,'(i2)') iparm
1350             enddo
1351           enddo
1352           close(ihist)
1353         endif
1354         endif
1355       enddo ! iparm
1356 #ifdef MPI
1357       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1358      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1359       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1360      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1361       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1362      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1363       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1364      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1365       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1366      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1367       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1368      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1369      &   WHAM_COMM,IERROR)
1370       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1371      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1372      &   WHAM_COMM,IERROR)
1373       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1374      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1375      &   WHAM_COMM,IERROR)
1376       if (me.eq.master) then
1377 #endif
1378       write (iout,'(/a)') 'Thermal characteristics of folding'
1379       if (nslice.eq.1) then
1380         nazwa=prefix
1381       else
1382         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1383       endif
1384       iln=ilen(nazwa)
1385       if (nparmset.eq.1 .and. .not.separate_parset) then
1386         nazwa=nazwa(:iln)//".thermal"
1387       else if (nparmset.eq.1 .and. separate_parset) then
1388         write(licz3,"(bz,i3.3)") myparm
1389         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1390       endif
1391       do iparm=1,nParmSet
1392       if (nparmset.gt.1) then
1393         write(licz3,"(bz,i3.3)") iparm
1394         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1395       endif
1396       open(34,file=nazwa)
1397       if (separate_parset) then
1398         write (iout,'(a,i3)') "Parameter set",myparm
1399       else
1400         write (iout,'(a,i3)') "Parameter set",iparm
1401       endif
1402 c      do i=0,NGridT
1403       do i=0,NGridT
1404         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1405         if (betaT.ge.beta_h(1,iparm)) then
1406           potEmin=potEmin_all(1,iparm)
1407         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1408           potEmin=potEmin_all(nT_h(iparm),iparm)
1409         else
1410           do l=1,nT_h(iparm)-1
1411             if (betaT.le.beta_h(l,iparm) .and.
1412      &          betaT.gt.beta_h(l+1,iparm)) then
1413               potEmin=potEmin_all(l,iparm)
1414               exit
1415             endif
1416           enddo
1417         endif
1418
1419 c        write (iout,*) "i",i," potEmin",potEmin
1420
1421         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1422         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1423      &    sumW(i,iparm)
1424         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1425      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1426         do j=1,nQ+2
1427           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1428           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1429      &     -sumQ(j,i,iparm)**2
1430           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1431      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1432         enddo
1433         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1434      &   (startGridT+i*delta_T))+potEmin
1435         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1436      &   sumW(i,iparm),sumE(i,iparm)
1437         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1438         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1439      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1440         write (iout,*)
1441         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1442      &   sumW(i,iparm),sumE(i,iparm)
1443         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1444         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1445      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1446         write (34,*)
1447       enddo
1448       close(34)
1449       enddo
1450       if (histout) then
1451       do t=0,tmax
1452         if (hfin_ent(t).gt.0.0d0) then
1453           liczba=t
1454           jj = mod(liczba,nbin1)
1455           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1456      &     hfin_ent(t)
1457           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1458      &      dmin+(jj+0.5d0)*delta,
1459      &     hfin_ent(t)
1460         endif
1461       enddo
1462       if (histfile) close(ihist)
1463       endif
1464
1465 #ifdef ZSCORE
1466 ! Write data for zscore
1467       if (nslice.eq.1) then
1468         zscname=prefix(:ilen(prefix))//".zsc"
1469       else
1470         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1471       endif
1472 #if defined(AIX) || defined(PGI)
1473       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1474 #else
1475       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1476 #endif
1477       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1478       do iparm=1,nParmSet
1479         write (izsc,'("NT=",i1)') nT_h(iparm)
1480       do ib=1,nT_h(iparm)
1481         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1482      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1483         jj = min0(nR(ib,iparm),7)
1484         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1485         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1486         write (izsc,'("&")')
1487         if (nR(ib,iparm).gt.7) then
1488           do ii=8,nR(ib,iparm),9
1489             jj = min0(nR(ib,iparm),ii+8)
1490             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1491             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1492             write (izsc,'("&")')
1493           enddo
1494         endif
1495         write (izsc,'("FI=",$)')
1496         jj=min0(nR(ib,iparm),7)
1497         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1498         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1499         write (izsc,'("&")')
1500         if (nR(ib,iparm).gt.7) then
1501           do ii=8,nR(ib,iparm),9
1502             jj = min0(nR(ib,iparm),ii+8)
1503             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1504             if (jj.eq.nR(ib,iparm)) then
1505               write (izsc,*)
1506             else
1507               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1508               write (izsc,'(t80,"&")')
1509             endif
1510           enddo
1511         endif
1512         do i=1,nR(ib,iparm)
1513           write (izsc,'("KH=",$)')
1514           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1515           write (izsc,'(" Q0=",$)')
1516           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1517           write (izsc,*)
1518         enddo
1519       enddo
1520       enddo
1521       close(izsc)
1522 #endif
1523 #ifdef MPI
1524       endif
1525 #endif
1526
1527       return
1528
1529       end