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