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