X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fwham_calc1.F.safe;fp=source%2Fwham%2Fsrc-NEWSC%2Fwham_calc1.F.safe;h=f51dcc45145d2b1ea230960066142abe887a92f4;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/wham_calc1.F.safe b/source/wham/src-NEWSC/wham_calc1.F.safe new file mode 100755 index 0000000..f51dcc4 --- /dev/null +++ b/source/wham/src-NEWSC/wham_calc1.F.safe @@ -0,0 +1,1195 @@ + subroutine WHAM_CALC(islice,*) +! Weighed Histogram Analysis Method (WHAM) code +! Written by A. Liwo based on the work of Kumar et al., +! J.Comput.Chem., 13, 1011 (1992) +! +! 2/1/05 Multiple temperatures allowed. +! 2/2/05 Free energies calculated directly from data points +! acc. to Eq. (21) of Kumar et al.; final histograms also +! constructed based on this equation. +! 2/12/05 Multiple parameter sets included +! +! 2/2/05 Parallel version + implicit none + include "DIMENSIONS" + include "DIMENSIONS.ZSCOPT" + include "DIMENSIONS.FREE" + integer nGridT + parameter (NGridT=400) + integer MaxBinRms,MaxBinRgy + parameter (MaxBinRms=100,MaxBinRgy=100) + integer MaxHdim +c parameter (MaxHdim=200000) + parameter (MaxHdim=200) + integer maxinde + parameter (maxinde=200) +#ifdef MPI + include "mpif.h" + include "COMMON.MPI" + integer ierror,errcode,status(MPI_STATUS_SIZE) +#endif + include "COMMON.CONTROL" + include "COMMON.IOUNITS" + include "COMMON.FREE" + include "COMMON.ENERGIES" + include "COMMON.FFIELD" + include "COMMON.SBRIDGE" + include "COMMON.PROT" + include "COMMON.ENEPS" + integer MaxPoint,MaxPointProc + parameter (MaxPoint=MaxStr, + & MaxPointProc=MaxStr_Proc) + double precision finorm_max,potfac,entmin,entmax,expfac,vf + parameter (finorm_max=1.0d0) + integer islice + integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln + integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy, + & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms + integer htot(0:MaxHdim),histent(0:2000) + double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm) + double precision energia(0:max_ene) +#ifdef MPI + integer tmax_t,upindE_p + double precision fi_p(MaxR,MaxT_h,Max_Parm) + double precision sumW_p(0:nGridT,Max_Parm), + & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm), + & sumQ_p(MaxQ1,0:nGridT,Max_Parm), + & sumQsq_p(MaxQ1,0:nGridT,Max_Parm), + & sumEQ_p(MaxQ1,0:nGridT,Max_Parm), + & sumEprim_p(MaxQ1,0:nGridT,Max_Parm), + & sumEbis_p(0:nGridT,Max_Parm) + double precision hfin_p(0:MaxHdim,maxT_h), + & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH, + & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h) + double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t + double precision potEmin_t,entmin_p,entmax_p + integer histent_p(0:2000) + logical lprint /.true./ +#endif + double precision delta_T /1.0d0/ + double precision rgymin,rmsmin,rgymax,rmsmax + double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm), + & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm), + & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm), + & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT, + & weight,econstr + double precision fi(MaxR,maxT_h,Max_Parm), + & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom, + & hfin(0:MaxHdim,maxT_h),histE(0:maxindE), + & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h), + & potEmin,ent, + & hfin_ent(0:MaxHdim),vmax,aux + double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, + & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/, + & eplus,eminus,logfac,tanhT,tt + double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, + & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, + & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor + + integer ind_point(maxpoint),upindE,indE + character*16 plik + character*1 licz1 + character*2 licz2 + character*3 licz3 + character*128 nazwa + integer ilen + external ilen + + write(licz2,'(bz,i2.2)') islice + nbin1 = 1.0d0/delta + write (iout,'(//80(1h-)/"Solving WHAM equations for slice", + & i2/80(1h-)//)') islice + write (iout,*) "delta",delta," nbin1",nbin1 + write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim + call flush(iout) + dmin=0.0d0 + tmax=0 + potEmin=1.0d10 + rgymin=1.0d10 + rmsmin=1.0d10 + rgymax=0.0d0 + rmsmax=0.0d0 + do t=0,MaxN + htot(t)=0 + enddo +#ifdef MPI + do i=1,scount(me1) +#else + do i=1,ntot(islice) +#endif + do j=1,nParmSet + if (potE(i,j).le.potEmin) potEmin=potE(i,j) + enddo + if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i) + if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i) + if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i) + if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i) + ind_point(i)=0 + do j=nQ,1,-1 + ind=(q(j,i)-dmin+1.0d-8)/delta + if (j.eq.1) then + ind_point(i)=ind_point(i)+ind + else + ind_point(i)=ind_point(i)+nbin1**(j-1)*ind + endif +c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point", +c & ind_point(i) + call flush(iout) + if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then + write (iout,*) "Error - index exceeds range for point",i, + & " q=",q(j,i)," ind",ind_point(i) +#ifdef MPI + write (iout,*) "Processor",me1 + call flush(iout) + call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode ) +#endif + stop + endif + enddo ! j + if (ind_point(i).gt.tmax) tmax=ind_point(i) + htot(ind_point(i))=htot(ind_point(i))+1 +#ifdef DEBUG + write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i), + & " htot",htot(ind_point(i)) + call flush(iout) +#endif + enddo ! i + call flush(iout) + + nbin=nbin1**nQ-1 + write (iout,'(a)') "Numbers of counts in Q bins" + do t=0,tmax + if (htot(t).gt.0) then + write (iout,'(i15,$)') t + liczba=t + do j=1,nQ + jj = mod(liczba,nbin1) + liczba=liczba/nbin1 + write (iout,'(i5,$)') jj + enddo + write (iout,'(i8)') htot(t) + endif + enddo + do iparm=1,nParmSet + write (iout,'(a,i3)') "Number of data points for parameter set", + & iparm + write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)), + & ib=1,nT_h(iparm)) + write (iout,'(i8)') stot(islice) + write (iout,'(a)') + enddo + call flush(iout) + +#ifdef MPI + call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX, + & WHAM_COMM,IERROR) + tmax=tmax_t + call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION, + & MPI_MIN,WHAM_COMM,IERROR) + call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION, + & MPI_MIN,WHAM_COMM,IERROR) + call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION, + & MPI_MAX,WHAM_COMM,IERROR) + call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION, + & MPI_MIN,WHAM_COMM,IERROR) + call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION, + & MPI_MAX,WHAM_COMM,IERROR) + potEmin=potEmin_t/2 + rgymin=rgymin_t + rgymax=rgymax_t + rmsmin=rmsmin_t + rmsmax=rmsmax_t + write (iout,*) "potEmin",potEmin +#endif + rmsmin=deltrms*dint(rmsmin/deltrms) + rmsmax=deltrms*dint(rmsmax/deltrms) + rgymin=deltrms*dint(rgymin/deltrgy) + rgymax=deltrms*dint(rgymax/deltrgy) + nbin_rms=(rmsmax-rmsmin)/deltrms + nbin_rgy=(rgymax-rgymin)/deltrgy + write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin, + & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy + nFi=0 + do i=1,nParmSet + do j=1,nT_h(i) + nFi=nFi+nR(j,i) + enddo + enddo + write (iout,*) "nFi",nFi +! Compute the Boltzmann factor corresponing to restrain potentials in different +! simulations. +#ifdef MPI + do i=1,scount(me1) +#else + do i=1,ntot(islice) +#endif +c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) + do iparm=1,nParmSet +#ifdef DEBUG + write (iout,'(2i5,21f8.2)') i,iparm, + & (enetb(k,i,iparm),k=1,21) +#endif + call restore_parm(iparm) +#ifdef DEBUG + write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, + & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, + & wtor_d,wsccor,wbond +#endif + do ib=1,nT_h(iparm) + if (rescale_mode.eq.1) then + quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) + quotl=1.0d0 + kfacl=1.0d0 + do l=1,5 + quotl1=quotl + quotl=quotl*quot + kfacl=kfacl*kfac + fT(l)=kfacl/(kfacl-1.0d0+quotl) + enddo +#if defined(FUNCTH) + tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3) + ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 +#elif defined(FUNCT) + ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) +#else + ft(6)=1.0d0 +#endif + else if (rescale_mode.eq.2) then + quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) + quotl=1.0d0 + do l=1,5 + quotl=quotl*quot + fT(l)=1.12692801104297249644d0/ + & dlog(dexp(quotl)+dexp(-quotl)) + enddo +#if defined(FUNCTH) + tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3) + ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 +#elif defined(FUNCT) + ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0) +#else + ft(6)=1.0d0 +#endif +c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft + else if (rescale_mode.eq.0) then + do l=1,6 + fT(l)=1.0d0 + enddo + else + write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", + & rescale_mode + call flush(iout) + return1 + endif + evdw=enetb(1,i,iparm) + evdw_t=enetb(21,i,iparm) +#ifdef SCP14 + evdw2_14=enetb(17,i,iparm) + evdw2=enetb(2,i,iparm)+evdw2_14 +#else + evdw2=enetb(2,i,iparm) + evdw2_14=0.0d0 +#endif +#ifdef SPLITELE + ees=enetb(3,i,iparm) + evdw1=enetb(16,i,iparm) +#else + ees=enetb(3,i,iparm) + evdw1=0.0d0 +#endif + ecorr=enetb(4,i,iparm) + ecorr5=enetb(5,i,iparm) + ecorr6=enetb(6,i,iparm) + eel_loc=enetb(7,i,iparm) + eello_turn3=enetb(8,i,iparm) + eello_turn4=enetb(9,i,iparm) + eturn6=enetb(10,i,iparm) + ebe=enetb(11,i,iparm) + escloc=enetb(12,i,iparm) + etors=enetb(13,i,iparm) + etors_d=enetb(14,i,iparm) + ehpb=enetb(15,i,iparm) + estr=enetb(18,i,iparm) + esccor=enetb(19,i,iparm) + edihcnstr=enetb(20,i,iparm) +#ifdef DEBUG + write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), + & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc, + & etors,etors_d,eello_turn3,eello_turn4,esccor +#endif + +#ifdef SPLITELE + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees + & +wvdwpp*evdw1 + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc + & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr +#else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 + & +ft(1)*welec*(ees+evdw1) + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr + & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr +#endif +#ifdef DEBUG + write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3), + & etot,potEmin +#endif +#ifdef DEBUG + if (iparm.eq.1 .and. ib.eq.1) then + write (iout,*)"Conformation",i + energia(0)=etot + do k=1,max_ene + energia(k)=enetb(k,i,iparm) + enddo + call enerprint(energia(0),fT) + endif +#endif + do kk=1,nR(ib,iparm) + Econstr=0.0d0 + do j=1,nQ + dd = q(j,i) + Econstr=Econstr+Kh(j,kk,ib,iparm) + & *(dd-q0(j,kk,ib,iparm))**2 + enddo + v(i,kk,ib,iparm)= + & -beta_h(ib,iparm)*(etot-potEmin+Econstr) +#ifdef DEBUG + write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, + & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm) +#endif + enddo ! kk + enddo ! ib + enddo ! iparm + enddo ! i +! Simple iteration to calculate free energies corresponding to all simulation +! runs. + do iter=1,maxit + +! Compute new free-energy values corresponding to the righ-hand side of the +! equation and their derivatives. + write (iout,*) "------------------------fi" +#ifdef MPI + do t=1,scount(me1) +#else + do t=1,ntot(islice) +#endif + vmax=-1.0d+20 + do i=1,nParmSet + do k=1,nT_h(i) + do l=1,nR(k,i) + vf=v(t,l,k,i)+f(l,k,i) + if (vf.gt.vmax) vmax=vf + enddo + enddo + enddo + denom=0.0d0 + do i=1,nParmSet + do k=1,nT_h(i) + do l=1,nR(k,i) + aux=f(l,k,i)+v(t,l,k,i)-vmax + if (aux.gt.-200.0d0) + & denom=denom+snk(l,k,i,islice)*dexp(aux) + enddo + enddo + enddo + entfac(t)=-dlog(denom)-vmax +#ifdef DEBUG + write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t) +#endif + enddo + do iparm=1,nParmSet + do iib=1,nT_h(iparm) + do ii=1,nR(iib,iparm) +#ifdef MPI + fi_p(ii,iib,iparm)=0.0d0 + do t=1,scount(me) + fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) + & +dexp(v(t,ii,iib,iparm)+entfac(t)) +#ifdef DEBUG + write (iout,'(4i5,3e15.5)') t,ii,iib,iparm, + & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm) +#endif + enddo +#else + fi(ii,iib,iparm)=0.0d0 + do t=1,ntot(islice) + fi(ii,iib,iparm)=fi(ii,iib,iparm) + & +dexp(v(t,ii,iib,iparm)+entfac(t)) + enddo +#endif + enddo ! ii + enddo ! iib + enddo ! iparm + +#ifdef MPI +#ifdef DEBUG + write (iout,*) "fi before MPI_Reduce me",me,' master',master + do iparm=1,nParmSet + do ib=1,nT_h(nparmset) + write (iout,*) "iparm",iparm," ib",ib + write (iout,*) "beta=",beta_h(ib,iparm) + write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm)) + enddo + enddo +#endif + write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, + & maxR*MaxT_h*nParmSet + write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, + & " WHAM_COMM",WHAM_COMM + call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet, + & MPI_DOUBLE_PRECISION, + & MPI_SUM,Master,WHAM_COMM,IERROR) +#ifdef DEBUG + write (iout,*) "fi after MPI_Reduce nparmset",nparmset + do iparm=1,nParmSet + write (iout,*) "iparm",iparm + do ib=1,nT_h(iparm) + write (iout,*) "beta=",beta_h(ib,iparm) + write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm)) + enddo + enddo +#endif + if (me1.eq.Master) then +#endif + avefi=0.0d0 + do iparm=1,nParmSet + do ib=1,nT_h(iparm) + do i=1,nR(ib,iparm) + fi(i,ib,iparm)=-dlog(fi(i,ib,iparm)) + avefi=avefi+fi(i,ib,iparm) + enddo + enddo + enddo + avefi=avefi/nFi + do iparm=1,nParmSet + write (iout,*) "Parameter set",iparm + do ib =1,nT_h(iparm) + write (iout,*) "beta=",beta_h(ib,iparm) + do i=1,nR(ib,iparm) + fi(i,ib,iparm)=fi(i,ib,iparm)-avefi + enddo + write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm)) + write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm)) + enddo + enddo + +! Compute the norm of free-energy increments. + finorm=0.0d0 + do iparm=1,nParmSet + do ib=1,nT_h(iparm) + do i=1,nR(ib,iparm) + finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm)) + f(i,ib,iparm)=fi(i,ib,iparm) + enddo + enddo + enddo + + write (iout,*) 'Iteration',iter,' finorm',finorm + +#ifdef MPI + endif + call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet, + & MPI_DOUBLE_PRECISION,Master, + & WHAM_COMM,IERROR) + call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master, + & WHAM_COMM,IERROR) +#endif +! Exit, if the increment norm is smaller than pre-assigned tolerance. + if (finorm.lt.fimin) then + write (iout,*) 'Iteration converged' + goto 20 + endif + + enddo ! iter + + 20 continue +! Now, put together the histograms from all simulations, in order to get the +! unbiased total histogram. +#ifdef MPI + do t=0,tmax + hfin_ent_p(t)=0.0d0 + enddo +#else + do t=0,tmax + hfin_ent(t)=0.0d0 + enddo +#endif + write (iout,*) "--------------hist" +#ifdef MPI + do iparm=1,nParmSet + do i=0,nGridT + sumW_p(i,iparm)=0.0d0 + sumE_p(i,iparm)=0.0d0 + sumEbis_p(i,iparm)=0.0d0 + sumEsq_p(i,iparm)=0.0d0 + do j=1,nQ+2 + sumQ_p(j,i,iparm)=0.0d0 + sumQsq_p(j,i,iparm)=0.0d0 + sumEQ_p(j,i,iparm)=0.0d0 + enddo + enddo + enddo + upindE_p=0 +#else + do iparm=1,nParmSet + do i=0,nGridT + sumW(i,iparm)=0.0d0 + sumE(i,iparm)=0.0d0 + sumEbis(i,iparm)=0.0d0 + sumEsq(i,iparm)=0.0d0 + do j=1,nQ+2 + sumQ(j,i,iparm)=0.0d0 + sumQsq(j,i,iparm)=0.0d0 + sumEQ(j,i,iparm)=0.0d0 + enddo + enddo + enddo + upindE=0 +#endif +c 8/26/05 entropy distribution +#ifdef MPI + entmin_p=1.0d10 + entmax_p=-1.0d10 + do t=1,scount(me1) +c ent=-dlog(entfac(t)) + ent=entfac(t) + if (ent.lt.entmin_p) entmin_p=ent + if (ent.gt.entmax_p) entmax_p=ent + enddo + write (iout,*) "entmin",entmin_p," entmax",entmax_p + call flush(iout) + call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN, + & WHAM_COMM,IERROR) + call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX, + & WHAM_COMM,IERROR) + ientmax=entmax-entmin + if (ientmax.gt.2000) ientmax=2000 + write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax + call flush(iout) + do t=1,scount(me1) +c ient=-dlog(entfac(t))-entmin + ient=entfac(t)-entmin + if (ient.le.2000) histent_p(ient)=histent_p(ient)+1 + enddo + call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER, + & MPI_SUM,WHAM_COMM,IERROR) + if (me1.eq.Master) then + write (iout,*) "Entropy histogram" + do i=0,ientmax + write(iout,'(f15.4,i10)') entmin+i,histent(i) + enddo + endif +#else + entmin=1.0d10 + entmax=-1.0d10 + do t=1,ntot(islice) + ent=entfac(t) + if (ent.lt.entmin) entmin=ent + if (ent.gt.entmax) entmax=ent + enddo + ientmax=-dlog(entmax)-entmin + if (ientmax.gt.2000) ientmax=2000 + do t=1,ntot(islice) + ient=entfac(t)-entmin + if (ient.le.2000) histent(ient)=histent(ient)+1 + enddo + write (iout,*) "Entropy histogram" + do i=0,ientmax + write(iout,'(2f15.4)') entmin+i,histent(i) + enddo +#endif + +#ifdef MPI +c write (iout,*) "me1",me1," scount",scount(me1) + + do iparm=1,nParmSet + +#ifdef MPI + do ib=1,nT_h(iparm) + do t=0,tmax + hfin_p(t,ib)=0.0d0 + enddo + enddo + do i=1,maxindE + histE_p(i)=0.0d0 + enddo +#else + do ib=1,nT_h(iparm) + do t=0,tmax + hfin(t,ib)=0.0d0 + enddo + enddo + do i=1,maxindE + histE(i)=0.0d0 + enddo +#endif + do ib=1,nT_h(iparm) + do i=0,MaxBinRms + do j=0,MaxBinRgy + hrmsrgy(j,i,ib)=0.0d0 +#ifdef MPI + hrmsrgy_p(j,i,ib)=0.0d0 +#endif + enddo + enddo + enddo + + do t=1,scount(me1) +#else + do t=1,ntot(islice) +#endif + ind = ind_point(t) +#ifdef MPI + hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t)) +#else + hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t)) +#endif +c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18) + call restore_parm(iparm) + evdw=enetb(21,t,iparm) + evdw_t=enetb(1,t,iparm) +#ifdef SCP14 + evdw2_14=enetb(17,t,iparm) + evdw2=enetb(2,t,iparm)+evdw2_14 +#else + evdw2=enetb(2,t,iparm) + evdw2_14=0.0d0 +#endif +#ifdef SPLITELE + ees=enetb(3,t,iparm) + evdw1=enetb(16,t,iparm) +#else + ees=enetb(3,t,iparm) + evdw1=0.0d0 +#endif + ecorr=enetb(4,t,iparm) + ecorr5=enetb(5,t,iparm) + ecorr6=enetb(6,t,iparm) + eel_loc=enetb(7,t,iparm) + eello_turn3=enetb(8,t,iparm) + eello_turn4=enetb(9,t,iparm) + eturn6=enetb(10,t,iparm) + ebe=enetb(11,t,iparm) + escloc=enetb(12,t,iparm) + etors=enetb(13,t,iparm) + etors_d=enetb(14,t,iparm) + ehpb=enetb(15,t,iparm) + estr=enetb(18,t,iparm) + esccor=enetb(19,t,iparm) + edihcnstr=enetb(20,t,iparm) + edihcnstr=0.0d0 + do k=0,nGridT + betaT=startGridT+k*delta_T + temper=betaT +c fT=T0/betaT +c ft=2*T0/(T0+betaT) + if (rescale_mode.eq.1) then + quot=betaT/T0 + quotl=1.0d0 + kfacl=1.0d0 + do l=1,5 + quotl1=quotl + quotl=quotl*quot + kfacl=kfacl*kfac + denom=kfacl-1.0d0+quotl + fT(l)=kfacl/denom + ftprim(l)=-l*ft(l)*quotl1/(T0*denom) + ftbis(l)=l*kfacl*quotl1* + & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3) + enddo +#if defined(FUNCTH) + ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ + & 320.0d0 + ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2) + ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) + & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3) +#elif defined(FUNCT) + fT(6)=betaT/T0 + ftprim(6)=1.0d0/T0 + ftbis(6)=0.0d0 +#else + fT(6)=1.0d0 + ftprim(6)=0.0d0 + ftbis(6)=0.0d0 +#endif + else if (rescale_mode.eq.2) then + quot=betaT/T0 + quotl=1.0d0 + do l=1,5 + quotl1=quotl + quotl=quotl*quot + eplus=dexp(quotl) + eminus=dexp(-quotl) + logfac=1.0d0/dlog(eplus+eminus) + tanhT=(eplus-eminus)/(eplus+eminus) + fT(l)=1.12692801104297249644d0*logfac + ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0 + ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- + & 2*l*quotl1/T0*logfac* + & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) + & +ftprim(l)*tanhT) + enddo +#if defined(FUNCTH) + ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ + & 320.0d0 + ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2) + ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) + & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3) +#elif defined(FUNCT) + fT(6)=betaT/T0 + ftprim(6)=1.0d0/T0 + ftbis(6)=0.0d0 +#else + fT(6)=1.0d0 + ftprim(6)=0.0d0 + ftbis(6)=0.0d0 +#endif + else if (rescale_mode.eq.0) then + do l=1,5 + fT(l)=1.0d0 + ftprim(l)=0.0d0 + enddo + else + write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", + & rescale_mode + call flush(iout) + return1 + endif +c write (iout,*) "ftprim",ftprim +c write (iout,*) "ftbis",ftbis + betaT=1.0d0/(1.987D-3*betaT) +#ifdef SPLITELE + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees + & +wvdwpp*evdw1 + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc + & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr + eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees + & +ftprim(1)*wtor*etors+ + & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ + & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ + & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+ + & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ + & ftprim(1)*wsccor*esccor + ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ + & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ + & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ + & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+ + & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ + & ftbis(1)*wsccor*esccor +#else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 + & +ft(1)*welec*(ees+evdw1) + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr + & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr + eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) + & +ftprim(1)*wtor*etors+ + & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ + & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ + & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+ + & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ + & ftprim(1)*wsccor*esccor + ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ + & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ + & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ + & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+ + & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ + & ftprim(1)*wsccor*esccor +#endif + weight=dexp(-betaT*(etot-potEmin)+entfac(t)) +#define DEBUG +#ifdef DEBUG + write (iout,*) "iparm",iparm," t",t," betaT",betaT, + & " etot",etot," entfac",entfac(t), + & " weight",weight," ebis",ebis +#endif +#undef DEBUG + etot=etot-temper*eprim +#ifdef MPI + sumW_p(k,iparm)=sumW_p(k,iparm)+weight + sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight + sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight + sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight + do j=1,nQ+2 + sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight + sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight + sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) + & +etot*q(j,t)*weight + enddo +#else + sumW(k,iparm)=sumW(k,iparm)+weight + sumE(k,iparm)=sumE(k,iparm)+etot*weight + sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight + sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight + do j=1,nQ+2 + sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight + sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight + sumEQ(j,k,iparm)=sumEQ(j,k,iparm) + & +etot*q(j,t)*weight + enddo +#endif + enddo + indE = aint(potE(t,iparm)-aint(potEmin)) + if (indE.ge.0 .and. indE.le.maxinde) then + if (indE.gt.upindE_p) upindE_p=indE + histE_p(indE)=histE_p(indE)+dexp(-entfac(t)) + endif +#ifdef MPI + do ib=1,nT_h(iparm) + expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) + hfin_p(ind,ib)=hfin_p(ind,ib)+ + & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) + if (rmsrgymap) then + indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) + indrms=dint((q(nQ+1,t)-rmsmin)/deltrms) + hrmsrgy_p(indrgy,indrms,ib)= + & hrmsrgy_p(indrgy,indrms,ib)+expfac + endif + enddo +#else + do ib=1,nT_h(iparm) + expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) + hfin(ind,ib)=hfin(ind,ib)+ + & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t)) + if (rmsrgymap) then + indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) + indrms=dint((q(nQ+1,t)-rmsmin)/deltrms) + hrmsrgy(indrgy,indrms,ib)= + & hrmsrgy(indrgy,indrms,ib)+expfac + endif + enddo +#endif + enddo ! t + do ib=1,nT_h(iparm) + if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + if (rmsrgymap) then + call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib), + & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master, + & WHAM_COMM,IERROR) + endif + enddo + call MPI_Reduce(upindE_p,upindE,1, + & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR) + call MPI_Reduce(histE_p(0),histE(0),maxindE, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + + if (me1.eq.master) then + + if (histout) then + + write (iout,'(6x,$)') + write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)), + & ib=1,nT_h(iparm)) + write (iout,*) + + write (iout,'(/a)') 'Final histograms' + if (histfile) then + if (nslice.eq.1) then + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist' + else + histname=prefix(:ilen(prefix))//'.hist' + endif + else + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//'_par'//licz3// + & '_slice_'//licz2//'.hist' + else + histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist' + endif + endif +#if defined(AIX) || defined(PGI) + open (ihist,file=histname,position='append') +#else + open (ihist,file=histname,access='append') +#endif + endif + + do t=0,tmax + liczba=t + sumH=0.0d0 + do ib=1,nT_h(iparm) + sumH=sumH+hfin(t,ib) + enddo + if (sumH.gt.0.0d0) then + do j=1,nQ + jj = mod(liczba,nbin1) + liczba=liczba/nbin1 + write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta + if (histfile) + & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta + enddo + do ib=1,nT_h(iparm) + write (iout,'(e20.10,$)') hfin(t,ib) + if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib) + enddo + write (iout,'(i5)') iparm + if (histfile) write (ihist,'(i5)') iparm + endif + enddo + + endif + + if (entfile) then + if (nslice.eq.1) then + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent' + else + histname=prefix(:ilen(prefix))//'.ent' + endif + else + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//'par_'//licz3// + & '_slice_'//licz2//'.ent' + else + histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent' + endif + endif +#if defined(AIX) || defined(PGI) + open (ihist,file=histname,position='append') +#else + open (ihist,file=histname,access='append') +#endif + write (ihist,'(a)') "# Microcanonical entropy" + do i=0,upindE + write (ihist,'(f8.0,$)') dint(potEmin)+i + if (histE(i).gt.0.0e0) then + write (ihist,'(f15.5,$)') dlog(histE(i)) + else + write (ihist,'(f15.5,$)') 0.0d0 + endif + enddo + write (ihist,*) + close(ihist) + endif + write (iout,*) "Microcanonical entropy" + do i=0,upindE + write (iout,'(f8.0,$)') dint(potEmin)+i + if (histE(i).gt.0.0e0) then + write (iout,'(f15.5,$)') dlog(histE(i)) + else + write (iout,'(f15.5,$)') 0.0d0 + endif + write (iout,*) + enddo + if (rmsrgymap) then + if (nslice.eq.1) then + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy' + else + histname=prefix(:ilen(prefix))//'.rmsrgy' + endif + else + if (separate_parset) then + write(licz3,"(bz,i3.3)") myparm + histname=prefix(:ilen(prefix))//'_par'//licz3// + & '_slice_'//licz2//'.rmsrgy' + else + histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy' + endif + endif +#if defined(AIX) || defined(PGI) + open (ihist,file=histname,position='append') +#else + open (ihist,file=histname,access='append') +#endif + do i=0,nbin_rms + do j=0,nbin_rgy + write(ihist,'(2f8.2,$)') + & rgymin+deltrgy*j,rmsmin+deltrms*i + do ib=1,nT_h(iparm) + if (hrmsrgy(j,i,ib).gt.0.0d0) then + write(ihist,'(e14.5,$)') + & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) + & +potEmin + else + write(ihist,'(e14.5,$)') 1.0d6 + endif + enddo + write (ihist,'(i2)') iparm + enddo + enddo + close(ihist) + endif + endif + enddo ! iparm +#ifdef MPI + call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet, + & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR) + call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1), + & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master, + & WHAM_COMM,IERROR) + call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1), + & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master, + & WHAM_COMM,IERROR) + call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1), + & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master, + & WHAM_COMM,IERROR) + if (me.eq.master) then +#endif + write (iout,'(/a)') 'Thermal characteristics of folding' + if (nslice.eq.1) then + nazwa=prefix + else + nazwa=prefix(:ilen(prefix))//"_slice_"//licz2 + endif + iln=ilen(nazwa) + if (nparmset.eq.1 .and. .not.separate_parset) then + nazwa=nazwa(:iln)//".thermal" + else if (nparmset.eq.1 .and. separate_parset) then + write(licz3,"(bz,i3.3)") myparm + nazwa=nazwa(:iln)//"_par_"//licz3//".thermal" + endif + do iparm=1,nParmSet + if (nparmset.gt.1) then + write(licz3,"(bz,i3.3)") iparm + nazwa=nazwa(:iln)//"_par_"//licz3//".thermal" + endif + open(34,file=nazwa) + if (separate_parset) then + write (iout,'(a,i3)') "Parameter set",myparm + else + write (iout,'(a,i3)') "Parameter set",iparm + endif + do i=0,NGridT + sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm) + sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ + & sumW(i,iparm) + sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) + & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2) + do j=1,nQ+2 + sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm) + sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) + & -sumQ(j,i,iparm)**2 + sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) + & -sumQ(j,i,iparm)*sumE(i,iparm) + enddo + sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* + & (startGridT+i*delta_T))+potEmin + write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T, + & sumW(i,iparm),sumE(i,iparm) + write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2) + write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm), + & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2) + write (iout,*) + write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T, + & sumW(i,iparm),sumE(i,iparm) + write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2) + write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm), + & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2) + write (34,*) + enddo + close(34) + enddo + if (histout) then + do t=0,tmax + if (hfin_ent(t).gt.0.0d0) then + liczba=t + jj = mod(liczba,nbin1) + write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta, + & hfin_ent(t) + if (histfile) write (ihist,'(f6.3,e20.10," ent")') + & dmin+(jj+0.5d0)*delta, + & hfin_ent(t) + endif + enddo + if (histfile) close(ihist) + endif + +#ifdef ZSCORE +! Write data for zscore + if (nslice.eq.1) then + zscname=prefix(:ilen(prefix))//".zsc" + else + zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc" + endif +#if defined(AIX) || defined(PGI) + open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append') +#else + open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append') +#endif + write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet + do iparm=1,nParmSet + write (izsc,'("NT=",i1)') nT_h(iparm) + do ib=1,nT_h(iparm) + write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') + & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm) + jj = min0(nR(ib,iparm),7) + write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj) + write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79) + write (izsc,'("&")') + if (nR(ib,iparm).gt.7) then + do ii=8,nR(ib,iparm),9 + jj = min0(nR(ib,iparm),ii+8) + write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj) + write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79) + write (izsc,'("&")') + enddo + endif + write (izsc,'("FI=",$)') + jj=min0(nR(ib,iparm),7) + write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj) + write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79) + write (izsc,'("&")') + if (nR(ib,iparm).gt.7) then + do ii=8,nR(ib,iparm),9 + jj = min0(nR(ib,iparm),ii+8) + write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj) + if (jj.eq.nR(ib,iparm)) then + write (izsc,*) + else + write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79) + write (izsc,'(t80,"&")') + endif + enddo + endif + do i=1,nR(ib,iparm) + write (izsc,'("KH=",$)') + write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ) + write (izsc,'(" Q0=",$)') + write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ) + write (izsc,*) + enddo + enddo + enddo + close(izsc) +#endif +#ifdef MPI + endif +#endif + + return + + end