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)
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
12 ! 2/2/05 Parallel version
15 include "DIMENSIONS.ZSCOPT"
16 include "DIMENSIONS.FREE"
17 integer MaxBinRms,MaxBinRgy
18 parameter (MaxBinRms=100,MaxBinRgy=100)
20 c parameter (MaxHdim=200000)
21 parameter (MaxHdim=200)
23 parameter (maxinde=200)
27 integer ierror,errcode,status(MPI_STATUS_SIZE)
29 include "COMMON.CONTROL"
30 include "COMMON.IOUNITS"
32 include "COMMON.ENERGIES"
33 include "COMMON.FFIELD"
34 include "COMMON.SBRIDGE"
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)
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)
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./
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,
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
86 integer ind_point(maxpoint),upindE,indE
95 write(licz2,'(bz,i2.2)') islice
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
106 potEmin_all(j,i)=1.0d10
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)
127 ind=(q(j,i)-dmin+1.0d-8)/delta
129 ind_point(i)=ind_point(i)+ind
131 ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
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)
137 write (iout,*) "Processor",me1
139 call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
144 if (ind_point(i).gt.tmax) tmax=ind_point(i)
145 htot(ind_point(i))=htot(ind_point(i))+1
147 write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
148 & " htot",htot(ind_point(i))
155 write (iout,'(a)') "Numbers of counts in Q bins"
157 if (htot(t).gt.0) then
158 write (iout,'(i15,$)') t
161 jj = mod(liczba,nbin1)
163 write (iout,'(i5,$)') jj
165 write (iout,'(i8)') htot(t)
169 write (iout,'(a,i3)') "Number of data points for parameter set",
171 write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
173 write (iout,'(i8)') stot(islice)
179 call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
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)
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
209 write (iout,*) "nFi",nFi
210 ! Compute the Boltzmann factor corresponing to restrain potentials in different
217 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
220 write (iout,'(2i5,21f8.2)') i,iparm,
221 & (enetb(k,i,iparm),k=1,21)
223 call restore_parm(iparm)
225 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
226 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
227 & wtor_d,wsccor,wbond
230 if (rescale_mode.eq.1) then
231 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
238 fT(l)=kfacl/(kfacl-1.0d0+quotl)
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
244 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
248 else if (rescale_mode.eq.2) then
249 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
253 fT(l)=1.12692801104297249644d0/
254 & dlog(dexp(quotl)+dexp(-quotl))
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
260 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
264 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
265 else if (rescale_mode.eq.0) then
270 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
275 evdw=enetb(1,i,iparm)
276 evdw_t=enetb(21,i,iparm)
278 evdw2_14=enetb(17,i,iparm)
279 evdw2=enetb(2,i,iparm)+evdw2_14
281 evdw2=enetb(2,i,iparm)
286 evdw1=enetb(16,i,iparm)
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)
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
313 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
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
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
334 write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
338 if (iparm.eq.1 .and. ib.eq.1) then
339 write (iout,*)"Conformation",i
342 energia(k)=enetb(k,i,iparm)
344 call enerprint(energia(0),fT)
351 Econstr=Econstr+Kh(j,kk,ib,iparm)
352 & *(dd-q0(j,kk,ib,iparm))**2
355 & -beta_h(ib,iparm)*(etot+Econstr)
357 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
358 & etot,v(i,kk,ib,iparm)
364 ! Simple iteration to calculate free energies corresponding to all simulation
368 ! Compute new free-energy values corresponding to the righ-hand side of the
369 ! equation and their derivatives.
370 write (iout,*) "------------------------fi"
380 vf=v(t,l,k,i)+f(l,k,i)
381 if (vf.gt.vmax) vmax=vf
389 aux=f(l,k,i)+v(t,l,k,i)-vmax
391 & denom=denom+snk(l,k,i,islice)*dexp(aux)
395 entfac(t)=-dlog(denom)-vmax
397 write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
402 do ii=1,nR(iib,iparm)
404 fi_p(ii,iib,iparm)=0.0d0
406 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
407 & +dexp(v(t,ii,iib,iparm)+entfac(t))
409 write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
410 & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
414 fi(ii,iib,iparm)=0.0d0
416 fi(ii,iib,iparm)=fi(ii,iib,iparm)
417 & +dexp(v(t,ii,iib,iparm)+entfac(t))
426 write (iout,*) "fi before MPI_Reduce me",me,' master',master
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))
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)
443 write (iout,*) "fi after MPI_Reduce nparmset",nparmset
445 write (iout,*) "iparm",iparm
447 write (iout,*) "beta=",beta_h(ib,iparm)
448 write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
452 if (me1.eq.Master) then
458 fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
459 avefi=avefi+fi(i,ib,iparm)
465 write (iout,*) "Parameter set",iparm
467 write (iout,*) "beta=",beta_h(ib,iparm)
469 fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
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))
476 ! Compute the norm of free-energy increments.
481 finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
482 f(i,ib,iparm)=fi(i,ib,iparm)
487 write (iout,*) 'Iteration',iter,' finorm',finorm
491 call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
492 & MPI_DOUBLE_PRECISION,Master,
494 call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
497 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
498 if (finorm.lt.fimin) then
499 write (iout,*) 'Iteration converged'
506 ! Now, put together the histograms from all simulations, in order to get the
507 ! unbiased total histogram.
509 C Determine the minimum free energies
515 c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
518 write (iout,'(2i5,21f8.2)') i,iparm,
519 & (enetb(k,i,iparm),k=1,21)
521 call restore_parm(iparm)
523 write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
524 & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
525 & wtor_d,wsccor,wbond
528 if (rescale_mode.eq.1) then
529 quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
536 fT(l)=kfacl/(kfacl-1.0d0+quotl)
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
542 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
546 else if (rescale_mode.eq.2) then
547 quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
551 fT(l)=1.12692801104297249644d0/
552 & dlog(dexp(quotl)+dexp(-quotl))
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
558 ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
562 c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
563 else if (rescale_mode.eq.0) then
568 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
573 evdw=enetb(1,i,iparm)
574 evdw_t=enetb(21,i,iparm)
576 evdw2_14=enetb(17,i,iparm)
577 evdw2=enetb(2,i,iparm)+evdw2_14
579 evdw2=enetb(2,i,iparm)
584 evdw1=enetb(16,i,iparm)
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)
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
611 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
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
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
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)
641 write (iout,*) "The potEmin array before reduction"
643 write (iout,*) "Parameter set",i
645 write (iout,*) j,PotEmin_all(j,i)
648 write (iout,*) "potEmin_min",potEmin_min
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)
656 potEmin_all(j,i)=potEmin_t_all(j,i)
660 potEmin_min=potEmin_all(1,1)
663 if (potEmin_all(j,i).lt.potEmin_min)
664 & potEmin_min=potEmin_all(j,i)
668 write (iout,*) "The potEmin array"
670 write (iout,*) "Parameter set",i
672 write (iout,*) j,PotEmin_all(j,i)
675 write (iout,*) "potEmin_min",potEmin_min
688 write (iout,*) "--------------hist"
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
697 sumQ_p(j,i,iparm)=0.0d0
698 sumQsq_p(j,i,iparm)=0.0d0
699 sumEQ_p(j,i,iparm)=0.0d0
709 sumEbis(i,iparm)=0.0d0
710 sumEsq(i,iparm)=0.0d0
712 sumQ(j,i,iparm)=0.0d0
713 sumQsq(j,i,iparm)=0.0d0
714 sumEQ(j,i,iparm)=0.0d0
720 c 8/26/05 entropy distribution
725 c ent=-dlog(entfac(t))
727 if (ent.lt.entmin_p) entmin_p=ent
728 if (ent.gt.entmax_p) entmax_p=ent
730 write (iout,*) "entmin",entmin_p," entmax",entmax_p
732 call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
734 call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
736 ientmax=entmax-entmin
737 if (ientmax.gt.2000) ientmax=2000
738 write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
741 c ient=-dlog(entfac(t))-entmin
742 ient=entfac(t)-entmin
743 if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
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"
750 write(iout,'(f15.4,i10)') entmin+i,histent(i)
758 if (ent.lt.entmin) entmin=ent
759 if (ent.gt.entmax) entmax=ent
761 ientmax=-dlog(entmax)-entmin
762 if (ientmax.gt.2000) ientmax=2000
764 ient=entfac(t)-entmin
765 if (ient.le.2000) histent(ient)=histent(ient)+1
767 write (iout,*) "Entropy histogram"
769 write(iout,'(2f15.4)') entmin+i,histent(i)
774 c write (iout,*) "me1",me1," scount",scount(me1)
800 hrmsrgy(j,i,ib)=0.0d0
802 hrmsrgy_p(j,i,ib)=0.0d0
814 hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
816 hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
818 call restore_parm(iparm)
819 evdw=enetb(21,t,iparm)
820 evdw_t=enetb(1,t,iparm)
822 evdw2_14=enetb(17,t,iparm)
823 evdw2=enetb(2,t,iparm)+evdw2_14
825 evdw2=enetb(2,t,iparm)
830 evdw1=enetb(16,t,iparm)
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)
851 betaT=startGridT+k*delta_T
855 if (rescale_mode.eq.1) then
863 denom=kfacl-1.0d0+quotl
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)
870 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.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)
884 else if (rescale_mode.eq.2) then
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)
902 ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.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)
916 else if (rescale_mode.eq.0) then
922 write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
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
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
948 c write (iout,*) ib," PotEmin",potEmin
950 etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
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
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
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
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
996 weight=dexp(-betaT*(etot-potEmin)+entfac(t))
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
1004 etot=etot-temper*eprim
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
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
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
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
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))
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))
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
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))
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
1063 if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1064 & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
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,
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)
1076 if (me1.eq.master) then
1080 write (iout,'(6x,$)')
1081 write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1085 write (iout,'(/a)') 'Final histograms'
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'
1092 histname=prefix(:ilen(prefix))//'.hist'
1095 if (separate_parset) then
1096 write(licz3,"(bz,i3.3)") myparm
1097 histname=prefix(:ilen(prefix))//'_par'//licz3//
1098 & '_slice_'//licz2//'.hist'
1100 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1103 #if defined(AIX) || defined(PGI)
1104 open (ihist,file=histname,position='append')
1106 open (ihist,file=histname,access='append')
1114 sumH=sumH+hfin(t,ib)
1116 if (sumH.gt.0.0d0) then
1118 jj = mod(liczba,nbin1)
1120 write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1122 & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1125 write (iout,'(e20.10,$)') hfin(t,ib)
1126 if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1128 write (iout,'(i5)') iparm
1129 if (histfile) write (ihist,'(i5)') iparm
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'
1141 histname=prefix(:ilen(prefix))//'.ent'
1144 if (separate_parset) then
1145 write(licz3,"(bz,i3.3)") myparm
1146 histname=prefix(:ilen(prefix))//'par_'//licz3//
1147 & '_slice_'//licz2//'.ent'
1149 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1152 #if defined(AIX) || defined(PGI)
1153 open (ihist,file=histname,position='append')
1155 open (ihist,file=histname,access='append')
1157 write (ihist,'(a)') "# Microcanonical entropy"
1159 write (ihist,'(f8.0,$)') dint(potEmin)+i
1160 if (histE(i).gt.0.0e0) then
1161 write (ihist,'(f15.5,$)') dlog(histE(i))
1163 write (ihist,'(f15.5,$)') 0.0d0
1169 write (iout,*) "Microcanonical entropy"
1171 write (iout,'(f8.0,$)') dint(potEmin)+i
1172 if (histE(i).gt.0.0e0) then
1173 write (iout,'(f15.5,$)') dlog(histE(i))
1175 write (iout,'(f15.5,$)') 0.0d0
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'
1185 histname=prefix(:ilen(prefix))//'.rmsrgy'
1188 if (separate_parset) then
1189 write(licz3,"(bz,i3.3)") myparm
1190 histname=prefix(:ilen(prefix))//'_par'//licz3//
1191 & '_slice_'//licz2//'.rmsrgy'
1193 histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1196 #if defined(AIX) || defined(PGI)
1197 open (ihist,file=histname,position='append')
1199 open (ihist,file=histname,access='append')
1203 write(ihist,'(2f8.2,$)')
1204 & rgymin+deltrgy*j,rmsmin+deltrms*i
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)
1211 write(ihist,'(e14.5,$)') 1.0d6
1214 write (ihist,'(i2)') iparm
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,
1235 call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1236 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1238 call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1239 & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1241 if (me.eq.master) then
1243 write (iout,'(/a)') 'Thermal characteristics of folding'
1244 if (nslice.eq.1) then
1247 nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
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"
1257 if (nparmset.gt.1) then
1258 write(licz3,"(bz,i3.3)") iparm
1259 nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1262 if (separate_parset) then
1263 write (iout,'(a,i3)') "Parameter set",myparm
1265 write (iout,'(a,i3)') "Parameter set",iparm
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)
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)
1282 sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1283 sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(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)
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)
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)
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)
1313 if (hfin_ent(t).gt.0.0d0) then
1315 jj = mod(liczba,nbin1)
1316 write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1318 if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1319 & dmin+(jj+0.5d0)*delta,
1323 if (histfile) close(ihist)
1327 ! Write data for zscore
1328 if (nslice.eq.1) then
1329 zscname=prefix(:ilen(prefix))//".zsc"
1331 zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1333 #if defined(AIX) || defined(PGI)
1334 open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1336 open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1338 write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1340 write (izsc,'("NT=",i1)') 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,'("&")')
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
1368 write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1369 write (izsc,'(t80,"&")')
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)