Merge branch 'master' of mmka:unres into devel
[unres.git] / source / wham / src / wham_calc1.F.safe
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       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 #define DEBUG
819 #ifdef DEBUG
820             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
821      &        " etot",etot," entfac",entfac(t),
822      &        " weight",weight," ebis",ebis
823 #endif
824 #undef DEBUG
825             etot=etot-temper*eprim
826 #ifdef MPI
827             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
828             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
829             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
830             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
831             do j=1,nQ+2
832               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
833               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
834               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
835      &         +etot*q(j,t)*weight
836             enddo
837 #else
838             sumW(k,iparm)=sumW(k,iparm)+weight
839             sumE(k,iparm)=sumE(k,iparm)+etot*weight
840             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
841             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
842             do j=1,nQ+2
843               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
844               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
845               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
846      &         +etot*q(j,t)*weight
847             enddo
848 #endif
849           enddo
850           indE = aint(potE(t,iparm)-aint(potEmin))
851           if (indE.ge.0 .and. indE.le.maxinde) then
852             if (indE.gt.upindE_p) upindE_p=indE
853             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
854           endif
855 #ifdef MPI
856           do ib=1,nT_h(iparm)
857             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
858             hfin_p(ind,ib)=hfin_p(ind,ib)+
859      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
860              if (rmsrgymap) then
861                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) 
862                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
863                hrmsrgy_p(indrgy,indrms,ib)=
864      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
865              endif
866           enddo
867 #else
868           do ib=1,nT_h(iparm)
869             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
870             hfin(ind,ib)=hfin(ind,ib)+
871      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
872              if (rmsrgymap) then
873                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
874                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
875                hrmsrgy(indrgy,indrms,ib)=
876      &           hrmsrgy(indrgy,indrms,ib)+expfac
877              endif
878           enddo
879 #endif
880         enddo ! t
881         do ib=1,nT_h(iparm)
882           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
883      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
884           if (rmsrgymap) then
885           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
886      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
887      &       WHAM_COMM,IERROR)
888           endif
889         enddo
890         call MPI_Reduce(upindE_p,upindE,1,
891      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
892         call MPI_Reduce(histE_p(0),histE(0),maxindE,
893      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
894
895         if (me1.eq.master) then
896
897         if (histout) then
898
899         write (iout,'(6x,$)')
900         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
901      &   ib=1,nT_h(iparm))
902         write (iout,*)
903
904         write (iout,'(/a)') 'Final histograms'
905         if (histfile) then
906           if (nslice.eq.1) then
907             if (separate_parset) then
908               write(licz3,"(bz,i3.3)") myparm
909               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
910             else
911               histname=prefix(:ilen(prefix))//'.hist'
912             endif
913           else
914             if (separate_parset) then
915               write(licz3,"(bz,i3.3)") myparm
916               histname=prefix(:ilen(prefix))//'_par'//licz3//
917      &         '_slice_'//licz2//'.hist'
918             else
919               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
920             endif
921           endif
922 #if defined(AIX) || defined(PGI)
923           open (ihist,file=histname,position='append')
924 #else
925           open (ihist,file=histname,access='append')
926 #endif
927         endif
928
929         do t=0,tmax
930           liczba=t
931           sumH=0.0d0
932           do ib=1,nT_h(iparm)
933             sumH=sumH+hfin(t,ib)
934           enddo
935           if (sumH.gt.0.0d0) then
936             do j=1,nQ
937               jj = mod(liczba,nbin1)
938               liczba=liczba/nbin1
939               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
940               if (histfile) 
941      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
942             enddo
943             do ib=1,nT_h(iparm)
944               write (iout,'(e20.10,$)') hfin(t,ib)
945               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
946             enddo
947             write (iout,'(i5)') iparm
948             if (histfile) write (ihist,'(i5)') iparm
949           endif
950         enddo
951
952         endif
953
954         if (entfile) then
955           if (nslice.eq.1) then
956             if (separate_parset) then
957               write(licz3,"(bz,i3.3)") myparm
958               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
959             else
960               histname=prefix(:ilen(prefix))//'.ent'
961             endif
962           else
963             if (separate_parset) then
964               write(licz3,"(bz,i3.3)") myparm
965               histname=prefix(:ilen(prefix))//'par_'//licz3//
966      &           '_slice_'//licz2//'.ent'
967             else
968               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
969             endif
970           endif
971 #if defined(AIX) || defined(PGI)
972           open (ihist,file=histname,position='append')
973 #else
974           open (ihist,file=histname,access='append')
975 #endif
976           write (ihist,'(a)') "# Microcanonical entropy"
977           do i=0,upindE
978             write (ihist,'(f8.0,$)') dint(potEmin)+i
979             if (histE(i).gt.0.0e0) then
980               write (ihist,'(f15.5,$)') dlog(histE(i))
981             else
982               write (ihist,'(f15.5,$)') 0.0d0
983             endif
984           enddo
985           write (ihist,*)
986           close(ihist)
987         endif
988         write (iout,*) "Microcanonical entropy"
989         do i=0,upindE
990           write (iout,'(f8.0,$)') dint(potEmin)+i
991           if (histE(i).gt.0.0e0) then
992             write (iout,'(f15.5,$)') dlog(histE(i))
993           else
994             write (iout,'(f15.5,$)') 0.0d0
995           endif
996           write (iout,*)
997         enddo
998         if (rmsrgymap) then
999           if (nslice.eq.1) then
1000             if (separate_parset) then
1001               write(licz3,"(bz,i3.3)") myparm
1002               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1003             else
1004               histname=prefix(:ilen(prefix))//'.rmsrgy'
1005             endif
1006           else
1007             if (separate_parset) then
1008               write(licz3,"(bz,i3.3)") myparm
1009               histname=prefix(:ilen(prefix))//'_par'//licz3//
1010      &         '_slice_'//licz2//'.rmsrgy'
1011             else
1012              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1013             endif
1014           endif
1015 #if defined(AIX) || defined(PGI)
1016           open (ihist,file=histname,position='append')
1017 #else
1018           open (ihist,file=histname,access='append')
1019 #endif
1020           do i=0,nbin_rms
1021             do j=0,nbin_rgy
1022               write(ihist,'(2f8.2,$)') 
1023      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1024               do ib=1,nT_h(iparm)
1025                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1026                   write(ihist,'(e14.5,$)') 
1027      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1028      &              +potEmin
1029                 else
1030                   write(ihist,'(e14.5,$)') 1.0d6
1031                 endif
1032               enddo
1033               write (ihist,'(i2)') iparm
1034             enddo
1035           enddo
1036           close(ihist)
1037         endif
1038         endif
1039       enddo ! iparm
1040 #ifdef MPI
1041       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1042      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1043       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1044      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1045       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1046      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1047       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1048      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1049       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1050      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1051       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1052      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1053      &   WHAM_COMM,IERROR)
1054       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1055      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1056      &   WHAM_COMM,IERROR)
1057       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1058      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1059      &   WHAM_COMM,IERROR)
1060       if (me.eq.master) then
1061 #endif
1062       write (iout,'(/a)') 'Thermal characteristics of folding'
1063       if (nslice.eq.1) then
1064         nazwa=prefix
1065       else
1066         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1067       endif
1068       iln=ilen(nazwa)
1069       if (nparmset.eq.1 .and. .not.separate_parset) then
1070         nazwa=nazwa(:iln)//".thermal"
1071       else if (nparmset.eq.1 .and. separate_parset) then
1072         write(licz3,"(bz,i3.3)") myparm
1073         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1074       endif
1075       do iparm=1,nParmSet
1076       if (nparmset.gt.1) then
1077         write(licz3,"(bz,i3.3)") iparm
1078         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1079       endif
1080       open(34,file=nazwa)
1081       if (separate_parset) then
1082         write (iout,'(a,i3)') "Parameter set",myparm
1083       else
1084         write (iout,'(a,i3)') "Parameter set",iparm
1085       endif
1086       do i=0,NGridT
1087         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1088         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1089      &    sumW(i,iparm)
1090         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1091      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1092         do j=1,nQ+2
1093           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1094           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1095      &     -sumQ(j,i,iparm)**2
1096           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1097      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1098         enddo
1099         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1100      &   (startGridT+i*delta_T))+potEmin
1101         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1102      &   sumW(i,iparm),sumE(i,iparm)
1103         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1104         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1105      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1106         write (iout,*) 
1107         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1108      &   sumW(i,iparm),sumE(i,iparm)
1109         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1110         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1111      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1112         write (34,*) 
1113       enddo
1114       close(34)
1115       enddo
1116       if (histout) then
1117       do t=0,tmax
1118         if (hfin_ent(t).gt.0.0d0) then
1119           liczba=t
1120           jj = mod(liczba,nbin1)
1121           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1122      &     hfin_ent(t)
1123           if (histfile) write (ihist,'(f6.3,e20.10," ent")') 
1124      &      dmin+(jj+0.5d0)*delta,
1125      &     hfin_ent(t)
1126         endif
1127       enddo
1128       if (histfile) close(ihist)
1129       endif
1130
1131 #ifdef ZSCORE
1132 ! Write data for zscore
1133       if (nslice.eq.1) then
1134         zscname=prefix(:ilen(prefix))//".zsc"
1135       else
1136         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1137       endif
1138 #if defined(AIX) || defined(PGI)
1139       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1140 #else
1141       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1142 #endif
1143       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1144       do iparm=1,nParmSet
1145         write (izsc,'("NT=",i1)') nT_h(iparm)
1146       do ib=1,nT_h(iparm)
1147         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') 
1148      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1149         jj = min0(nR(ib,iparm),7)
1150         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1151         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1152         write (izsc,'("&")')
1153         if (nR(ib,iparm).gt.7) then
1154           do ii=8,nR(ib,iparm),9
1155             jj = min0(nR(ib,iparm),ii+8)
1156             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj) 
1157             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1158             write (izsc,'("&")')
1159           enddo
1160         endif
1161         write (izsc,'("FI=",$)')
1162         jj=min0(nR(ib,iparm),7)
1163         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1164         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1165         write (izsc,'("&")')
1166         if (nR(ib,iparm).gt.7) then
1167           do ii=8,nR(ib,iparm),9
1168             jj = min0(nR(ib,iparm),ii+8)
1169             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj) 
1170             if (jj.eq.nR(ib,iparm)) then
1171               write (izsc,*) 
1172             else
1173               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1174               write (izsc,'(t80,"&")')
1175             endif
1176           enddo
1177         endif
1178         do i=1,nR(ib,iparm)
1179           write (izsc,'("KH=",$)') 
1180           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1181           write (izsc,'(" Q0=",$)')
1182           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1183           write (izsc,*)
1184         enddo
1185       enddo
1186       enddo
1187       close(izsc)
1188 #endif
1189 #ifdef MPI
1190       endif
1191 #endif
1192
1193       return
1194
1195       end