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