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 CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol
20 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/
27 c print *,"calling WRTCLUST",ncon
28 c write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
31 temper=1.0d0/(beta_h(ib)*1.987d-3)
32 if (temper.lt.100.0d0) then
33 write(ctemper,'(f3.0)') temper
35 else if (temper.lt.1000.0) then
36 write (ctemper,'(f4.0)') temper
39 write (ctemper,'(f5.0)') temper
43 do i=1,ncon*(ncon-1)/2
46 close(80,status='delete')
48 C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
50 ii1= index(intinname,'/')
55 ii2=index(intinname(ii1:),'/')
57 ii = ii1+index(intinname(ii1:),'.')-1
63 prefixp=intinname(ii1:ii)
64 cd print *,icut,printang(icut),printpdb(icut),printmol2(icut)
65 cd print *,'ecut=',ecut
68 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
69 NRECORD=LICZ(IGR)/num_in_line
71 DO 63 IRECORD=1,NRECORD
72 IND2=IND1+num_in_line-1
73 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
74 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
77 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),
78 & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
80 ICON=list_conf(NCONF(IGR,1))
81 c WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
82 C 12/8/93 Estimation of "diameters" of the subsequent families.
85 c write (iout,*) "ecut",ecut
86 emin=totfree(nconf(igr,1))
89 if (totfree(ii)-emin .gt. ecut) goto 10
94 ind=ioffset(ncon,ii,jj)
96 ind=ioffset(ncon,jj,ii)
98 c write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
101 curr_dist=dabs(diss(ind)+0.0d0)
102 c write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
103 c & list_conf(jj),curr_dist
104 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
105 ave_dim=ave_dim+curr_dist**2
108 10 if (licz(igr) .gt. 1)
109 & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
110 write (iout,'(/A,F8.1,A,F8.1)')
111 & 'Max. distance in the family:',amax_dim,
112 & '; average distance in the family:',ave_dim
114 gdt_ts_ave(igr)=0.0d0
115 gdt_ha_ave(igr)=0.0d0
116 tmscore_ave(igr)=0.0d0
118 e1=totfree(nconf(igr,1))
121 boltz=dexp(-(totfree(icon)-e1))
122 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
123 gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon)
124 gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon)
125 tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon)
128 rmsave(igr)=rmsave(igr)/qpart
129 gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart
130 gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart
131 tmscore_ave(igr)=tmscore_ave(igr)/qpart
132 write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD",
133 & rmsave(igr)," A, ",
134 & "TMscore",tmscore_ave(igr),
135 & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA",
139 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
140 c print *,icut,printang(icut)
141 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
143 c print *,'emin',emin,' ngr',ngr
144 if (lprint_cart) then
145 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
148 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
153 if (totfree_gr(igr)-emin.le.ecut) then
154 if (lprint_cart) then
155 call cartout(igr,icon,totfree(icon)/beta_h(ib),
156 & totfree_gr(igr)/beta_h(ib),
157 & rmstb(icon),cfname)
159 c print '(a)','calling briefout'
162 c(j,i)=allcart(j,i,icon)
165 call int_from_cart1(.false.)
166 call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),
167 & totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),
168 & jhpb_all(1,icon),cfname)
169 c print '(a)','exit briefout'
175 IF (PRINTPDB(ICUT).gt.0) THEN
176 c Write out a number of conformations from each family in PDB format and
177 c create InsightII command file for their displaying in different colors
178 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
179 & //"K_"//'ave'//exten
180 write (iout,*) "cfname",cfname
181 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
182 write (ipdb,'(a,f8.2)')
183 & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
189 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
190 c write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
191 write (NUMM,'(bz,i4.4)') i
192 ncon_lim=min0(licz(i),printpdb(icut))
193 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
194 & //"K_"//numm(:ilen(numm))//exten
195 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
196 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5,
197 & " AVE RMSD",0pf5.2)')
198 & i,totfree_gr(i)/beta_h(ib),rmsave(i)
199 c Write conformations of the family i to PDB files
201 do while (ncon_out.lt.printpdb(icut) .and.
202 & ncon_out.lt.licz(i).and.
203 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
205 c write (iout,*) i,ncon_out,nconf(i,ncon_out),
206 c & totfree(nconf(i,ncon_out)),emin1,ecut
208 c write (iout,*) "ncon_out",ncon_out
218 c(k,ii)=allcart(k,ii,icon)
222 call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel)
223 write (ipdb,'("TER")')
226 c Average structures and structures closest to average
227 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
228 & //"K_"//'ave'//exten
229 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',
232 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
234 call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
235 write (ipdb,'("TER")')
236 call closest_coord(i)
237 c write (iout,*) "Calling rmsnat"
238 rms_closest(i) = rmsnat(i)
239 call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i),
240 & tmscore_closest(i),cfname,.true.)
241 c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
242 c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
243 c & "tmscore",tmscore_closest(i)
245 call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
246 write (ipdb,'("TER")')
253 IF (printmol2(icut).gt.0) THEN
254 c Write out a number of conformations from each family in PDB format and
255 c create InsightII command file for their displaying in different colors
260 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
261 write (NUMM,'(bz,i4.4)') i
262 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
263 & //"K_"//numm(:ilen(numm))//extmol
264 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
266 do while (ncon_out.lt.printmol2(icut) .and.
267 & ncon_out.lt.licz(i).and.
268 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
275 c(k,ii)=allcart(k,ii,icon)
278 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
286 call WRITE_STATS(ICUT,NCON,IB)
287 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
288 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
289 & ' CONTAINS ',I4,' CONFORMATION(S): ')
290 c 300 FORMAT ( 8(I4,F6.1))
291 300 FORMAT (5(I4,1pe12.3))
292 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
293 500 FORMAT (8(2I4,2X))
294 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
297 c------------------------------------------------------------------------------
298 subroutine ave_coord(igr)
301 include 'sizesclu.dat'
302 include 'COMMON.CONTROL'
303 include 'COMMON.CLUSTER'
304 include 'COMMON.CHAIN'
305 include 'COMMON.INTERACT'
307 include 'COMMON.TEMPFAC'
308 include 'COMMON.IOUNITS'
310 double precision przes(3),obrot(3,3)
311 double precision xx(3,maxres2),yy(3,maxres2),csq(3,maxres2)
312 double precision eref
313 integer i,ii,j,k,icon,jcon,igr
314 double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
315 c write (iout,*) "AVE_COORD: igr",igr
318 boltz = dexp(-totfree(jcon)+eref)
322 c(j,i)=allcart(j,i,jcon)*boltz
323 cref1(j,i)=allcart(j,i,jcon)
324 csq(j,i)=allcart(j,i,jcon)**2*boltz
334 xx(j,ii)=allcart(j,i,jcon)
339 c if (itype(i).ne.10) then
342 xx(j,ii)=allcart(j,i+nres,jcon)
343 yy(j,ii)=cref1(j,i+nres)
347 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
351 cwork(j,i)=allcart(j,i,jcon)
354 call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot
357 c write (iout,*) "rms",rms
359 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
362 print *,'error, rms^2 = ',rms,icon,jcon
365 if (non_conv) print *,non_conv,icon,jcon
366 boltz=dexp(-totfree(jcon)+eref)
367 qpart = qpart + boltz
370 xx(j,i)=allcart(j,i,jcon)
373 call matvec(cwork,obrot,xx,2*nres)
375 c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
376 c & (allcart(j,i,jcon),j=1,3)
378 cwork(j,i)=cwork(j,i)+przes(j)
379 c(j,i)=c(j,i)+cwork(j,i)*boltz
380 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
387 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
389 c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
395 tempfac(1,i)=tempfac(1,i)+csq(j,i)
396 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
398 tempfac(1,i)=dsqrt(tempfac(1,i))
399 tempfac(2,i)=dsqrt(tempfac(2,i))
403 c------------------------------------------------------------------------------
404 subroutine closest_coord(igr)
407 include 'sizesclu.dat'
408 include 'COMMON.IOUNITS'
409 include 'COMMON.CONTROL'
410 include 'COMMON.CLUSTER'
411 include 'COMMON.CHAIN'
412 include 'COMMON.INTERACT'
415 double precision przes(3),obrot(3,3)
416 double precision xx(3,maxres2),yy(3,maxres2)
417 integer i,ii,j,k,icon,jcon,jconmin,igr
418 double precision rms,rmsmin,cwork(3,maxres2)
428 xx(j,ii)=allcart(j,i,jcon)
433 c if (itype(i).ne.10) then
436 xx(j,ii)=allcart(j,i+nres,jcon)
441 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
445 cwork(j,i)=allcart(j,i,jcon)
448 call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot
452 print *,'error, rms^2 = ',rms,icon,jcon
455 c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
456 if (non_conv) print *,non_conv,icon,jcon
457 if (rms.lt.rmsmin) then
462 c write (iout,*) "rmsmin",rmsmin," rms",rms
466 c(j,i)=allcart(j,i,jconmin)
471 c------------------------------------------------------------------------------
475 include 'sizesclu.dat'
476 include 'COMMON.IOUNITS'
477 include 'COMMON.CONTROL'
478 include 'COMMON.CLUSTER'
479 include 'COMMON.CHAIN'
480 include 'COMMON.INTERACT'
482 double precision przes(3)
483 integer i,ii,j,k,icon,jcon,jconmin,igr
487 przes(j)=przes(j)+c(j,i)
491 przes(j)=przes(j)/nres
495 c(j,i)=c(j,i)-przes(j)