1 SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
2 implicit real*8 (a-h,o-z)
5 parameter (num_in_line=5)
6 LOGICAL PRINTANG(max_cut)
7 integer PRINTPDB(max_cut),printmol2(max_cut)
8 include 'COMMON.CONTROL'
9 include 'COMMON.HEADER'
10 include 'COMMON.CHAIN'
12 include 'COMMON.CLUSTER'
13 include 'COMMON.IOUNITS'
16 include 'COMMON.TEMPFAC'
17 include 'COMMON.FFIELD'
18 include 'COMMON.SBRIDGE'
19 include 'COMMON.TORCNSTR'
21 CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
24 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
27 integer ib,list_peak_viol(maxdim_cont)
28 double precision Esaxs_all(maxgr),Pcalc_all(maxsaxs,maxgr)
33 c print *,"calling WRTCLUST",ncon
34 c write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
37 temper=1.0d0/(beta_h(ib)*1.987d-3)
38 if (temper.lt.100.0d0) then
39 write(ctemper,'(f3.0)') temper
41 else if (temper.lt.1000.0) then
42 write (ctemper,'(f4.0)') temper
45 write (ctemper,'(f5.0)') temper
49 do i=1,ncon*(ncon-1)/2
52 close(80,status='delete')
54 C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
56 ii1= index(intinname,'/')
61 ii2=index(intinname(ii1:),'/')
63 ii = ii1+index(intinname(ii1:),'.')-1
69 prefixp=intinname(ii1:ii)
70 cd print *,icut,printang(icut),printpdb(icut),printmol2(icut)
71 cd print *,'ecut=',ecut
74 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
75 NRECORD=LICZ(IGR)/num_in_line
77 DO 63 IRECORD=1,NRECORD
78 IND2=IND1+num_in_line-1
79 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
80 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
83 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
84 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
86 ICON=list_conf(NCONF(IGR,1))
87 c WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
88 C 12/8/93 Estimation of "diameters" of the subsequent families.
91 c write (iout,*) "ecut",ecut
92 emin=totfree(nconf(igr,1))
95 if (totfree(ii)-emin .gt. ecut) goto 10
100 ind=ioffset(ncon,ii,jj)
102 ind=ioffset(ncon,jj,ii)
104 c write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
105 c & " ind",ind," diss",diss(ind)
107 curr_dist=dabs(diss(ind)+0.0d0)
108 c write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
109 c & list_conf(jj),curr_dist
110 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
111 ave_dim=ave_dim+curr_dist**2
114 10 if (licz(igr) .gt. 1)
115 & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
116 write (iout,'(/A,F8.1,A,F8.1)')
117 & 'Max. distance in the family:',amax_dim,
118 & '; average distance in the family:',ave_dim
120 gdt_ts_ave(igr)=0.0d0
121 gdt_ha_ave(igr)=0.0d0
122 tmscore_ave(igr)=0.0d0
124 e1=totfree(nconf(igr,1))
127 boltz=dexp(-(totfree(icon)-e1))
128 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
129 gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
130 gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
131 tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
133 c write (iout,'(2i5,10f10.5)') i,icon,boltz,rmstb(icon),
134 c & gdt_ts_tb(icon),gdt_ha_tb(icon),tmscore_tb(icon)
136 c write (iout,*) "qpart",qpart
137 rmsave(igr)=rmsave(igr)/qpart
138 gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
139 gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
140 tmscore_ave(igr)=tmscore_ave(igr)/qpart
141 write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
142 & rmsave(igr)," A, ",
143 & "TMscore",tmscore_ave(igr),
144 & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
148 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
149 c print *,icut,printang(icut)
150 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
152 c print *,'emin',emin,' ngr',ngr
153 if (lprint_cart) then
154 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
157 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
162 if (totfree_gr(igr)-emin.le.ecut) then
163 if (lprint_cart) then
164 call cartout(igr,icon,totfree(icon)/beta_h(ib),
165 & totfree_gr(igr)/beta_h(ib),
166 & rmstb(icon),cfname)
168 c print '(a)','calling briefout'
171 c(j,i)=allcart(j,i,icon)
174 call int_from_cart1(.false.)
175 call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
176 & totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
177 & jhpb_all(1,icon),cfname)
178 c print '(a)','exit briefout'
184 IF (PRINTPDB(ICUT).gt.0) THEN
185 c Write out a number of conformations from each family in PDB format and
186 c create InsightII command file for their displaying in different colors
187 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
188 & //"K_"//'ave'//exten
189 write (iout,*) "cfname",cfname
190 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
191 write (ipdb,'(a,f8.2)')
192 & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
198 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
199 c write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
200 write (NUMM,'(bz,i4.4)') i
201 ncon_lim=min0(licz(i),printpdb(icut))
202 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
203 & //"K_"//numm(:ilen(numm))//exten
204 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
205 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
206 & " AVE RMSD",0pf5.2)')
207 & i,totfree_gr(i)/beta_h(ib),rmsave(i)
208 c Write conformations of the family i to PDB files
210 do while (ncon_out.lt.printpdb(icut) .and.
211 & ncon_out.lt.licz(i).and.
212 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
214 c write (iout,*) i,ncon_out,nconf(i,ncon_out),
215 c & totfree(nconf(i,ncon_out)),emin1,ecut
217 c write (iout,*) "ncon_out",ncon_out
227 c(k,ii)=allcart(k,ii,icon)
231 call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
232 write (ipdb,'("TER")')
235 c Average structures and structures closest to average
236 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
237 & //"K_"//'ave'//exten
238 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
241 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
243 call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
244 write (ipdb,'("TER")')
245 if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0
246 & .or.npeak.gt.0)) then
247 call fittest_coord(i)
249 call closest_coord(i)
251 c write (iout,*) "Calling rmsnat"
252 rms_closest(i) = rmsnat(i)
254 write (iout,*) "Cluster",i
255 call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
256 & tmscore_closest(i),cfname,.true.)
257 c write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i
258 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
259 call e_saxs(Esaxs_constr)
262 Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
263 & (Pcalc(j+1)+Pcalc(j))/2
266 Pcalc_all(j,i)=Pcalc(j)/Cnorm
268 c write (iout,*) "Pcalc"
269 c write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
270 Esaxs_all(i)=Esaxs_constr
271 write (iout,*) "Esaxs",Esaxs_constr
274 if (link_start.gt.0) then
275 do j=link_start,link_end
276 if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
277 dxlink=dist(ihpb(j),jhpb(j))
278 if (dxlink.le.25.0d0) then
279 write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
280 & irestr_type(j),ihpb(j),jhpb(j),
283 nviolxlink=nviolxlink+1
284 write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
285 & irestr_type(j),ihpb(j),jhpb(j),
291 & write (iout,*) nviolxlink," crosslink violations."
292 c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
293 c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
294 c & "tmscore",tmscore_closest(i)
296 c Determine # violated NMR restraints
297 if (link_end_peak.gt.0) then
299 write (NUMM,'(bz,i4.4)') i
300 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
301 & //"K_"//NUMM(:ilen(NUMM))//'.nmr'
302 open(jrms,file=cfname)
303 do j=link_start_peak,link_end_peak
305 do ip=ipeak(1,j),ipeak(2,j)
309 c iip=ip-ipeak(1,j)+1
310 C iii and jjj point to the residues for which the distance is assigned.
320 if (dd.lt.dhpb1_peak(ip)) then
322 c write (iout,*) j,iii,jjj,iiib
323 write (jrms,'(4i6)') j,iii,jjj,iiib
327 nviolpeak=nviolpeak+1
328 list_peak_viol(nviolpeak)=j
331 if (nviolpeak.gt.0) then
332 write (iout,'(a,i5,2h (f8.4,2h%))')
333 & "Number of violated NMR restraints:",
334 & nviolpeak,100*(nviolpeak+0.)/npeak
335 write (iout,'(a)')"List of violated restraints:"
336 write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak)
340 if (.not.raw_psipred .and. idihconstr_end.gt.0) then
341 cfname=prefixp(:ilen(prefixp))//"_T"
342 & //ctemper(:ilen(ctemper))
343 & //"K_"//NUMM(:ilen(NUMM))//'.angle'
344 open(jrms,file=cfname)
345 call int_from_cart1(.false.)
347 do j=idihconstr_start,idihconstr_end
350 difi=pinorm(phii-phi0(j))
351 if (difi.gt.drange(j) .or. difi.lt.-drange(j))
352 & nangviol=nangviol+1
353 write (jrms,'(i5,3f10.3)') itori,phii*rad2deg,
354 & phi0(j)*rad2deg,rad2deg*drange(j)
356 write (iout,'(a,i5)')"Number of angle-restraint violations:"
361 call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
362 write (ipdb,'("TER")')
369 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
370 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
371 & //"K_"//'ave'//'.dist'
372 OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
373 write (99,'(5h# ,10f10.5)')
374 & (Esaxs_all(i)*wsaxs,i=1,ngr_print)
376 write (99,'(f6.2,10f10.5)') distsaxs(j),
377 & (Pcalc_all(j,i),i=1,ngr_print)
382 IF (printmol2(icut).gt.0) THEN
383 c Write out a number of conformations from each family in PDB format and
384 c create InsightII command file for their displaying in different colors
389 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
390 write (NUMM,'(bz,i4.4)') i
391 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
392 & //"K_"//numm(:ilen(numm))//extmol
393 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
395 do while (ncon_out.lt.printmol2(icut) .and.
396 & ncon_out.lt.licz(i).and.
397 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
404 c(k,ii)=allcart(k,ii,icon)
407 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
415 call WRITE_STATS(ICUT,NCON,IB)
416 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
417 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
418 & ' CONTAINS ',I4,' CONFORMATION(S): ')
419 c 300 FORMAT ( 8(I4,F6.1))
420 300 FORMAT (5(I4,1pe12.3))
421 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
422 500 FORMAT (8(2I4,2X))
423 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
426 c------------------------------------------------------------------------------
427 subroutine ave_coord(igr)
430 include 'sizesclu.dat'
431 include 'COMMON.CONTROL'
432 include 'COMMON.CLUSTER'
433 include 'COMMON.CHAIN'
434 include 'COMMON.INTERACT'
436 include 'COMMON.TEMPFAC'
437 include 'COMMON.IOUNITS'
439 double precision przes(3),obrot(3,3)
440 double precision xx(3,maxres2),csq(3,maxres2)
441 double precision eref
442 double precision rmscalc
443 c double precision rmscheck
444 integer i,ii,j,k,icon,jcon,igr,ipermmin
445 double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
446 c write (iout,*) "AVE_COORD: igr",igr
449 boltz = dexp(-totfree(jcon)+eref)
453 c(j,i)=allcart(j,i,jcon)*boltz
454 cref1(j,i)=allcart(j,i,jcon)
455 csq(j,i)=allcart(j,i,jcon)**2*boltz
460 c write (iout,*) "k",k," jcon",jcon
463 cwork(j,i)=allcart(j,i,jcon)
466 rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin)
467 c write (iout,*) "rms",rms," ipermmin",ipermmin
469 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),
470 c & (obrot(i,j),j=1,3)
472 c if (rms.lt.0.0) then
473 c print *,'error, rms^2 = ',rms,icon,jcon
476 c if (non_conv) print *,non_conv,icon,jcon
477 boltz=dexp(-totfree(jcon)+eref)
478 qpart = qpart + boltz
481 xx(j,i)=allcart(j,i,jcon)
484 call matvec(cwork,obrot,xx,2*nres)
486 c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
487 c & (allcart(j,i,jcon),j=1,3)
489 cwork(j,i)=cwork(j,i)+przes(j)
490 c(j,i)=c(j,i)+cwork(j,i)*boltz
491 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
498 c rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2
501 c write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms
506 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
508 c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
514 tempfac(1,i)=tempfac(1,i)+csq(j,i)
515 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
517 tempfac(1,i)=dsqrt(tempfac(1,i))
518 tempfac(2,i)=dsqrt(tempfac(2,i))
522 c------------------------------------------------------------------------------
523 subroutine fittest_coord(igr)
526 include 'sizesclu.dat'
527 include 'COMMON.IOUNITS'
528 include 'COMMON.CONTROL'
529 include 'COMMON.CLUSTER'
530 include 'COMMON.CHAIN'
531 include 'COMMON.INTERACT'
533 include 'COMMON.FFIELD'
534 include 'COMMON.TORCNSTR'
535 include 'COMMON.SAXS'
537 double precision przes(3),obrot(3,3)
538 double precision xx(3,maxres2),yy(3,maxres2)
539 integer i,ii,j,k,icon,jcon,jconmin,igr
540 double precision rms,rmsmin,cwork(3,maxres2)
541 double precision ehpb,Esaxs_constr,edihcnstr
548 c(j,i)=allcart(j,i,jcon)
551 call int_from_cart1(.false.)
555 if (nsaxs.gt.0) call e_saxs(Esaxs_constr)
557 if (ndih_constr.gt.0) call etor_constr(edihcnstr)
558 rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr
559 c write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb,
560 c & " Edihcnstr",edihcnstr
561 if (rms.lt.rmsmin) then
566 write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
569 c(j,i)=allcart(j,i,jconmin)
574 c------------------------------------------------------------------------------
575 subroutine closest_coord(igr)
578 include 'sizesclu.dat'
579 include 'COMMON.IOUNITS'
580 include 'COMMON.CONTROL'
581 include 'COMMON.CLUSTER'
582 include 'COMMON.CHAIN'
583 include 'COMMON.INTERACT'
586 double precision przes(3),obrot(3,3)
587 integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin
588 double precision rms,rmsmin,cwork(3,maxres2)
589 double precision xx(3,maxres2),yy(3,maxres2)
590 double precision rmscalc
598 yy(j,i)=allcart(j,i,jcon)
601 rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin)
602 c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
603 if (non_conv) print *,non_conv,icon,jcon
604 if (rms.lt.rmsmin) then
609 c write (iout,*) "rmsmin",rmsmin," rms",rms
613 c(j,i)=allcart(j,i,jconmin)
618 c------------------------------------------------------------------------------
622 include 'sizesclu.dat'
623 include 'COMMON.IOUNITS'
624 include 'COMMON.CONTROL'
625 include 'COMMON.CLUSTER'
626 include 'COMMON.CHAIN'
627 include 'COMMON.INTERACT'
629 double precision przes(3)
630 integer i,ii,j,k,icon,jcon,jconmin,igr
634 przes(j)=przes(j)+c(j,i)
638 przes(j)=przes(j)/nres
642 c(j,i)=c(j,i)-przes(j)