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