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