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