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