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'
20 CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
23 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
26 integer ib,list_peak_viol(maxdim)
27 double precision Esaxs_all(maxgr),Pcalc_all(maxsaxs,maxgr)
32 c print *,"calling WRTCLUST",ncon
33 c write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
36 temper=1.0d0/(beta_h(ib)*1.987d-3)
37 if (temper.lt.100.0d0) then
38 write(ctemper,'(f3.0)') temper
40 else if (temper.lt.1000.0) then
41 write (ctemper,'(f4.0)') temper
44 write (ctemper,'(f5.0)') temper
48 do i=1,ncon*(ncon-1)/2
51 close(80,status='delete')
53 C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
55 ii1= index(intinname,'/')
60 ii2=index(intinname(ii1:),'/')
62 ii = ii1+index(intinname(ii1:),'.')-1
68 prefixp=intinname(ii1:ii)
69 cd print *,icut,printang(icut),printpdb(icut),printmol2(icut)
70 cd print *,'ecut=',ecut
73 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
74 NRECORD=LICZ(IGR)/num_in_line
76 DO 63 IRECORD=1,NRECORD
77 IND2=IND1+num_in_line-1
78 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
79 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
82 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
83 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
85 ICON=list_conf(NCONF(IGR,1))
86 c WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
87 C 12/8/93 Estimation of "diameters" of the subsequent families.
90 c write (iout,*) "ecut",ecut
91 emin=totfree(nconf(igr,1))
94 if (totfree(ii)-emin .gt. ecut) goto 10
99 ind=ioffset(ncon,ii,jj)
101 ind=ioffset(ncon,jj,ii)
103 c write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
104 c & " ind",ind," diss",diss(ind)
106 curr_dist=dabs(diss(ind)+0.0d0)
107 c write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
108 c & list_conf(jj),curr_dist
109 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
110 ave_dim=ave_dim+curr_dist**2
113 10 if (licz(igr) .gt. 1)
114 & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
115 write (iout,'(/A,F8.1,A,F8.1)')
116 & 'Max. distance in the family:',amax_dim,
117 & '; average distance in the family:',ave_dim
119 gdt_ts_ave(igr)=0.0d0
120 gdt_ha_ave(igr)=0.0d0
121 tmscore_ave(igr)=0.0d0
123 e1=totfree(nconf(igr,1))
126 boltz=dexp(-(totfree(icon)-e1))
127 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
128 gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
129 gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
130 tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
132 c write (iout,'(2i5,10f10.5)') i,icon,boltz,rmstb(icon),
133 c & gdt_ts_tb(icon),gdt_ha_tb(icon),tmscore_tb(icon)
135 c write (iout,*) "qpart",qpart
136 rmsave(igr)=rmsave(igr)/qpart
137 gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
138 gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
139 tmscore_ave(igr)=tmscore_ave(igr)/qpart
140 write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
141 & rmsave(igr)," A, ",
142 & "TMscore",tmscore_ave(igr),
143 & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
147 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
148 c print *,icut,printang(icut)
149 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
151 c print *,'emin',emin,' ngr',ngr
152 if (lprint_cart) then
153 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
156 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
161 if (totfree_gr(igr)-emin.le.ecut) then
162 if (lprint_cart) then
163 call cartout(igr,icon,totfree(icon)/beta_h(ib),
164 & totfree_gr(igr)/beta_h(ib),
165 & rmstb(icon),cfname)
167 c print '(a)','calling briefout'
170 c(j,i)=allcart(j,i,icon)
173 call int_from_cart1(.false.)
174 call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
175 & totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
176 & jhpb_all(1,icon),cfname)
177 c print '(a)','exit briefout'
183 IF (PRINTPDB(ICUT).gt.0) THEN
184 c Write out a number of conformations from each family in PDB format and
185 c create InsightII command file for their displaying in different colors
186 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
187 & //"K_"//'ave'//exten
188 write (iout,*) "cfname",cfname
189 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
190 write (ipdb,'(a,f8.2)')
191 & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
197 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
198 c write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
199 write (NUMM,'(bz,i4.4)') i
200 ncon_lim=min0(licz(i),printpdb(icut))
201 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
202 & //"K_"//numm(:ilen(numm))//exten
203 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
204 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
205 & " AVE RMSD",0pf5.2)')
206 & i,totfree_gr(i)/beta_h(ib),rmsave(i)
207 c Write conformations of the family i to PDB files
209 do while (ncon_out.lt.printpdb(icut) .and.
210 & ncon_out.lt.licz(i).and.
211 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
213 c write (iout,*) i,ncon_out,nconf(i,ncon_out),
214 c & totfree(nconf(i,ncon_out)),emin1,ecut
216 c write (iout,*) "ncon_out",ncon_out
226 c(k,ii)=allcart(k,ii,icon)
230 call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
231 write (ipdb,'("TER")')
234 c Average structures and structures closest to average
235 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
236 & //"K_"//'ave'//exten
237 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
240 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
242 call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
243 write (ipdb,'("TER")')
244 if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0
245 & .or.npeak.gt.0)) then
246 call fittest_coord(i)
248 call closest_coord(i)
250 c write (iout,*) "Calling rmsnat"
251 rms_closest(i) = rmsnat(i)
253 write (iout,*) "Cluster",i
254 call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
255 & tmscore_closest(i),cfname,.true.)
256 c write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i
257 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
258 call e_saxs(Esaxs_constr)
261 Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
262 & (Pcalc(j+1)+Pcalc(j))/2
265 Pcalc_all(j,i)=Pcalc(j)/Cnorm
267 c write (iout,*) "Pcalc"
268 c write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
269 Esaxs_all(i)=Esaxs_constr
270 write (iout,*) "Esaxs",Esaxs_constr
273 if (link_start.gt.0) then
274 do j=link_start,link_end
275 if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
276 dxlink=dist(ihpb(j),jhpb(j))
277 if (dxlink.le.25.0d0) then
278 write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
279 & irestr_type(j),ihpb(j),jhpb(j),
282 nviolxlink=nviolxlink+1
283 write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
284 & irestr_type(j),ihpb(j),jhpb(j),
290 & write (iout,*) nviolxlink," crosslink violations."
291 c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
292 c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
293 c & "tmscore",tmscore_closest(i)
295 c Determine # violated NMR restraints
296 if (link_end_peak.gt.0) then
298 write (NUMM,'(bz,i4.4)') i
299 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
300 & //"K_"//NUMM(:ilen(NUMM))//'.nmr'
301 open(jrms,file=cfname)
302 do j=link_start_peak,link_end_peak
304 do ip=ipeak(1,j),ipeak(2,j)
308 c iip=ip-ipeak(1,j)+1
309 C iii and jjj point to the residues for which the distance is assigned.
319 if (dd.lt.dhpb1_peak(ip)) then
321 c write (iout,*) j,iii,jjj,iiib
322 write (jrms,'(4i6)') j,iii,jjj,iiib
326 nviolpeak=nviolpeak+1
327 list_peak_viol(nviolpeak)=j
330 if (nviolpeak.gt.0) then
331 write (iout,'(a,i5,2h (f8.4,2h%))')
332 & "Number of violated NMR restraints:",
333 & nviolpeak,100*(nviolpeak+0.)/npeak
334 write (iout,'(a)')"List of violated restraints:"
335 write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak)
339 if (.not.raw_psipred .and. idihconstr_end.gt.0) then
340 cfname=prefixp(:ilen(prefixp))//"_T"
341 & //ctemper(:ilen(ctemper))
342 & //"K_"//NUMM(:ilen(NUMM))//'.angle'
343 open(jrms,file=cfname)
344 call int_from_cart1(.false.)
346 do j=idihconstr_start,idihconstr_end
349 difi=pinorm(phii-phi0(j))
350 if (difi.gt.drange(j) .or. difi.lt.-drange(j))
351 & nangviol=nangviol+1
352 write (jrms,'(i5,3f10.3)') itori,phii*rad2deg,
353 & phi0(j)*rad2deg,rad2deg*drange(j)
355 write (iout,'(a,i5)')"Number of angle-restraint violations:"
360 call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
361 write (ipdb,'("TER")')
368 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
369 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
370 & //"K_"//'ave'//'.dist'
371 OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
372 write (99,'(5h# ,10f10.5)')
373 & (Esaxs_all(i)*wsaxs,i=1,ngr_print)
375 write (99,'(f6.2,10f10.5)') distsaxs(j),
376 & (Pcalc_all(j,i),i=1,ngr_print)
381 IF (printmol2(icut).gt.0) THEN
382 c Write out a number of conformations from each family in PDB format and
383 c create InsightII command file for their displaying in different colors
388 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
389 write (NUMM,'(bz,i4.4)') i
390 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
391 & //"K_"//numm(:ilen(numm))//extmol
392 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
394 do while (ncon_out.lt.printmol2(icut) .and.
395 & ncon_out.lt.licz(i).and.
396 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
403 c(k,ii)=allcart(k,ii,icon)
406 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
414 call WRITE_STATS(ICUT,NCON,IB)
415 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
416 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
417 & ' CONTAINS ',I4,' CONFORMATION(S): ')
418 c 300 FORMAT ( 8(I4,F6.1))
419 300 FORMAT (5(I4,1pe12.3))
420 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
421 500 FORMAT (8(2I4,2X))
422 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
425 c------------------------------------------------------------------------------
426 subroutine ave_coord(igr)
429 include 'sizesclu.dat'
430 include 'COMMON.CONTROL'
431 include 'COMMON.CLUSTER'
432 include 'COMMON.CHAIN'
433 include 'COMMON.INTERACT'
435 include 'COMMON.TEMPFAC'
436 include 'COMMON.IOUNITS'
438 double precision przes(3),obrot(3,3)
439 double precision xx(3,maxres2),csq(3,maxres2)
440 double precision eref
441 double precision rmscalc
442 c double precision rmscheck
443 integer i,ii,j,k,icon,jcon,igr,ipermmin
444 double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
445 c write (iout,*) "AVE_COORD: igr",igr
448 boltz = dexp(-totfree(jcon)+eref)
452 c(j,i)=allcart(j,i,jcon)*boltz
453 cref1(j,i)=allcart(j,i,jcon)
454 csq(j,i)=allcart(j,i,jcon)**2*boltz
459 c write (iout,*) "k",k," jcon",jcon
462 cwork(j,i)=allcart(j,i,jcon)
465 rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin)
466 c write (iout,*) "rms",rms," ipermmin",ipermmin
468 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),
469 c & (obrot(i,j),j=1,3)
471 c if (rms.lt.0.0) then
472 c print *,'error, rms^2 = ',rms,icon,jcon
475 c if (non_conv) print *,non_conv,icon,jcon
476 boltz=dexp(-totfree(jcon)+eref)
477 qpart = qpart + boltz
480 xx(j,i)=allcart(j,i,jcon)
483 call matvec(cwork,obrot,xx,2*nres)
485 c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
486 c & (allcart(j,i,jcon),j=1,3)
488 cwork(j,i)=cwork(j,i)+przes(j)
489 c(j,i)=c(j,i)+cwork(j,i)*boltz
490 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
497 c rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2
500 c write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms
505 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
507 c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
513 tempfac(1,i)=tempfac(1,i)+csq(j,i)
514 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
516 tempfac(1,i)=dsqrt(tempfac(1,i))
517 tempfac(2,i)=dsqrt(tempfac(2,i))
521 c------------------------------------------------------------------------------
522 subroutine fittest_coord(igr)
525 include 'sizesclu.dat'
526 include 'COMMON.IOUNITS'
527 include 'COMMON.CONTROL'
528 include 'COMMON.CLUSTER'
529 include 'COMMON.CHAIN'
530 include 'COMMON.INTERACT'
532 include 'COMMON.FFIELD'
533 include 'COMMON.TORCNSTR'
535 double precision przes(3),obrot(3,3)
536 double precision xx(3,maxres2),yy(3,maxres2)
537 integer i,ii,j,k,icon,jcon,jconmin,igr
538 double precision rms,rmsmin,cwork(3,maxres2)
539 double precision ehpb,Esaxs_constr,edihcnstr
546 c(j,i)=allcart(j,i,jcon)
549 call int_from_cart1(.false.)
553 if (nsaxs.gt.0) call e_saxs(Esaxs_constr)
555 if (ndih_constr.gt.0) call etor_constr(edihcnstr)
556 rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr
557 c write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb,
558 c & " Edihcnstr",edihcnstr
559 if (rms.lt.rmsmin) then
564 write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
567 c(j,i)=allcart(j,i,jconmin)
572 c------------------------------------------------------------------------------
573 subroutine closest_coord(igr)
576 include 'sizesclu.dat'
577 include 'COMMON.IOUNITS'
578 include 'COMMON.CONTROL'
579 include 'COMMON.CLUSTER'
580 include 'COMMON.CHAIN'
581 include 'COMMON.INTERACT'
584 double precision przes(3),obrot(3,3)
585 integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin
586 double precision rms,rmsmin,cwork(3,maxres2)
587 double precision xx(3,maxres2),yy(3,maxres2)
588 double precision rmscalc
596 yy(j,i)=allcart(j,i,jcon)
599 rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin)
600 c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
601 if (non_conv) print *,non_conv,icon,jcon
602 if (rms.lt.rmsmin) then
607 c write (iout,*) "rmsmin",rmsmin," rms",rms
611 c(j,i)=allcart(j,i,jconmin)
616 c------------------------------------------------------------------------------
620 include 'sizesclu.dat'
621 include 'COMMON.IOUNITS'
622 include 'COMMON.CONTROL'
623 include 'COMMON.CLUSTER'
624 include 'COMMON.CHAIN'
625 include 'COMMON.INTERACT'
627 double precision przes(3)
628 integer i,ii,j,k,icon,jcon,jconmin,igr
632 przes(j)=przes(j)+c(j,i)
636 przes(j)=przes(j)/nres
640 c(j,i)=c(j,i)-przes(j)