1 subroutine proc_data(nvarr,x,*)
4 include "DIMENSIONS.ZSCOPT"
7 integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
10 include "COMMON.CONTROL"
11 include "COMMON.ENERGIES"
12 include "COMMON.VMCPAR"
13 include "COMMON.OPTIM"
14 include "COMMON.WEIGHTS"
15 include "COMMON.PROTNAME"
16 include "COMMON.IOUNITS"
17 include "COMMON.FFIELD"
18 include "COMMON.PROTFILES"
19 include "COMMON.COMPAR"
20 include "COMMON.CHAIN"
21 include "COMMON.CLASSES"
22 include "COMMON.THERMAL"
23 include "COMMON.WEIGHTDER"
25 include "COMMON.SBRIDGE"
26 include "COMMON.INTERACT"
28 double precision fac0,t0
29 double precision fT(5),ftprim(5),ftbis(5),quot,quotl,quotl1,
30 & denom,kfac,kfacl,logfac,eplus,eminus,tanhT,entfacmax,
32 real csingle(3,maxres2)
33 double precision creff(3,maxres2),cc(3,maxres2)
34 double precision thetareff(maxres2),phireff(maxres2),
35 & xxreff(maxres2),yyreff(maxres2),zzreff(maxres2)
36 double precision przes(3),obrot(3,3),etoti
39 integer iprot,i,ind,ii,j,k,kk,l,ncl,is,ie,nsum,nng,nvarr,ng,
40 & inn,in,iref,ib,ibatch,num,itmp,imin
41 integer l1,l2,ncon_all
42 double precision rcyfr
43 integer jj,istart_conf,iend_conf,ipass_conf,iii
46 double precision ej,emaxj,sumP,sumP_t
48 double precision x(max_paropt)
49 double precision wconf(maxstr_proc)
50 double precision dx,dy,dz,dxref,dyref,dzref,
51 & dtheta,dphi,rmsthet,rmsphi,rmsside
52 double precision pinorm
55 double precision e_total_T_(maxstr_proc)
59 double precision rms,rmsave,rmsmin,rmsnat,qwolynes,qwolyness
63 write (iout,*) "fac0",fac0
64 write(iout,*) "torsional and valence angle mode",tor_mode
66 c Determine the number of variables and the initial vector of optimizable
70 if (restart_flag) then
71 call w2x(nvarr,x_orig,*11)
75 open (88,file=prefix(:ilen(prefix))//'.restin',status='unknown')
76 read(88,*,end=99,err=99) (x(i),i=1,nvarr)
80 & "Variables replaced with values from the restart file."
82 write (iout,'(i5,f15.5)') i,x(i)
88 & "Error - restart file nonexistent, incompatible or damaged."
90 call MPI_Finalize(Ierror)
93 & "Error - restart file nonexistent, incompatible or damaged."
96 write (iout,*) "PROC: after W2X nvarr",nvarr
106 c Setup weights for UNRES
135 call restore_molinfo(iprot)
137 DO IBATCH=1,natlike(iprot)+2
139 c Calculate temperature-dependent weight scaling factors
140 do ib=1,nbeta(ibatch,iprot)
141 fac=betaT(ib,ibatch,iprot)
143 if (rescale_mode.eq.0) then
149 else if (rescale_mode.eq.1) then
156 denom=kfacl-1.0d0+quotl
158 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
159 ftbis(l)=l*kfacl*quotl1*
160 & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
161 write (iout,*) "ib",ib," beta",betaT(ib,ibatch,iprot),
162 & " quot",quot," l",l,fT(l),fTprim(l),fTbis(l)
164 else if (rescale_mode.eq.2) then
171 logfac=1.0d0/dlog(eplus+eminus)
172 tanhT=(eplus-eminus)/(eplus+eminus)
173 fT(l)=1.12692801104297249644d0*logfac
174 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
175 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
176 & 2*l*quotl1/T0*logfac*
177 & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
181 write (iout,*) "Unknown RESCALE_MODE",rescale_mode
186 escal(k,ib,ibatch,iprot)=1.0d0
187 escalprim(k,ib,ibatch,iprot)=0.0d0
188 escalbis(k,ib,ibatch,iprot)=0.0d0
190 if (shield_mode.gt.0) then
191 escal(1,ib,ibatch,iprot)=ft(1)
192 escal(2,ib,ibatch,iprot)=ft(1)
193 escal(16,ib,ibatch,iprot)=ft(1)
194 escalprim(1,ib,ibatch,iprot)=ft(1)
195 escalprim(2,ib,ibatch,iprot)=ft(1)
196 escalprim(16,ib,ibatch,iprot)=ft(1)
197 escalbis(1,ib,ibatch,iprot)=ft(1)
198 escalbis(2,ib,ibatch,iprot)=ft(1)
199 escalbis(16,ib,ibatch,iprot)=ft(1)
201 escal(3,ib,ibatch,iprot)=ft(1)
202 escal(4,ib,ibatch,iprot)=ft(3)
203 escal(5,ib,ibatch,iprot)=ft(4)
204 escal(6,ib,ibatch,iprot)=ft(5)
205 escal(7,ib,ibatch,iprot)=ft(2)
206 escal(8,ib,ibatch,iprot)=ft(2)
207 escal(9,ib,ibatch,iprot)=ft(3)
208 escal(10,ib,ibatch,iprot)=ft(5)
209 escal(13,ib,ibatch,iprot)=ft(1)
210 escal(14,ib,ibatch,iprot)=ft(2)
211 escal(19,ib,ibatch,iprot)=ft(1)
212 escalprim(3,ib,ibatch,iprot)=ftprim(1)
213 escalprim(4,ib,ibatch,iprot)=ftprim(3)
214 escalprim(5,ib,ibatch,iprot)=ftprim(4)
215 escalprim(6,ib,ibatch,iprot)=ftprim(5)
216 escalprim(7,ib,ibatch,iprot)=ftprim(2)
217 escalprim(8,ib,ibatch,iprot)=ftprim(2)
218 escalprim(9,ib,ibatch,iprot)=ftprim(3)
219 escalprim(10,ib,ibatch,iprot)=ftprim(5)
220 escalprim(13,ib,ibatch,iprot)=ftprim(1)
221 escalprim(14,ib,ibatch,iprot)=ftprim(2)
222 escalprim(19,ib,ibatch,iprot)=ftprim(1)
223 escalbis(3,ib,ibatch,iprot)=ftbis(1)
224 escalbis(4,ib,ibatch,iprot)=ftbis(3)
225 escalbis(5,ib,ibatch,iprot)=ftbis(4)
226 escalbis(6,ib,ibatch,iprot)=ftbis(5)
227 escalbis(7,ib,ibatch,iprot)=ftbis(2)
228 escalbis(8,ib,ibatch,iprot)=ftbis(2)
229 escalbis(9,ib,ibatch,iprot)=ftbis(3)
230 escalbis(10,ib,ibatch,iprot)=ftbis(5)
231 escalbis(13,ib,ibatch,iprot)=ftbis(1)
232 escalbis(14,ib,ibatch,iprot)=ftbis(2)
233 escalbis(19,ib,ibatch,iprot)=ftbis(1)
236 betmin(ibatch,iprot)=betaT(1,ibatch,iprot)
237 betmax(ibatch,iprot)=betaT(1,ibatch,iprot)
238 do ib=2,nbeta(ibatch,iprot)
239 if (betaT(ib,ibatch,iprot).lt.betmin(ibatch,iprot))
240 & betmin(ibatch,iprot)=betaT(ib,ibatch,iprot)
241 if (betaT(ib,ibatch,iprot).gt.betmax(ibatch,iprot))
242 & betmax(ibatch,iprot)=betaT(ib,ibatch,iprot)
244 write (iout,'(2(a,i5,2x),2(a,f7.2,2x))')
245 & "Protein",iprot," batch",ibatch,
246 & " betmin",betmin(ibatch,iprot)," betmax",betmax(ibatch,iprot)
252 c Compute energy and entropy histograms to filter out conformations
253 c with potntially too low or too high weights.
255 call determine_histograms
257 c Make the working list of conformations
259 call make_list(.true.,.true.,nvarr,x)
260 c 3/8/2013 Calculate RMSDs/Qs and the Ptab array
263 call restore_molinfo(iprot)
265 open (icbase,file=bprotfiles(iprot),status="old",
266 & form="unformatted",access="direct",recl=lenrec(iprot))
267 entfacmax=entfac(1,iprot)
268 sument=entfac(1,iprot)
269 sumentsq=entfac(1,iprot)**2
270 do i=2,ntot_work(iprot)
271 if (entfac(i,iprot).gt.entfacmax) entfacmax=entfac(i,iprot)
272 sument=sument+entfac(i,iprot)
273 sumentsq=sumentsq+entfac(i,iprot)**2
275 sument=sument/ntot_work(iprot)
276 sumentsq=dsqrt(sumentsq/ntot_work(iprot)-sument**2)
277 write (2,*) "entfacmax",entfacmax," entave",sument,
281 nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
285 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
286 ipass_conf=ipass_conf+1
287 write (iout,*) "Pass",ipass_conf
289 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
291 nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
294 do i=1,ntot(iprot),maxstr_proc
295 ipass_conf=ipass_conf+1
296 write (iout,*) "Pass",ipass_conf
298 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
301 c Read the chunk of conformations off a DA scratchfile.
303 call daread_ccoords(iprot,istart_conf,iend_conf)
305 IF (GEOM_AND_WHAM_WEIGHTS) THEN
308 do ib=1,nbeta(1,iprot)
309 write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
311 do iii=istart_conf,iend_conf
317 etoti=etoti+ww(k)*escal(k,ib,1,iprot)
320 e_total_T(iii,ib)=entfac(kk,iprot)
321 & +betaT(ib,1,iprot)*(etoti-elowest(ib,1,iprot))
325 write (iout,*) "Energies before Gather"
326 do iii=1,ntot_work(iprot)
327 write (iout,'(i5,100E12.5)') iii,
328 & (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
332 do ib=1,nbeta(1,iprot)
333 do k=1,scount(me1,iprot)
334 e_total_T_(k)=e_total_T(indstart(me1,iprot)+k-1,ib)
336 c call MPI_AllGatherv(e_total_T(indstart(me1,iprot),ib),
337 c & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
338 c & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
340 call MPI_AllGatherv(e_total_T_(1),
341 & scount(me1,iprot),MPI_DOUBLE_PRECISION,e_total_T(1,ib),
342 & scount(0,iprot),idispl(0,iprot),MPI_DOUBLE_PRECISION,
346 write (iout,*) "Energies after Gather"
347 do iii=1,ntot_work(iprot)
348 write (iout,'(i5,100E12.5)') iii,
349 & (e_total_T(iii,ib),ib=1,nbeta(1,iprot))
354 ENDIF ! GEOM_AND_WHAM_WEIGHTS
357 c write (iout,*) "Start calculating weirms",istart_conf,iend_conf
358 do iii=istart_conf,iend_conf
361 c write (iout,*) "Calculating weirms iii",iii," ii",ii," k",k
363 call restore_ccoords(iprot,ii)
364 c Calculate the rmsds of the current conformations and convert them into weights
365 c to approximate equal weighting of all points of the conformational space
371 c 7/6/16 AL Add angle comparison
372 if (anglecomp(iprot) .or. sidecomp(iprot)) then
373 call int_from_cart(.true.,.false.)
374 call sc_loc_geom(.false.)
377 thetareff(k)=theta(k)
383 do ib=1,nbeta(1,iprot)
386 do iref=1,ntot_work(iprot)
387 if (iref.ne.iii) then
388 read(icbase,rec=iref) ((csingle(l,k),l=1,3),k=1,nres),
389 & ((csingle(l,k),l=1,3),k=nnt+nres,nct+nres),
390 & nss,(ihpb(k),jhpb(k),k=1,nss),eini(iref,iprot),
391 & entfac(iref,iprot),
392 & ((nu(k,l,ii,iprot),k=1,maxdimnat),l=1,natlike(iprot))
400 write (iout,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
405 if (rmscomp(iprot)) then
406 call fitsq(rms,c(1,nstart_sup(iprot)),
407 & creff(1,nstart_sup(iprot)),
408 & nsup(iprot),przes,obrot,non_conv)
410 print *,'Error: FITSQ non-convergent, iref',iref
412 else if (rms.lt.-1.0d-6) then
413 print *,'Error: rms^2 = ',rms,iref
415 else if (rms.lt.0) then
421 write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
422 & e_total_T(iref,ib),
423 & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
424 & " rms",rms," fac",dexp(-0.5d0*(rms/sigma2(iprot))**2),
425 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
426 & dexp(-0.5d0*(rms/sigma2(iprot))**2)*
427 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
429 aux=dexp(-0.5d0*(rms/sigma2(iprot))**2)
431 rms=qwolyness(creff,iprot)
433 write (2,*) "iii",iii," iref=",iref," ib",ib," etotal",
434 & e_total_T(iref,ib),
435 & e_total_T(iii,ib),-e_total_T(iref,ib)+e_total_T(iii,ib),
436 & " rms",-dlog(rms)," fac",rms,
437 & dexp(-e_total_T(iref,ib)+e_total_T(iii,ib)),
438 & rms*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
442 c 7/6/16 AL Add angle comparison
443 if (anglecomp(iprot).or.sidecomp(iprot))
444 & call int_from_cart(.true.,.false.)
445 if (anglecomp(iprot)) then
448 dtheta = theta(k)-thetareff(k)
449 rmsthet = rmsthet+dtheta*dtheta
451 rmsthet=rmsthet/(nct-nnt-2)
454 dphi = pinorm(phi(k)-phireff(k))
455 rmsphi = rmsphi + dphi*dphi
457 rmsphi=rmsphi/(nct-nnt-3)
458 aux=aux*dexp(-0.5d0*(rmsthet+rmsphi)/sigmaang2(iprot)**2)
460 if (sidecomp(iprot)) then
461 call sc_loc_geom(.false.)
464 dxref = xxtab(k)-xxreff(k)
465 dyref = yytab(k)-yyreff(k)
466 dzref = zztab(k)-zzreff(k)
467 rmsside = rmsside + dx*dx+dy*dy+dz*dz
469 rmsside=rmsside/(nct-nnt+1)
470 aux=aux*dexp(-0.5d0*rmsside/sigmaside2(iprot)**2)
475 do ib=1,nbeta(1,iprot)
476 if (GEOM_AND_WHAM_WEIGHTS) then
477 weirms(iii,ib) = weirms(iii,ib)
478 & +aux*dexp(-e_total_T(iref,ib)+e_total_T(iii,ib))
480 weirms(iii,ib) = weirms(iii,ib)+aux
486 write (iout,*) "weirms"
487 do iii=istart_conf,iend_conf
488 write (iout,"(i5,20e12.5)") iii,(weirms(iii,ib),ib=1,
496 nchunk_conf(iprot) = iroof(scount(me1,iprot),maxstr_proc)
497 DO IB=1,NBETA(1,IPROT)
499 write(csa_bank,'(a,f4.0,4hPtab,i4.4)')
500 & protname(iprot)(:ilen(protname(iprot))),
501 & 1.0d0/(0.001987*betaT(ib,1,iprot)),me
502 open (icsa_bank,file=csa_bank,status="unknown")
507 do i=indstart(me1,iprot),indend(me1,iprot),maxstr_proc
508 ipass_conf=ipass_conf+1
509 write (iout,*) "Pass",ipass_conf
511 iend_conf=min0(i+maxstr_proc-1,indend(me1,iprot))
513 nchunk_conf(iprot) = iroof(ntot(iprot),maxstr_proc)
516 do i=1,ntot(iprot),maxstr_proc
517 ipass_conf=ipass_conf+1
518 write (iout,*) "Pass",ipass_conf
520 iend_conf=min0(i+maxstr_proc-1,ntot(iprot))
524 c Read the chunk of conformations off a DA scratchfile.
526 call daread_ccoords(iprot,istart_conf,iend_conf)
527 write (iout,*) "IB",ib," IPROT",iprot," NREF",nref(ib,iprot)
528 do iii=istart_conf,iend_conf
536 call restore_ccoords(iprot,ii)
537 c write (iout,*) "zeroing Ptab",iii,ii,jj
538 Ptab(jj,ib,iprot)=0.0d0
539 if (rmscomp(iprot)) then
540 do iref=1,nref(ib,iprot)
544 call pdbout(0.0d0,"conf")
545 c write (2,*) "nstart_sup",nstart_sup(iprot),
546 c & " nend_sup",nend_sup(iprot)
547 c do k=nstart_sup(iprot),nend_sup(iprot)
548 c write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k),l=1,3),
549 c & (cref(l,k,iref,ib,iprot),l=1,3)
551 c do k=nstart_sup(iprot),nend_sup(iprot)
552 c write(2,'(i5,3f10.5,5x,3f10.5)') k,(c(l,k+nres),l=1,3),
553 c & (cref(l,k+nres,iref,ib,iprot),l=1,3)
558 c(l,k)=cref(l,k,iref,ib,iprot)
561 call pdbout(0.0d0,"ref")
569 rms=rmsnat(ib,jj,iref,iprot)
570 c 7/6/16 AL Add angle comparison
571 if (anglecomp(iprot).or.sidecomp(iprot))
572 & call int_from_cart(.true.,.false.)
573 if (anglecomp(iprot)) then
574 c write (iout,*) "jj",jj," iref",iref
577 dtheta = theta(k)-theta_ref(k,iref,ib,iprot)
578 c write (iout,*) k,theta(k),theta_ref(k,iref,ib,iprot),
580 rmsthet = rmsthet+dtheta*dtheta
582 rmsthet=rmsthet/(nct-nnt-2)
585 dphi = pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
586 c write (iout,*) k,phi(k),phi_ref(k,iref,ib,iprot),
587 c & pinorm(phi(k)-phi_ref(k,iref,ib,iprot))
588 rmsphi = rmsphi + dphi*dphi
590 rmsphi=rmsphi/(nct-nnt-3)
592 if (sidecomp(iprot)) then
593 call sc_loc_geom(.false.)
596 dxref = xxtab(k)-xxref(k,iref,ib,iprot)
597 dyref = yytab(k)-yyref(k,iref,ib,iprot)
598 dzref = zztab(k)-zzref(k,iref,ib,iprot)
599 rmsside = rmsside + dx*dx+dy*dy+dz*dz
601 rmsside=rmsside/(nct-nnt+1)
604 write (iout,*) "ii",ii," jj",jj," iref",iref,
606 & "efree",efreeref(iref,ib,iprot)
610 if (rms.lt.rmsmin) then
615 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)
616 & +dexp(-efreeref(iref,ib,iprot)
617 & -0.5d0*(rms/sigma2(iprot))**2)
618 if (anglecomp(iprot))
619 & Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
620 & dexp(-0.5d0*(rmsthet+rmsphi)/sigmaang2(iprot)**2)
622 & Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*
623 & dexp(-0.5d0*rmsside/sigmaside2(iprot)**2)
626 rmsave=rmsave/nref(ib,iprot)
629 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
630 c search of conformational space.
631 aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
632 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
634 write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
635 c & " entfac",entfac(k,iprot),"aux",aux," Ptab",Ptab(jj,ib,iprot)
636 & "entfac",entfac(k,iprot)," weirms",weirms(iii,ib),
637 & " Ptab",Ptab(jj,ib,iprot)
639 #elif defined(WEIDIST)
640 c write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
641 c & " weirms",weirms(iii,ib)," Ptab/weirms",
642 c & Ptab(jj,ib,iprot)/weirms(iii,ib)
644 write (icsa_bank,'(2i6,3f10.5,2f8.4)')
645 & iii,imin,-dlog(Ptab(jj,ib,iprot)),
646 & -dlog(weirms(iii,ib)),
647 & -dlog(Ptab(jj,ib,iprot)/weirms(iii,ib)),
650 if (ib.eq.1) rmstb(iii,iprot)=rmsmin
651 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
652 #elif defined(WEICLUST)
653 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,iprot)
655 sumP=sumP+Ptab(jj,ib,iprot)
657 do iref=1,nref(ib,iprot)
658 rms = qwolynes(ib,iref,iprot)
660 write (iout,*) "iii",iii," ii",ii," jj",jj,
661 & "iref",iref," rms",rms,-dlog(rms),
662 & "efree",efreeref(iref,ib,iprot)
664 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)+rms*
665 & *dexp(-efreeref(iref,ib,iprot))
668 c AL 1/24/2014 Correct the weight of a conformation to represent energy-free
669 c search of conformational space.
670 aux=dmin1(entfacmax-entfac(k,iprot),60.0d0)
671 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)*dexp(aux)
673 write (iout,*) "iprot",iprot," ib",ib," iii",iii," k",k,
674 & " entfac",entfac(k,iprot)," aux",aux," Ptab",Ptab(jj,ib,iprot)
676 #elif defined(WEICLUST)
677 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(k,iprot)
678 #elif defined(WEIDIST)
679 write (iout,*) "iii",iii," Ptab",Ptab(jj,ib,iprot),
680 & " weirms",weirms(iii,ib)," Ptab/weirms",
681 & Ptab(jj,ib,iprot)/weirms(iii,ib)
682 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/weirms(iii,ib)
684 sumP=sumP+Ptab(jj,ib,iprot)
687 write (iout,*) "ii",ii," jj",jj," ib",ib," iprot",iprot,
688 & " Ptab",Ptab(jj,ib,iprot)," sumP",sumP
693 write (iout,*) "iprot",iprot," ib",ib," sumP before REDUCE",sumP
696 call MPI_Allreduce(sumP,sumP_t,1,MPI_DOUBLE_PRECISION,
697 & MPI_SUM,Comm1,IERROR)
701 write (iout,*) "iprot",iprot," ib",ib," sumP after REDUCE",sumP
704 do iii=indstart(me1,iprot),indend(me1,iprot)
707 write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
709 Ptab(jj,ib,iprot)=Ptab(jj,ib,iprot)/sumP
711 write (iout,*) "iii",iii," jj",jj," Ptab",Ptab(jj,ib,iprot)
721 write (iout,*) "ELOWEST at the end of PROC_DATA"
723 do ibatch=1,natlike(iprot)+2
724 do ib=1,nbeta(ibatch,iprot)
725 write (iout,*) "iprot",iprot," ibatch",ibatch," ib",ib,
726 & " elowest",elowest(ib,ibatch,iprot)
730 write (iout,*) "Ptab at the end of PROC"
732 do ib=1,nbeta(1,iprot)
733 write (iout,*) "Protein",iprot," beta",ib
735 do i=indstart(me1,iprot),indend(me1,iprot)
737 write (iout,*) iprot,ib,ii,jj,Ptab(jj,ib,iprot)
746 c-------------------------------------------------------------------------------
747 subroutine determine_histograms
750 include "DIMENSIONS.ZSCOPT"
752 parameter (maxind=1000)
755 integer Ierror,Errcode,status(MPI_STATUS_SIZE),tag
758 include "COMMON.COMPAR"
759 include "COMMON.ENERGIES"
760 include "COMMON.VMCPAR"
761 include "COMMON.WEIGHTS"
762 include "COMMON.PROTNAME"
763 include "COMMON.IOUNITS"
764 include "COMMON.FFIELD"
765 include "COMMON.PROTFILES"
766 include "COMMON.CHAIN"
767 include "COMMON.CONTROL"
768 include "COMMON.CLASSES"
769 integer iprot,i,ii,iii,j,k,l,ibatch,inde,indt,idelte,ideltt
770 double precision eminh,emaxh,tminh,tmaxh
771 integer he(0:maxind),
772 & ht(0:maxind),hemax,htmax
774 c Construct energy and entropy histograms
776 write (iout,*) "Calling DETERMINE HISTOGRAMS"
789 if (eini(iii,iprot).lt.eminh)
790 & eminh=eini(iii,iprot)
791 if (eini(iii,iprot).gt.emaxh)
792 & emaxh=eini(iii,iprot)
793 if (entfac(iii,iprot).lt.tminh)
794 & tminh=entfac(iii,iprot)
795 if (entfac(iii,iprot).gt.tmaxh)
796 & tmaxh=entfac(iii,iprot)
798 c write (2,*) "DETERMINE HISTOGRAMS 1"
802 inde=eini(iii,iprot)-eminh
803 indt=entfac(iii,iprot)-tminh
804 if (inde.le.maxind) he(inde)=he(inde)+1
805 if (indt.le.maxind) ht(indt)=ht(indt)+1
807 c write (2,*) "DETERMINE HISTOGRAMS 2"
811 c write (iout,*) "idelte",idelte," ideltt",ideltt
814 c write (iout,*) "i",i," he",he(i)," hemax",hemax
816 if (he(i).gt.hemax) hemax=he(i)
820 c write (iout,*) "i",i," ht",ht(i)," htmax",htmax
822 if (ht(i).gt.htmax) htmax=ht(i)
824 c write (2,*) "DETERMINE HISTOGRAMS 3"
826 hemax=hemax/hefac(iprot)
827 if (hemax.lt.hemax_low(iprot))
828 & hemax=hemax_low(iprot)
829 htmax=htmax/htfac(iprot)
830 if (htmax.lt.htmax_low(iprot))
831 & htmax=htmax_low(iprot)
833 do while (i.lt.idelte .and. he(i).lt.hemax)
836 c write (2,*) "DETERMINE HISTOGRAMS 4"
838 e_lowb(iprot)=eminh+i
840 do while (i.lt.ideltt .and. ht(i).lt.htmax)
843 t_lowb(iprot)=tminh+i
845 write (iout,*) "protein",iprot
846 write (iout,*) "energy histogram"
848 write (iout,'(f15.5,i10)') eminh+i,he(i)
850 write (iout,*) "entropy histogram"
852 write (iout,'(f15.5,i10)') tminh+i,ht(i)
854 write (iout,*) "emin",eminh," tmin",tminh,
855 & " hemax",hemax," htmax",htmax,
856 & " e_lowb",e_lowb(iprot),
857 & " t_lowb",t_lowb(iprot)
863 c------------------------------------------------------------------------------
864 subroutine imysort(n, x, ipermut)
874 if (x(j).lt.xtemp) then
882 ipermut(imax)=ipermut(i)
887 c--------------------------------------------------------------------------
888 subroutine write_conf_count
891 include "DIMENSIONS.ZSCOPT"
892 include "COMMON.CLASSES"
893 include "COMMON.COMPAR"
894 include "COMMON.VMCPAR"
895 include "COMMON.PROTNAME"
896 include "COMMON.IOUNITS"
897 include "COMMON.PROTFILES"
898 integer iprot,ibatch,i,ii,j,kk
903 write (iout,*)"Numbers of conformations after applying the cutoff"
906 write (iout,'(a,2x,a,i10)') "Protein",
907 & protname(iprot)(:ilen(protname(iprot))),
912 c-------------------------------------------------------------------------
913 double precision function fdum()