small fix in .gitingnore
[unres.git] / 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 nGridT
18       parameter (NGridT=400)
19       integer MaxBinRms,MaxBinRgy
20       parameter (MaxBinRms=1000,MaxBinRgy=1000)
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       double precision sumW_p(0:nGridT,Max_Parm),
55      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
56      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
57      & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
58      & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
59      & sumEprim_p(MaxQ1,0:nGridT,Max_Parm), 
60      & sumEbis_p(0:nGridT,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,entmin_p,entmax_p
66       integer histent_p(0:2000)
67       logical lprint /.true./
68 #endif
69       double precision delta_T /1.0d0/
70       double precision rgymin,rmsmin,rgymax,rmsmax
71       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
72      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
73      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
74      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
75      & weight,econstr
76       double precision fi(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,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/,startGridT/200.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
89       integer ind_point(maxpoint),upindE,indE
90       character*16 plik
91       character*1 licz1
92       character*2 licz2
93       character*3 licz3
94       character*128 nazwa
95       integer ilen
96       external ilen
97  
98       write(licz2,'(bz,i2.2)') islice
99       nbin1 = 1.0d0/delta
100       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
101      &  i2/80(1h-)//)') islice
102       write (iout,*) "delta",delta," nbin1",nbin1
103       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
104       call flush(iout)
105       dmin=0.0d0
106       tmax=0
107       potEmin=1.0d10
108       rgymin=1.0d10
109       rmsmin=1.0d10
110       rgymax=0.0d0
111       rmsmax=0.0d0
112       do t=0,MaxN
113         htot(t)=0
114       enddo
115 #ifdef MPI
116       do i=1,scount(me1)
117 #else
118       do i=1,ntot(islice)
119 #endif
120         do j=1,nParmSet
121           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
122         enddo
123         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
124         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
125         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
126         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
127         ind_point(i)=0
128         do j=nQ,1,-1
129           ind=(q(j,i)-dmin+1.0d-8)/delta
130           if (j.eq.1) then
131             ind_point(i)=ind_point(i)+ind
132           else 
133             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
134           endif
135 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
136 c     &      ind_point(i)
137           call flush(iout)
138           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
139             write (iout,*) "Error - index exceeds range for point",i,
140      &      " q=",q(j,i)," ind",ind_point(i)
141 #ifdef MPI 
142             write (iout,*) "Processor",me1
143             call flush(iout)
144             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
145 #endif
146             stop
147           endif
148         enddo ! j
149         if (ind_point(i).gt.tmax) tmax=ind_point(i)
150         htot(ind_point(i))=htot(ind_point(i))+1
151 #ifdef DEBUG
152         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
153      &   " htot",htot(ind_point(i))
154         call flush(iout)
155 #endif
156       enddo ! i
157       call flush(iout)
158
159       nbin=nbin1**nQ-1
160       write (iout,'(a)') "Numbers of counts in Q bins"
161       do t=0,tmax
162         if (htot(t).gt.0) then
163         write (iout,'(i15,$)') t
164         liczba=t
165         do j=1,nQ
166           jj = mod(liczba,nbin1)
167           liczba=liczba/nbin1
168           write (iout,'(i5,$)') jj
169         enddo
170         write (iout,'(i8)') htot(t)
171         endif
172       enddo
173       do iparm=1,nParmSet
174       write (iout,'(a,i3)') "Number of data points for parameter set",
175      & iparm
176       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
177      & ib=1,nT_h(iparm))
178       write (iout,'(i8)') stot(islice)
179       write (iout,'(a)')
180       enddo
181       call flush(iout)
182
183 #ifdef MPI
184       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
185      &  WHAM_COMM,IERROR)
186       tmax=tmax_t
187       call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
188      &  MPI_MIN,WHAM_COMM,IERROR)
189       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
190      &  MPI_MIN,WHAM_COMM,IERROR)
191       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
192      &  MPI_MAX,WHAM_COMM,IERROR)
193       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
194      &  MPI_MIN,WHAM_COMM,IERROR)
195       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
196      &  MPI_MAX,WHAM_COMM,IERROR)
197       potEmin=potEmin_t/2
198       rgymin=rgymin_t
199       rgymax=rgymax_t
200       rmsmin=rmsmin_t
201       rmsmax=rmsmax_t
202       write (iout,*) "potEmin",potEmin
203 #endif
204       rmsmin=deltrms*dint(rmsmin/deltrms)
205       rmsmax=deltrms*dint(rmsmax/deltrms)
206       rgymin=deltrms*dint(rgymin/deltrgy)
207       rgymax=deltrms*dint(rgymax/deltrgy)
208       nbin_rms=(rmsmax-rmsmin)/deltrms
209       nbin_rgy=(rgymax-rgymin)/deltrgy
210       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
211      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
212       nFi=0
213       do i=1,nParmSet
214         do j=1,nT_h(i)
215           nFi=nFi+nR(j,i)
216         enddo
217       enddo
218       write (iout,*) "nFi",nFi
219 ! Compute the Boltzmann factor corresponing to restrain potentials in different
220 ! simulations.
221 #ifdef MPI
222       do i=1,scount(me1)
223 #else
224       do i=1,ntot(islice)
225 #endif
226 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
227         do iparm=1,nParmSet
228 #ifdef DEBUG
229           write (iout,'(2i5,21f8.2)') i,iparm,
230      &     (enetb(k,i,iparm),k=1,21)
231 #endif
232           call restore_parm(iparm)
233 #ifdef DEBUG
234           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
235      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
236      &      wtor_d,wsccor,wbond
237 #endif
238           do ib=1,nT_h(iparm)
239             if (rescale_mode.eq.1) then
240               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
241               quotl=1.0d0
242               kfacl=1.0d0
243               do l=1,5
244                 quotl1=quotl
245                 quotl=quotl*quot
246                 kfacl=kfacl*kfac
247                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
248               enddo
249 #if defined(FUNCTH)
250               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
251               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
252 #elif defined(FUNCT)
253               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
254 #else
255               ft(6)=1.0d0
256 #endif
257             else if (rescale_mode.eq.2) then
258               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
259               quotl=1.0d0
260               do l=1,5
261                 quotl=quotl*quot
262                 fT(l)=1.12692801104297249644d0/
263      &             dlog(dexp(quotl)+dexp(-quotl))
264               enddo
265 #if defined(FUNCTH)
266               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
267               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
268 #elif defined(FUNCT)
269               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
270 #else
271               ft(6)=1.0d0
272 #endif
273 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
274             else if (rescale_mode.eq.0) then
275               do l=1,6
276                 fT(l)=1.0d0
277               enddo
278             else
279               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
280      &          rescale_mode
281               call flush(iout)
282               return1
283             endif
284             evdw=enetb(1,i,iparm)
285             evdw_t=enetb(21,i,iparm)
286 #ifdef SCP14
287             evdw2_14=enetb(17,i,iparm)
288             evdw2=enetb(2,i,iparm)+evdw2_14
289 #else
290             evdw2=enetb(2,i,iparm)
291             evdw2_14=0.0d0
292 #endif
293 #ifdef SPLITELE
294             ees=enetb(3,i,iparm)
295             evdw1=enetb(16,i,iparm)
296 #else
297             ees=enetb(3,i,iparm)
298             evdw1=0.0d0
299 #endif
300             ecorr=enetb(4,i,iparm)
301             ecorr5=enetb(5,i,iparm)
302             ecorr6=enetb(6,i,iparm)
303             eel_loc=enetb(7,i,iparm)
304             eello_turn3=enetb(8,i,iparm)
305             eello_turn4=enetb(9,i,iparm)
306             eturn6=enetb(10,i,iparm)
307             ebe=enetb(11,i,iparm)
308             escloc=enetb(12,i,iparm)
309             etors=enetb(13,i,iparm)
310             etors_d=enetb(14,i,iparm)
311             ehpb=enetb(15,i,iparm)
312             estr=enetb(18,i,iparm)
313             esccor=enetb(19,i,iparm)
314             edihcnstr=enetb(20,i,iparm)
315 #ifdef DEBUG
316             write (iout,'(3i5,6f5.2,14f12.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 #endif
320
321 #ifdef SPLITELE
322             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
323      &      +wvdwpp*evdw1
324      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
325      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
326      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
327      &      +ft(2)*wturn3*eello_turn3
328      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
329      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
330      &      +wbond*estr
331 #else
332             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
333      &      +ft(1)*welec*(ees+evdw1)
334      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
335      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
336      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
337      &      +ft(2)*wturn3*eello_turn3
338      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
339      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
340      &      +wbond*estr
341 #endif
342 #ifdef DEBUG
343             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
344      &        etot,potEmin
345 #endif
346 #ifdef DEBUG
347             if (iparm.eq.1 .and. ib.eq.1) then
348             write (iout,*)"Conformation",i
349             energia(0)=etot
350             do k=1,max_ene
351               energia(k)=enetb(k,i,iparm)
352             enddo
353             call enerprint(energia(0),fT)
354             endif
355 #endif
356             do kk=1,nR(ib,iparm)
357               Econstr=0.0d0
358               do j=1,nQ
359                 dd = q(j,i)
360                 Econstr=Econstr+Kh(j,kk,ib,iparm)
361      &           *(dd-q0(j,kk,ib,iparm))**2
362               enddo
363               v(i,kk,ib,iparm)=
364      &          -beta_h(ib,iparm)*(etot-potEmin+Econstr)
365 #ifdef DEBUG
366               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
367      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
368 #endif
369             enddo ! kk
370           enddo   ! ib
371         enddo     ! iparm
372       enddo       ! i
373 ! Simple iteration to calculate free energies corresponding to all simulation
374 ! runs.
375       do iter=1,maxit 
376         
377 ! Compute new free-energy values corresponding to the righ-hand side of the 
378 ! equation and their derivatives.
379         write (iout,*) "------------------------fi"
380 #ifdef MPI
381         do t=1,scount(me1)
382 #else
383         do t=1,ntot(islice)
384 #endif
385           vmax=-1.0d+20
386           do i=1,nParmSet
387             do k=1,nT_h(i)
388               do l=1,nR(k,i)
389                 vf=v(t,l,k,i)+f(l,k,i)
390                 if (vf.gt.vmax) vmax=vf
391               enddo
392             enddo
393           enddo        
394           denom=0.0d0
395           do i=1,nParmSet
396             do k=1,nT_h(i)
397               do l=1,nR(k,i)
398                 aux=f(l,k,i)+v(t,l,k,i)-vmax
399                 if (aux.gt.-200.0d0)
400      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
401               enddo
402             enddo
403           enddo
404           entfac(t)=-dlog(denom)-vmax
405 #ifdef DEBUG
406           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
407 #endif
408         enddo
409         do iparm=1,nParmSet
410           do iib=1,nT_h(iparm)
411             do ii=1,nR(iib,iparm)
412 #ifdef MPI
413               fi_p(ii,iib,iparm)=0.0d0
414               do t=1,scount(me)
415                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
416      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
417 #ifdef DEBUG
418               write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
419      &         v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
420 #endif
421               enddo
422 #else
423               fi(ii,iib,iparm)=0.0d0
424               do t=1,ntot(islice)
425                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
426      &            +dexp(v(t,ii,iib,iparm)+entfac(t))
427               enddo
428 #endif
429             enddo ! ii
430           enddo ! iib
431         enddo ! iparm
432
433 #ifdef MPI
434 #ifdef DEBUG
435         write (iout,*) "fi before MPI_Reduce me",me,' master',master
436         do iparm=1,nParmSet
437           do ib=1,nT_h(nparmset)
438             write (iout,*) "iparm",iparm," ib",ib
439             write (iout,*) "beta=",beta_h(ib,iparm)
440             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
441           enddo
442         enddo
443 #endif
444         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
445      &   maxR*MaxT_h*nParmSet
446         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
447      &   " WHAM_COMM",WHAM_COMM
448         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
449      &   MPI_DOUBLE_PRECISION,
450      &   MPI_SUM,Master,WHAM_COMM,IERROR)
451 #ifdef DEBUG
452         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
453         do iparm=1,nParmSet
454           write (iout,*) "iparm",iparm
455           do ib=1,nT_h(iparm)
456             write (iout,*) "beta=",beta_h(ib,iparm)
457             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
458           enddo
459         enddo
460 #endif
461         if (me1.eq.Master) then
462 #endif
463         avefi=0.0d0
464         do iparm=1,nParmSet
465           do ib=1,nT_h(iparm)
466             do i=1,nR(ib,iparm)
467               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
468               avefi=avefi+fi(i,ib,iparm)
469             enddo
470           enddo
471         enddo
472         avefi=avefi/nFi
473         do iparm=1,nParmSet
474           write (iout,*) "Parameter set",iparm
475           do ib =1,nT_h(iparm)
476             write (iout,*) "beta=",beta_h(ib,iparm)
477             do i=1,nR(ib,iparm)
478               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
479             enddo
480             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
481             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
482           enddo
483         enddo
484
485 ! Compute the norm of free-energy increments.
486         finorm=0.0d0
487         do iparm=1,nParmSet
488           do ib=1,nT_h(iparm)
489             do i=1,nR(ib,iparm)
490               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
491               f(i,ib,iparm)=fi(i,ib,iparm)
492             enddo  
493           enddo
494         enddo
495
496         write (iout,*) 'Iteration',iter,' finorm',finorm
497
498 #ifdef MPI
499         endif
500         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
501      &   MPI_DOUBLE_PRECISION,Master,
502      &   WHAM_COMM,IERROR)
503         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
504      &   WHAM_COMM,IERROR)
505 #endif
506 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
507         if (finorm.lt.fimin) then
508           write (iout,*) 'Iteration converged'
509           goto 20
510         endif
511
512       enddo ! iter
513
514    20 continue
515 ! Now, put together the histograms from all simulations, in order to get the
516 ! unbiased total histogram.
517 #ifdef MPI
518       do t=0,tmax
519         hfin_ent_p(t)=0.0d0
520       enddo
521 #else
522       do t=0,tmax
523         hfin_ent(t)=0.0d0
524       enddo
525 #endif
526       write (iout,*) "--------------hist"
527 #ifdef MPI
528       do iparm=1,nParmSet
529         do i=0,nGridT
530           sumW_p(i,iparm)=0.0d0
531           sumE_p(i,iparm)=0.0d0
532           sumEbis_p(i,iparm)=0.0d0
533           sumEsq_p(i,iparm)=0.0d0
534           do j=1,nQ+2
535             sumQ_p(j,i,iparm)=0.0d0
536             sumQsq_p(j,i,iparm)=0.0d0
537             sumEQ_p(j,i,iparm)=0.0d0
538           enddo
539         enddo
540       enddo
541       upindE_p=0
542 #else
543       do iparm=1,nParmSet
544         do i=0,nGridT
545           sumW(i,iparm)=0.0d0
546           sumE(i,iparm)=0.0d0
547           sumEbis(i,iparm)=0.0d0
548           sumEsq(i,iparm)=0.0d0
549           do j=1,nQ+2
550             sumQ(j,i,iparm)=0.0d0
551             sumQsq(j,i,iparm)=0.0d0
552             sumEQ(j,i,iparm)=0.0d0
553           enddo
554         enddo
555       enddo
556       upindE=0
557 #endif
558 c 8/26/05 entropy distribution
559 #ifdef MPI
560       entmin_p=1.0d10
561       entmax_p=-1.0d10
562       do t=1,scount(me1)
563 c        ent=-dlog(entfac(t))
564         ent=entfac(t)
565         if (ent.lt.entmin_p) entmin_p=ent
566         if (ent.gt.entmax_p) entmax_p=ent
567       enddo
568       write (iout,*) "entmin",entmin_p," entmax",entmax_p
569       call flush(iout)
570       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
571      &  WHAM_COMM,IERROR)
572       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
573      &  WHAM_COMM,IERROR)
574       ientmax=entmax-entmin 
575       if (ientmax.gt.2000) ientmax=2000
576       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
577       call flush(iout)
578       do t=1,scount(me1)
579 c        ient=-dlog(entfac(t))-entmin
580         ient=entfac(t)-entmin
581         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
582       enddo
583       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
584      &  MPI_SUM,WHAM_COMM,IERROR)
585       if (me1.eq.Master) then
586         write (iout,*) "Entropy histogram"
587         do i=0,ientmax
588           write(iout,'(f15.4,i10)') entmin+i,histent(i)
589         enddo
590       endif
591 #else
592       entmin=1.0d10
593       entmax=-1.0d10
594       do t=1,ntot(islice)
595         ent=entfac(t)
596         if (ent.lt.entmin) entmin=ent
597         if (ent.gt.entmax) entmax=ent
598       enddo
599       ientmax=-dlog(entmax)-entmin
600       if (ientmax.gt.2000) ientmax=2000
601       do t=1,ntot(islice)
602         ient=entfac(t)-entmin
603         if (ient.le.2000) histent(ient)=histent(ient)+1
604       enddo
605       write (iout,*) "Entropy histogram"
606       do i=0,ientmax
607         write(iout,'(2f15.4)') entmin+i,histent(i)
608       enddo
609 #endif
610       
611 #ifdef MPI
612 c      write (iout,*) "me1",me1," scount",scount(me1)
613
614       do iparm=1,nParmSet
615
616 #ifdef MPI
617         do ib=1,nT_h(iparm)
618           do t=0,tmax
619             hfin_p(t,ib)=0.0d0
620           enddo
621         enddo
622         do i=1,maxindE
623           histE_p(i)=0.0d0
624         enddo
625 #else
626         do ib=1,nT_h(iparm)
627           do t=0,tmax
628             hfin(t,ib)=0.0d0
629           enddo
630         enddo
631         do i=1,maxindE
632           histE(i)=0.0d0
633         enddo
634 #endif
635         do ib=1,nT_h(iparm)
636           do i=0,MaxBinRms
637             do j=0,MaxBinRgy
638               hrmsrgy(j,i,ib)=0.0d0
639 #ifdef MPI
640               hrmsrgy_p(j,i,ib)=0.0d0
641 #endif
642             enddo
643           enddo
644         enddo
645
646         do t=1,scount(me1)
647 #else
648         do t=1,ntot(islice)
649 #endif
650           ind = ind_point(t)
651 #ifdef MPI
652           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
653 #else
654           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
655 #endif
656 c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
657           call restore_parm(iparm)
658           evdw=enetb(21,t,iparm)
659           evdw_t=enetb(1,t,iparm)
660 #ifdef SCP14
661           evdw2_14=enetb(17,t,iparm)
662           evdw2=enetb(2,t,iparm)+evdw2_14
663 #else
664           evdw2=enetb(2,t,iparm)
665           evdw2_14=0.0d0
666 #endif
667 #ifdef SPLITELE
668           ees=enetb(3,t,iparm)
669           evdw1=enetb(16,t,iparm)
670 #else
671           ees=enetb(3,t,iparm)
672           evdw1=0.0d0
673 #endif
674           ecorr=enetb(4,t,iparm)
675           ecorr5=enetb(5,t,iparm)
676           ecorr6=enetb(6,t,iparm)
677           eel_loc=enetb(7,t,iparm)
678           eello_turn3=enetb(8,t,iparm)
679           eello_turn4=enetb(9,t,iparm)
680           eturn6=enetb(10,t,iparm)
681           ebe=enetb(11,t,iparm)
682           escloc=enetb(12,t,iparm)
683           etors=enetb(13,t,iparm)
684           etors_d=enetb(14,t,iparm)
685           ehpb=enetb(15,t,iparm)
686           estr=enetb(18,t,iparm)
687           esccor=enetb(19,t,iparm)
688           edihcnstr=enetb(20,t,iparm)
689           edihcnstr=0.0d0
690           do k=0,nGridT
691             betaT=startGridT+k*delta_T
692             temper=betaT
693 c            fT=T0/betaT
694 c            ft=2*T0/(T0+betaT)
695             if (rescale_mode.eq.1) then
696               quot=betaT/T0
697               quotl=1.0d0
698               kfacl=1.0d0
699               do l=1,5
700                 quotl1=quotl
701                 quotl=quotl*quot
702                 kfacl=kfacl*kfac
703                 denom=kfacl-1.0d0+quotl
704                 fT(l)=kfacl/denom
705                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
706                 ftbis(l)=l*kfacl*quotl1*
707      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
708               enddo
709 #if defined(FUNCTH)
710               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
711      &                  320.0d0
712               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
713              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
714      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
715 #elif defined(FUNCT)
716               fT(6)=betaT/T0
717               ftprim(6)=1.0d0/T0
718               ftbis(6)=0.0d0
719 #else
720               fT(6)=1.0d0
721               ftprim(6)=0.0d0
722               ftbis(6)=0.0d0
723 #endif
724             else if (rescale_mode.eq.2) then
725               quot=betaT/T0
726               quotl=1.0d0
727               do l=1,5
728                 quotl1=quotl
729                 quotl=quotl*quot
730                 eplus=dexp(quotl)
731                 eminus=dexp(-quotl)
732                 logfac=1.0d0/dlog(eplus+eminus)
733                 tanhT=(eplus-eminus)/(eplus+eminus)
734                 fT(l)=1.12692801104297249644d0*logfac
735                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
736                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
737      &          2*l*quotl1/T0*logfac*
738      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
739      &          +ftprim(l)*tanhT)
740               enddo
741 #if defined(FUNCTH)
742               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
743      &                 320.0d0
744               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
745              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
746      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
747 #elif defined(FUNCT)
748               fT(6)=betaT/T0
749               ftprim(6)=1.0d0/T0
750               ftbis(6)=0.0d0
751 #else
752               fT(6)=1.0d0
753               ftprim(6)=0.0d0
754               ftbis(6)=0.0d0
755 #endif
756             else if (rescale_mode.eq.0) then
757               do l=1,5
758                 fT(l)=1.0d0
759                 ftprim(l)=0.0d0
760               enddo
761             else
762               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
763      &          rescale_mode
764               call flush(iout)
765               return1
766             endif
767 c            write (iout,*) "ftprim",ftprim
768 c            write (iout,*) "ftbis",ftbis
769             betaT=1.0d0/(1.987D-3*betaT)
770 #ifdef SPLITELE
771             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
772      &      +wvdwpp*evdw1
773      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
774      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
775      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
776      &      +ft(2)*wturn3*eello_turn3
777      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
778      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
779      &      +wbond*estr
780             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
781      &            +ftprim(1)*wtor*etors+
782      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
783      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
784      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
785      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
786      &            ftprim(1)*wsccor*esccor
787             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
788      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
789      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
790      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
791      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
792      &            ftbis(1)*wsccor*esccor
793 #else
794             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
795      &      +ft(1)*welec*(ees+evdw1)
796      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
797      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
798      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
799      &      +ft(2)*wturn3*eello_turn3
800      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
801      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
802      &      +wbond*estr
803             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
804      &           +ftprim(1)*wtor*etors+
805      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
806      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
807      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
808      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
809      &            ftprim(1)*wsccor*esccor
810             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
811      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
812      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
813      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
814      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
815      &            ftprim(1)*wsccor*esccor
816 #endif
817             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
818 #ifdef DEBUG
819             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
820      &        " etot",etot," entfac",entfac(t),
821      &        " weight",weight," ebis",ebis
822 #endif
823             etot=etot-temper*eprim
824 #ifdef MPI
825             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
826             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
827             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
828             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
829             do j=1,nQ+2
830               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
831               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
832               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
833      &         +etot*q(j,t)*weight
834             enddo
835 #else
836             sumW(k,iparm)=sumW(k,iparm)+weight
837             sumE(k,iparm)=sumE(k,iparm)+etot*weight
838             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
839             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
840             do j=1,nQ+2
841               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
842               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
843               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
844      &         +etot*q(j,t)*weight
845             enddo
846 #endif
847           enddo
848           indE = aint(potE(t,iparm)-aint(potEmin))
849           if (indE.ge.0 .and. indE.le.maxinde) then
850             if (indE.gt.upindE_p) upindE_p=indE
851             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
852           endif
853 #ifdef MPI
854           do ib=1,nT_h(iparm)
855             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
856             hfin_p(ind,ib)=hfin_p(ind,ib)+
857      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
858              if (rmsrgymap) then
859                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) 
860                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
861                hrmsrgy_p(indrgy,indrms,ib)=
862      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
863              endif
864           enddo
865 #else
866           do ib=1,nT_h(iparm)
867             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
868             hfin(ind,ib)=hfin(ind,ib)+
869      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
870              if (rmsrgymap) then
871                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
872                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
873                hrmsrgy(indrgy,indrms,ib)=
874      &           hrmsrgy(indrgy,indrms,ib)+expfac
875              endif
876           enddo
877 #endif
878         enddo ! t
879         do ib=1,nT_h(iparm)
880           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
881      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
882           if (rmsrgymap) then
883           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
884      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
885      &       WHAM_COMM,IERROR)
886           endif
887         enddo
888         call MPI_Reduce(upindE_p,upindE,1,
889      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
890         call MPI_Reduce(histE_p(0),histE(0),maxindE,
891      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
892
893         if (me1.eq.master) then
894
895         if (histout) then
896
897         write (iout,'(6x,$)')
898         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
899      &   ib=1,nT_h(iparm))
900         write (iout,*)
901
902         write (iout,'(/a)') 'Final histograms'
903         if (histfile) then
904           if (nslice.eq.1) then
905             if (separate_parset) then
906               write(licz3,"(bz,i3.3)") myparm
907               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
908             else
909               histname=prefix(:ilen(prefix))//'.hist'
910             endif
911           else
912             if (separate_parset) then
913               write(licz3,"(bz,i3.3)") myparm
914               histname=prefix(:ilen(prefix))//'_par'//licz3//
915      &         '_slice_'//licz2//'.hist'
916             else
917               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
918             endif
919           endif
920 #if defined(AIX) || defined(PGI)
921           open (ihist,file=histname,position='append')
922 #else
923           open (ihist,file=histname,access='append')
924 #endif
925         endif
926
927         do t=0,tmax
928           liczba=t
929           sumH=0.0d0
930           do ib=1,nT_h(iparm)
931             sumH=sumH+hfin(t,ib)
932           enddo
933           if (sumH.gt.0.0d0) then
934             do j=1,nQ
935               jj = mod(liczba,nbin1)
936               liczba=liczba/nbin1
937               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
938               if (histfile) 
939      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
940             enddo
941             do ib=1,nT_h(iparm)
942               write (iout,'(e20.10,$)') hfin(t,ib)
943               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
944             enddo
945             write (iout,'(i5)') iparm
946             if (histfile) write (ihist,'(i5)') iparm
947           endif
948         enddo
949
950         endif
951
952         if (entfile) then
953           if (nslice.eq.1) then
954             if (separate_parset) then
955               write(licz3,"(bz,i3.3)") myparm
956               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
957             else
958               histname=prefix(:ilen(prefix))//'.ent'
959             endif
960           else
961             if (separate_parset) then
962               write(licz3,"(bz,i3.3)") myparm
963               histname=prefix(:ilen(prefix))//'par_'//licz3//
964      &           '_slice_'//licz2//'.ent'
965             else
966               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
967             endif
968           endif
969 #if defined(AIX) || defined(PGI)
970           open (ihist,file=histname,position='append')
971 #else
972           open (ihist,file=histname,access='append')
973 #endif
974           write (ihist,'(a)') "# Microcanonical entropy"
975           do i=0,upindE
976             write (ihist,'(f8.0,$)') dint(potEmin)+i
977             if (histE(i).gt.0.0e0) then
978               write (ihist,'(f15.5,$)') dlog(histE(i))
979             else
980               write (ihist,'(f15.5,$)') 0.0d0
981             endif
982           enddo
983           write (ihist,*)
984           close(ihist)
985         endif
986         write (iout,*) "Microcanonical entropy"
987         do i=0,upindE
988           write (iout,'(f8.0,$)') dint(potEmin)+i
989           if (histE(i).gt.0.0e0) then
990             write (iout,'(f15.5,$)') dlog(histE(i))
991           else
992             write (iout,'(f15.5,$)') 0.0d0
993           endif
994           write (iout,*)
995         enddo
996         if (rmsrgymap) then
997           if (nslice.eq.1) then
998             if (separate_parset) then
999               write(licz3,"(bz,i3.3)") myparm
1000               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1001             else
1002               histname=prefix(:ilen(prefix))//'.rmsrgy'
1003             endif
1004           else
1005             if (separate_parset) then
1006               write(licz3,"(bz,i3.3)") myparm
1007               histname=prefix(:ilen(prefix))//'_par'//licz3//
1008      &         '_slice_'//licz2//'.rmsrgy'
1009             else
1010              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1011             endif
1012           endif
1013 #if defined(AIX) || defined(PGI)
1014           open (ihist,file=histname,position='append')
1015 #else
1016           open (ihist,file=histname,access='append')
1017 #endif
1018           do i=0,nbin_rms
1019             do j=0,nbin_rgy
1020               write(ihist,'(2f8.2,$)') 
1021      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1022               do ib=1,nT_h(iparm)
1023                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1024                   write(ihist,'(e14.5,$)') 
1025      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1026      &              +potEmin
1027                 else
1028                   write(ihist,'(e14.5,$)') 1.0d6
1029                 endif
1030               enddo
1031               write (ihist,'(i2)') iparm
1032             enddo
1033           enddo
1034           close(ihist)
1035         endif
1036         endif
1037       enddo ! iparm
1038 #ifdef MPI
1039       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1040      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1041       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1042      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1043       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1044      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1045       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1046      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1047       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1048      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1049       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1050      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1051      &   WHAM_COMM,IERROR)
1052       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1053      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1054      &   WHAM_COMM,IERROR)
1055       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1056      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1057      &   WHAM_COMM,IERROR)
1058       if (me.eq.master) then
1059 #endif
1060       write (iout,'(/a)') 'Thermal characteristics of folding'
1061       if (nslice.eq.1) then
1062         nazwa=prefix
1063       else
1064         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1065       endif
1066       iln=ilen(nazwa)
1067       if (nparmset.eq.1 .and. .not.separate_parset) then
1068         nazwa=nazwa(:iln)//".thermal"
1069       else if (nparmset.eq.1 .and. separate_parset) then
1070         write(licz3,"(bz,i3.3)") myparm
1071         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1072       endif
1073       do iparm=1,nParmSet
1074       if (nparmset.gt.1) then
1075         write(licz3,"(bz,i3.3)") iparm
1076         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1077       endif
1078       open(34,file=nazwa)
1079       if (separate_parset) then
1080         write (iout,'(a,i3)') "Parameter set",myparm
1081       else
1082         write (iout,'(a,i3)') "Parameter set",iparm
1083       endif
1084       do i=0,NGridT
1085         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1086         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1087      &    sumW(i,iparm)
1088         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1089      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1090         do j=1,nQ+2
1091           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1092           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1093      &     -sumQ(j,i,iparm)**2
1094           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1095      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1096         enddo
1097         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1098      &   (startGridT+i*delta_T))+potEmin
1099         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1100      &   sumW(i,iparm),sumE(i,iparm)
1101         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1102         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1103      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1104         write (iout,*) 
1105         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1106      &   sumW(i,iparm),sumE(i,iparm)
1107         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1108         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1109      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1110         write (34,*) 
1111       enddo
1112       close(34)
1113       enddo
1114       if (histout) then
1115       do t=0,tmax
1116         if (hfin_ent(t).gt.0.0d0) then
1117           liczba=t
1118           jj = mod(liczba,nbin1)
1119           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1120      &     hfin_ent(t)
1121           if (histfile) write (ihist,'(f6.3,e20.10," ent")') 
1122      &      dmin+(jj+0.5d0)*delta,
1123      &     hfin_ent(t)
1124         endif
1125       enddo
1126       if (histfile) close(ihist)
1127       endif
1128
1129 #ifdef ZSCORE
1130 ! Write data for zscore
1131       if (nslice.eq.1) then
1132         zscname=prefix(:ilen(prefix))//".zsc"
1133       else
1134         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1135       endif
1136 #if defined(AIX) || defined(PGI)
1137       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1138 #else
1139       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1140 #endif
1141       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1142       do iparm=1,nParmSet
1143         write (izsc,'("NT=",i1)') nT_h(iparm)
1144       do ib=1,nT_h(iparm)
1145         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') 
1146      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1147         jj = min0(nR(ib,iparm),7)
1148         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1149         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1150         write (izsc,'("&")')
1151         if (nR(ib,iparm).gt.7) then
1152           do ii=8,nR(ib,iparm),9
1153             jj = min0(nR(ib,iparm),ii+8)
1154             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj) 
1155             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1156             write (izsc,'("&")')
1157           enddo
1158         endif
1159         write (izsc,'("FI=",$)')
1160         jj=min0(nR(ib,iparm),7)
1161         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1162         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1163         write (izsc,'("&")')
1164         if (nR(ib,iparm).gt.7) then
1165           do ii=8,nR(ib,iparm),9
1166             jj = min0(nR(ib,iparm),ii+8)
1167             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj) 
1168             if (jj.eq.nR(ib,iparm)) then
1169               write (izsc,*) 
1170             else
1171               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1172               write (izsc,'(t80,"&")')
1173             endif
1174           enddo
1175         endif
1176         do i=1,nR(ib,iparm)
1177           write (izsc,'("KH=",$)') 
1178           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1179           write (izsc,'(" Q0=",$)')
1180           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1181           write (izsc,*)
1182         enddo
1183       enddo
1184       enddo
1185       close(izsc)
1186 #endif
1187 #ifdef MPI
1188       endif
1189 #endif
1190
1191       return
1192
1193       end