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 CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
22 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
25 double precision Esaxs_all(maxgr),Pcalc_all(maxsaxs,maxgr)
30 c print *,"calling WRTCLUST",ncon
31 c write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
34 temper=1.0d0/(beta_h(ib)*1.987d-3)
35 if (temper.lt.100.0d0) then
36 write(ctemper,'(f3.0)') temper
38 else if (temper.lt.1000.0) then
39 write (ctemper,'(f4.0)') temper
42 write (ctemper,'(f5.0)') temper
46 do i=1,ncon*(ncon-1)/2
49 close(80,status='delete')
51 C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
53 ii1= index(intinname,'/')
58 ii2=index(intinname(ii1:),'/')
60 ii = ii1+index(intinname(ii1:),'.')-1
66 prefixp=intinname(ii1:ii)
67 cd print *,icut,printang(icut),printpdb(icut),printmol2(icut)
68 cd print *,'ecut=',ecut
71 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
72 NRECORD=LICZ(IGR)/num_in_line
74 DO 63 IRECORD=1,NRECORD
75 IND2=IND1+num_in_line-1
76 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
77 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
80 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
81 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
83 ICON=list_conf(NCONF(IGR,1))
84 c WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
85 C 12/8/93 Estimation of "diameters" of the subsequent families.
88 c write (iout,*) "ecut",ecut
91 if (totfree(ii)-emin .gt. ecut) goto 10
96 ind=ioffset(ncon,ii,jj)
98 ind=ioffset(ncon,jj,ii)
100 c write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
103 curr_dist=dabs(diss(ind)+0.0d0)
104 c write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
105 c & list_conf(jj),curr_dist
106 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
107 ave_dim=ave_dim+curr_dist**2
110 10 if (licz(igr) .gt. 1)
111 & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
112 write (iout,'(/A,F8.1,A,F8.1)')
113 & 'Max. distance in the family:',amax_dim,
114 & '; average distance in the family:',ave_dim
116 gdt_ts_ave(igr)=0.0d0
117 gdt_ha_ave(igr)=0.0d0
118 tmscore_ave(igr)=0.0d0
120 e1=totfree(nconf(igr,1))
123 boltz=dexp(-(totfree(icon)-e1))
124 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
125 gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
126 gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
127 tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
129 c write (iout,'(2i5,10f10.5)') i,icon,boltz,rmstb(icon),
130 c & gdt_ts_tb(icon),gdt_ha_tb(icon),tmscore_tb(icon)
132 c write (iout,*) "qpart",qpart
133 rmsave(igr)=rmsave(igr)/qpart
134 gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
135 gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
136 tmscore_ave(igr)=tmscore_ave(igr)/qpart
137 write (iout,'(a,f5.2,a,3(a,f6.4))') "Average RMSD",
139 & "average TMscore",tmscore_ave(igr),
140 & " average GDT_TS",gdt_ts_ave(igr)," average GDT_HA",
144 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
145 c print *,icut,printang(icut)
146 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
148 c print *,'emin',emin,' ngr',ngr
149 if (lprint_cart) then
150 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
153 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
158 if (totfree_gr(igr)-emin.le.ecut) then
159 if (lprint_cart) then
160 call cartout(igr,icon,totfree(icon)/beta_h(ib),
161 & totfree_gr(igr)/beta_h(ib),
162 & rmstb(icon),cfname)
164 c print '(a)','calling briefout'
167 c(j,i)=allcart(j,i,icon)
170 call int_from_cart1(.false.)
171 call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
172 & totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
173 & jhpb_all(1,icon),cfname)
174 c print '(a)','exit briefout'
180 IF (PRINTPDB(ICUT).gt.0) THEN
181 c Write out a number of conformations from each family in PDB format and
182 c create InsightII command file for their displaying in different colors
183 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
184 & //"K_"//'ave'//exten
185 write (iout,*) "cfname",cfname
186 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
187 write (ipdb,'(a,f8.2)')
188 & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
194 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
195 c write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
196 write (NUMM,'(bz,i4.4)') i
197 ncon_lim=min0(licz(i),printpdb(icut))
198 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
199 & //"K_"//numm(:ilen(numm))//exten
200 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
201 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
202 & " AVE RMSD",0pf5.2)')
203 & i,totfree_gr(i)/beta_h(ib),rmsave(i)
204 c Write conformations of the family i to PDB files
206 do while (ncon_out.lt.printpdb(icut) .and.
207 & ncon_out.lt.licz(i).and.
208 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
210 c write (iout,*) i,ncon_out,nconf(i,ncon_out),
211 c & totfree(nconf(i,ncon_out)),emin1,ecut
213 c write (iout,*) "ncon_out",ncon_out
223 c(k,ii)=allcart(k,ii,icon)
227 call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
228 write (ipdb,'("TER")')
231 c Average structures and structures closest to average
232 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
233 & //"K_"//'ave'//exten
234 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
237 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
239 call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
240 write (ipdb,'("TER")')
241 if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0)) then
242 call fittest_coord(i)
244 call closest_coord(i)
246 c write (iout,*) "Calling rmsnat"
247 rms_closest(i) = rmsnat(i)
249 write (iout,*) "Cluster",i
250 call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
251 & tmscore_closest(i),cfname,.true.)
252 c write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i
253 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
254 call e_saxs(Esaxs_constr)
257 Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
258 & (Pcalc(j+1)+Pcalc(j))/2
261 Pcalc_all(j,i)=Pcalc(j)/Cnorm
263 c write (iout,*) "Pcalc"
264 c write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
265 Esaxs_all(i)=Esaxs_constr
266 write (iout,*) "Esaxs",Esaxs_constr
269 if (link_start.gt.0) then
270 do j=link_start,link_end
271 if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
272 dxlink=dist(ihpb(j),jhpb(j))
273 if (dxlink.le.25.0d0) then
274 write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
275 & irestr_type(j),ihpb(j),jhpb(j),
278 nviolxlink=nviolxlink+1
279 write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
280 & irestr_type(j),ihpb(j),jhpb(j),
286 & write (iout,*) nviolxlink," crosslink violations."
287 c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
288 c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
289 c & "tmscore",tmscore_closest(i)
292 call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
293 write (ipdb,'("TER")')
300 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
301 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
302 & //"K_"//'ave'//'.dist'
303 OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
304 write (99,'(5h# ,10f10.5)')
305 & (Esaxs_all(i)*wsaxs,i=1,ngr_print)
307 write (99,'(f6.2,10f10.5)') distsaxs(j),
308 & (Pcalc_all(j,i),i=1,ngr_print)
313 IF (printmol2(icut).gt.0) THEN
314 c Write out a number of conformations from each family in PDB format and
315 c create InsightII command file for their displaying in different colors
320 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
321 write (NUMM,'(bz,i4.4)') i
322 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
323 & //"K_"//numm(:ilen(numm))//extmol
324 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
326 do while (ncon_out.lt.printmol2(icut) .and.
327 & ncon_out.lt.licz(i).and.
328 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
335 c(k,ii)=allcart(k,ii,icon)
338 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
346 call WRITE_STATS(ICUT,NCON,IB)
347 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
348 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
349 & ' CONTAINS ',I4,' CONFORMATION(S): ')
350 c 300 FORMAT ( 8(I4,F6.1))
351 300 FORMAT (5(I4,1pe12.3))
352 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
353 500 FORMAT (8(2I4,2X))
354 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
357 c------------------------------------------------------------------------------
358 subroutine ave_coord(igr)
361 include 'sizesclu.dat'
362 include 'COMMON.CONTROL'
363 include 'COMMON.CLUSTER'
364 include 'COMMON.CHAIN'
365 include 'COMMON.INTERACT'
367 include 'COMMON.TEMPFAC'
368 include 'COMMON.IOUNITS'
370 double precision przes(3),obrot(3,3)
371 double precision xx(3,maxres2),yy(3,maxres2),csq(3,maxres2)
372 double precision eref
373 integer i,ii,j,k,icon,jcon,igr
374 double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
375 c write (iout,*) "AVE_COORD: igr",igr
378 boltz = dexp(-totfree(jcon)+eref)
382 c(j,i)=allcart(j,i,jcon)*boltz
383 cref1(j,i)=allcart(j,i,jcon)
384 csq(j,i)=allcart(j,i,jcon)**2*boltz
394 xx(j,ii)=allcart(j,i,jcon)
399 c if (itype(i).ne.10) then
402 xx(j,ii)=allcart(j,i+nres,jcon)
403 yy(j,ii)=cref1(j,i+nres)
407 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
411 cwork(j,i)=allcart(j,i,jcon)
414 call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot
417 c write (iout,*) "rms",rms
419 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
422 print *,'error, rms^2 = ',rms,icon,jcon
425 if (non_conv) print *,non_conv,icon,jcon
426 boltz=dexp(-totfree(jcon)+eref)
427 qpart = qpart + boltz
430 xx(j,i)=allcart(j,i,jcon)
433 call matvec(cwork,obrot,xx,2*nres)
435 c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
436 c & (allcart(j,i,jcon),j=1,3)
438 cwork(j,i)=cwork(j,i)+przes(j)
439 c(j,i)=c(j,i)+cwork(j,i)*boltz
440 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
447 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
449 c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
455 tempfac(1,i)=tempfac(1,i)+csq(j,i)
456 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
458 tempfac(1,i)=dsqrt(tempfac(1,i))
459 tempfac(2,i)=dsqrt(tempfac(2,i))
463 c------------------------------------------------------------------------------
464 subroutine fittest_coord(igr)
467 include 'sizesclu.dat'
468 include 'COMMON.IOUNITS'
469 include 'COMMON.CONTROL'
470 include 'COMMON.CLUSTER'
471 include 'COMMON.CHAIN'
472 include 'COMMON.INTERACT'
474 include 'COMMON.FFIELD'
476 double precision przes(3),obrot(3,3)
477 double precision xx(3,maxres2),yy(3,maxres2)
478 integer i,ii,j,k,icon,jcon,jconmin,igr
479 double precision rms,rmsmin,cwork(3,maxres2)
480 double precision ehpb,Esaxs_constr
487 c(j,i)=allcart(j,i,jcon)
490 call e_saxs(Esaxs_constr)
492 rms=wsaxs*esaxs_constr+wstrain*ehpb
493 if (rms.lt.rmsmin) then
498 write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
501 c(j,i)=allcart(j,i,jconmin)
506 c------------------------------------------------------------------------------
507 subroutine closest_coord(igr)
510 include 'sizesclu.dat'
511 include 'COMMON.IOUNITS'
512 include 'COMMON.CONTROL'
513 include 'COMMON.CLUSTER'
514 include 'COMMON.CHAIN'
515 include 'COMMON.INTERACT'
518 double precision przes(3),obrot(3,3)
519 double precision xx(3,maxres2),yy(3,maxres2)
520 integer i,ii,j,k,icon,jcon,jconmin,igr
521 double precision rms,rmsmin,cwork(3,maxres2)
531 xx(j,ii)=allcart(j,i,jcon)
536 c if (itype(i).ne.10) then
539 xx(j,ii)=allcart(j,i+nres,jcon)
544 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
548 cwork(j,i)=allcart(j,i,jcon)
551 call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot
555 print *,'error, rms^2 = ',rms,icon,jcon
558 c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
559 if (non_conv) print *,non_conv,icon,jcon
560 if (rms.lt.rmsmin) then
565 c write (iout,*) "rmsmin",rmsmin," rms",rms
569 c(j,i)=allcart(j,i,jconmin)
574 c------------------------------------------------------------------------------
578 include 'sizesclu.dat'
579 include 'COMMON.IOUNITS'
580 include 'COMMON.CONTROL'
581 include 'COMMON.CLUSTER'
582 include 'COMMON.CHAIN'
583 include 'COMMON.INTERACT'
585 double precision przes(3)
586 integer i,ii,j,k,icon,jcon,jconmin,igr
590 przes(j)=przes(j)+c(j,i)
594 przes(j)=przes(j)/nres
598 c(j,i)=c(j,i)-przes(j)