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      &      +wdfa_dist*edfadis
1090      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
1091             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1092      &            +ftprim(1)*wtor*etors+
1093      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1094      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1095      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1096      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1097      &            ftprim(1)*wsccor*esccor
1098             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1099      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1100      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1101      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1102      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1103      &            ftbis(1)*wsccor*esccor
1104 #else
1105             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1106      &      +ft(1)*welec*(ees+evdw1)
1107      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1108      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1109      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1110      &      +ft(2)*wturn3*eello_turn3
1111      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1112      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1113      &      +wbond*estr+ehomology_constr
1114      &      +wdfa_dist*edfadis
1115      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
1116             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1117      &           +ftprim(1)*wtor*etors+
1118      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1119      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1120      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1121      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1122      &            ftprim(1)*wsccor*esccor
1123             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1124      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1125      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1126      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1127      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1128      &            ftprim(1)*wsccor*esccor
1129 #endif
1130             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1131 #ifdef DEBUG
1132             write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
1133      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
1134      &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
1135      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
1136      &       wdfa_nei*edfanei+wdfa_beta*edfabet
1137             write (iout,*) "iparm",iparm," t",t," temper",temper,
1138      &        " etot",etot," entfac",entfac(t),
1139      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1140      &        " boltz",-betaT*(etot-potEmin)+entfac(t),
1141      &        " weight",weight," ebis",ebis
1142 #endif
1143             etot=etot-temper*eprim
1144 #ifdef MPI
1145             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1146             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1147             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1148             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1149             do j=1,nQ+2
1150               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1151               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1152               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1153      &         +etot*q(j,t)*weight
1154             enddo
1155 #else
1156             sumW(k,iparm)=sumW(k,iparm)+weight
1157             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1158             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1159             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1160             do j=1,nQ+2
1161               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1162               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1163               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1164      &         +etot*q(j,t)*weight
1165             enddo
1166 #endif
1167           enddo
1168           indE = aint(potE(t,iparm)-aint(potEmin))
1169           if (indE.ge.0 .and. indE.le.maxinde) then
1170             if (indE.gt.upindE_p) upindE_p=indE
1171             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1172           endif
1173 #ifdef MPI
1174           do ib=1,nT_h(iparm)
1175             potEmin=potEmin_all(ib,iparm)
1176             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1177             hfin_p(ind,ib)=hfin_p(ind,ib)+
1178      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1179              if (rmsrgymap) then
1180                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1181                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1182                hrmsrgy_p(indrgy,indrms,ib)=
1183      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1184              endif
1185           enddo
1186 #else
1187           do ib=1,nT_h(iparm)
1188             potEmin=potEmin_all(ib,iparm)
1189             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1190             hfin(ind,ib)=hfin(ind,ib)+
1191      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1192              if (rmsrgymap) then
1193                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1194                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1195                hrmsrgy(indrgy,indrms,ib)=
1196      &           hrmsrgy(indrgy,indrms,ib)+expfac
1197              endif
1198           enddo
1199 #endif
1200         enddo ! t
1201         do ib=1,nT_h(iparm)
1202           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1203      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1204           if (rmsrgymap) then
1205           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1206      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1207      &       WHAM_COMM,IERROR)
1208           endif
1209         enddo
1210         call MPI_Reduce(upindE_p,upindE,1,
1211      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1212         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1213      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1214
1215         if (me1.eq.master) then
1216
1217         if (histout) then
1218
1219         write (iout,'(6x,$)')
1220         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1221      &   ib=1,nT_h(iparm))
1222         write (iout,*)
1223
1224         write (iout,'(/a)') 'Final histograms'
1225         if (histfile) then
1226           if (nslice.eq.1) then
1227             if (separate_parset) then
1228               write(licz3,"(bz,i3.3)") myparm
1229               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1230             else
1231               histname=prefix(:ilen(prefix))//'.hist'
1232             endif
1233           else
1234             if (separate_parset) then
1235               write(licz3,"(bz,i3.3)") myparm
1236               histname=prefix(:ilen(prefix))//'_par'//licz3//
1237      &         '_slice_'//licz2//'.hist'
1238             else
1239               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1240             endif
1241           endif
1242 #if defined(AIX) || defined(PGI)
1243           open (ihist,file=histname,position='append')
1244 #else
1245           open (ihist,file=histname,access='append')
1246 #endif
1247         endif
1248
1249         do t=0,tmax
1250           liczba=t
1251           sumH=0.0d0
1252           do ib=1,nT_h(iparm)
1253             sumH=sumH+hfin(t,ib)
1254           enddo
1255           if (sumH.gt.0.0d0) then
1256             do j=1,nQ
1257               jj = mod(liczba,nbin1)
1258               liczba=liczba/nbin1
1259               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1260               if (histfile)
1261      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1262             enddo
1263             do ib=1,nT_h(iparm)
1264               write (iout,'(e20.10,$)') hfin(t,ib)
1265               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1266             enddo
1267             write (iout,'(i5)') iparm
1268             if (histfile) write (ihist,'(i5)') iparm
1269           endif
1270         enddo
1271
1272         endif
1273
1274         if (entfile) then
1275           if (nslice.eq.1) then
1276             if (separate_parset) then
1277               write(licz3,"(bz,i3.3)") myparm
1278               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1279             else
1280               histname=prefix(:ilen(prefix))//'.ent'
1281             endif
1282           else
1283             if (separate_parset) then
1284               write(licz3,"(bz,i3.3)") myparm
1285               histname=prefix(:ilen(prefix))//'par_'//licz3//
1286      &           '_slice_'//licz2//'.ent'
1287             else
1288               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1289             endif
1290           endif
1291 #if defined(AIX) || defined(PGI)
1292           open (ihist,file=histname,position='append')
1293 #else
1294           open (ihist,file=histname,access='append')
1295 #endif
1296           write (ihist,'(a)') "# Microcanonical entropy"
1297           do i=0,upindE
1298             write (ihist,'(f8.0,$)') dint(potEmin)+i
1299             if (histE(i).gt.0.0e0) then
1300               write (ihist,'(f15.5,$)') dlog(histE(i))
1301             else
1302               write (ihist,'(f15.5,$)') 0.0d0
1303             endif
1304           enddo
1305           write (ihist,*)
1306           close(ihist)
1307         endif
1308         write (iout,*) "Microcanonical entropy"
1309         do i=0,upindE
1310           write (iout,'(f8.0,$)') dint(potEmin)+i
1311           if (histE(i).gt.0.0e0) then
1312             write (iout,'(f15.5,$)') dlog(histE(i))
1313           else
1314             write (iout,'(f15.5,$)') 0.0d0
1315           endif
1316           write (iout,*)
1317         enddo
1318         if (rmsrgymap) then
1319           if (nslice.eq.1) then
1320             if (separate_parset) then
1321               write(licz3,"(bz,i3.3)") myparm
1322               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1323             else
1324               histname=prefix(:ilen(prefix))//'.rmsrgy'
1325             endif
1326           else
1327             if (separate_parset) then
1328               write(licz3,"(bz,i3.3)") myparm
1329               histname=prefix(:ilen(prefix))//'_par'//licz3//
1330      &         '_slice_'//licz2//'.rmsrgy'
1331             else
1332              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1333             endif
1334           endif
1335 #if defined(AIX) || defined(PGI)
1336           open (ihist,file=histname,position='append')
1337 #else
1338           open (ihist,file=histname,access='append')
1339 #endif
1340           do i=0,nbin_rms
1341             do j=0,nbin_rgy
1342               write(ihist,'(2f8.2,$)')
1343      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1344               do ib=1,nT_h(iparm)
1345                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1346                   write(ihist,'(e14.5,$)')
1347      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1348      &              +potEmin
1349                 else
1350                   write(ihist,'(e14.5,$)') 1.0d6
1351                 endif
1352               enddo
1353               write (ihist,'(i2)') iparm
1354             enddo
1355           enddo
1356           close(ihist)
1357         endif
1358         endif
1359       enddo ! iparm
1360 #ifdef MPI
1361       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1362      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1363       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1364      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1365       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1366      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1367       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1368      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1369       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1370      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1371       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1372      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1373      &   WHAM_COMM,IERROR)
1374       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1375      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1376      &   WHAM_COMM,IERROR)
1377       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1378      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1379      &   WHAM_COMM,IERROR)
1380       if (me.eq.master) then
1381 #endif
1382       write (iout,'(/a)') 'Thermal characteristics of folding'
1383       if (nslice.eq.1) then
1384         nazwa=prefix
1385       else
1386         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1387       endif
1388       iln=ilen(nazwa)
1389       if (nparmset.eq.1 .and. .not.separate_parset) then
1390         nazwa=nazwa(:iln)//".thermal"
1391       else if (nparmset.eq.1 .and. separate_parset) then
1392         write(licz3,"(bz,i3.3)") myparm
1393         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1394       endif
1395       do iparm=1,nParmSet
1396       if (nparmset.gt.1) then
1397         write(licz3,"(bz,i3.3)") iparm
1398         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1399       endif
1400       open(34,file=nazwa)
1401       if (separate_parset) then
1402         write (iout,'(a,i3)') "Parameter set",myparm
1403       else
1404         write (iout,'(a,i3)') "Parameter set",iparm
1405       endif
1406 c      do i=0,NGridT
1407       do i=0,NGridT
1408         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1409         if (betaT.ge.beta_h(1,iparm)) then
1410           potEmin=potEmin_all(1,iparm)
1411         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1412           potEmin=potEmin_all(nT_h(iparm),iparm)
1413         else
1414           do l=1,nT_h(iparm)-1
1415             if (betaT.le.beta_h(l,iparm) .and.
1416      &          betaT.gt.beta_h(l+1,iparm)) then
1417               potEmin=potEmin_all(l,iparm)
1418               exit
1419             endif
1420           enddo
1421         endif
1422
1423 c        write (iout,*) "i",i," potEmin",potEmin
1424
1425         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1426         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1427      &    sumW(i,iparm)
1428         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1429      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1430         do j=1,nQ+2
1431           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1432           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1433      &     -sumQ(j,i,iparm)**2
1434           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1435      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1436         enddo
1437         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1438      &   (startGridT+i*delta_T))+potEmin
1439         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1440      &   sumW(i,iparm),sumE(i,iparm)
1441         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1442         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1443      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1444         write (iout,*)
1445         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1446      &   sumW(i,iparm),sumE(i,iparm)
1447         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1448         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1449      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1450         write (34,*)
1451       enddo
1452       close(34)
1453       enddo
1454       if (histout) then
1455       do t=0,tmax
1456         if (hfin_ent(t).gt.0.0d0) then
1457           liczba=t
1458           jj = mod(liczba,nbin1)
1459           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1460      &     hfin_ent(t)
1461           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1462      &      dmin+(jj+0.5d0)*delta,
1463      &     hfin_ent(t)
1464         endif
1465       enddo
1466       if (histfile) close(ihist)
1467       endif
1468
1469 #ifdef ZSCORE
1470 ! Write data for zscore
1471       if (nslice.eq.1) then
1472         zscname=prefix(:ilen(prefix))//".zsc"
1473       else
1474         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1475       endif
1476 #if defined(AIX) || defined(PGI)
1477       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1478 #else
1479       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1480 #endif
1481       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1482       do iparm=1,nParmSet
1483         write (izsc,'("NT=",i1)') nT_h(iparm)
1484       do ib=1,nT_h(iparm)
1485         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1486      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1487         jj = min0(nR(ib,iparm),7)
1488         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1489         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1490         write (izsc,'("&")')
1491         if (nR(ib,iparm).gt.7) then
1492           do ii=8,nR(ib,iparm),9
1493             jj = min0(nR(ib,iparm),ii+8)
1494             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1495             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1496             write (izsc,'("&")')
1497           enddo
1498         endif
1499         write (izsc,'("FI=",$)')
1500         jj=min0(nR(ib,iparm),7)
1501         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1502         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1503         write (izsc,'("&")')
1504         if (nR(ib,iparm).gt.7) then
1505           do ii=8,nR(ib,iparm),9
1506             jj = min0(nR(ib,iparm),ii+8)
1507             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1508             if (jj.eq.nR(ib,iparm)) then
1509               write (izsc,*)
1510             else
1511               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1512               write (izsc,'(t80,"&")')
1513             endif
1514           enddo
1515         endif
1516         do i=1,nR(ib,iparm)
1517           write (izsc,'("KH=",$)')
1518           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1519           write (izsc,'(" Q0=",$)')
1520           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1521           write (izsc,*)
1522         enddo
1523       enddo
1524       enddo
1525       close(izsc)
1526 #endif
1527 #ifdef MPI
1528       endif
1529 #endif
1530
1531       return
1532
1533       end