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 c 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)
252 c write (iout,*) "Calling rmsnat"
253 rms_closest(i) = rmsnat(i)
254 c 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
259 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
260 call e_saxs(Esaxs_constr)
263 Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))*
264 & (Pcalc(j+1)+Pcalc(j))/2
267 Pcalc_all(j,i)=Pcalc(j)/Cnorm
269 c write (iout,*) "Pcalc"
270 c write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs)
271 Esaxs_all(i)=Esaxs_constr
272 write (iout,*) "Esaxs",Esaxs_constr
275 if (link_start.gt.0) then
276 do j=link_start,link_end
277 if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then
278 dxlink=dist(ihpb(j),jhpb(j))
279 if (dxlink.le.25.0d0) then
280 write (iout,'(a,i2,2i5,f8.2)') "XLINK-",
281 & irestr_type(j),ihpb(j),jhpb(j),
284 nviolxlink=nviolxlink+1
285 write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-",
286 & irestr_type(j),ihpb(j),jhpb(j),
292 & write (iout,*) nviolxlink," crosslink violations."
293 c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts",
294 c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i),
295 c & "tmscore",tmscore_closest(i)
297 c Determine # violated NMR restraints
298 if (link_end_peak.gt.0) then
300 write (NUMM,'(bz,i4.4)') i
301 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
302 & //"K_"//NUMM(:ilen(NUMM))//'.nmr'
303 open(jrms,file=cfname)
304 do j=link_start_peak,link_end_peak
306 do ip=ipeak(1,j),ipeak(2,j)
310 c iip=ip-ipeak(1,j)+1
311 C iii and jjj point to the residues for which the distance is assigned.
321 if (dd.lt.dhpb1_peak(ip)) then
323 c write (iout,*) j,iii,jjj,iiib
324 write (jrms,'(4i6)') j,iii,jjj,iiib
328 nviolpeak=nviolpeak+1
329 list_peak_viol(nviolpeak)=j
332 if (nviolpeak.gt.0) then
333 write (iout,'(a,i5,2h (f8.4,2h%))')
334 & "Number of violated NMR restraints:",
335 & nviolpeak,100*(nviolpeak+0.)/npeak
336 write (iout,'(a)')"List of violated restraints:"
337 write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak)
341 if (.not.raw_psipred .and. idihconstr_end.gt.0) then
342 cfname=prefixp(:ilen(prefixp))//"_T"
343 & //ctemper(:ilen(ctemper))
344 & //"K_"//NUMM(:ilen(NUMM))//'.angle'
345 open(jrms,file=cfname)
346 call int_from_cart1(.false.)
348 do j=idihconstr_start,idihconstr_end
351 difi=pinorm(phii-phi0(j))
352 if (difi.gt.drange(j) .or. difi.lt.-drange(j))
353 & nangviol=nangviol+1
354 write (jrms,'(i5,3f10.3)') itori,phii*rad2deg,
355 & phi0(j)*rad2deg,rad2deg*drange(j)
357 write (iout,'(a,i5)')"Number of angle-restraint violations:"
362 call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel)
363 write (ipdb,'("TER")')
370 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then
371 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
372 & //"K_"//'ave'//'.dist'
373 OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
374 write (99,'(5h# ,10f10.5)')
375 & (Esaxs_all(i)*wsaxs,i=1,ngr_print)
377 write (99,'(f6.2,10f10.5)') distsaxs(j),
378 & (Pcalc_all(j,i),i=1,ngr_print)
383 IF (printmol2(icut).gt.0) THEN
384 c Write out a number of conformations from each family in PDB format and
385 c create InsightII command file for their displaying in different colors
390 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
391 write (NUMM,'(bz,i4.4)') i
392 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))
393 & //"K_"//numm(:ilen(numm))//extmol
394 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
396 do while (ncon_out.lt.printmol2(icut) .and.
397 & ncon_out.lt.licz(i).and.
398 & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
405 c(k,ii)=allcart(k,ii,icon)
408 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
416 call WRITE_STATS(ICUT,NCON,IB)
417 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
418 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,
419 & ' CONTAINS ',I4,' CONFORMATION(S): ')
420 c 300 FORMAT ( 8(I4,F6.1))
421 300 FORMAT (5(I4,1pe12.3))
422 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
423 500 FORMAT (8(2I4,2X))
424 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
427 c------------------------------------------------------------------------------
428 subroutine ave_coord(igr)
431 include 'sizesclu.dat'
432 include 'COMMON.CONTROL'
433 include 'COMMON.CLUSTER'
434 include 'COMMON.CHAIN'
435 include 'COMMON.INTERACT'
437 include 'COMMON.TEMPFAC'
438 include 'COMMON.IOUNITS'
440 double precision przes(3),obrot(3,3)
441 double precision xx(3,maxres2),csq(3,maxres2)
442 double precision eref
443 double precision rmscalc
444 c double precision rmscheck
445 integer i,ii,j,k,icon,jcon,igr,ipermmin
446 double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2)
447 c write (iout,*) "AVE_COORD: igr",igr
450 boltz = dexp(-totfree(jcon)+eref)
454 c(j,i)=allcart(j,i,jcon)*boltz
455 cref1(j,i)=allcart(j,i,jcon)
456 csq(j,i)=allcart(j,i,jcon)**2*boltz
461 c write (iout,*) "k",k," jcon",jcon
464 cwork(j,i)=allcart(j,i,jcon)
467 rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin)
468 c write (iout,*) "rms",rms," ipermmin",ipermmin
470 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),
471 c & (obrot(i,j),j=1,3)
473 c if (rms.lt.0.0) then
474 c print *,'error, rms^2 = ',rms,icon,jcon
477 c if (non_conv) print *,non_conv,icon,jcon
478 boltz=dexp(-totfree(jcon)+eref)
479 qpart = qpart + boltz
482 xx(j,i)=allcart(j,i,jcon)
485 call matvec(cwork,obrot,xx,2*nres)
487 c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
488 c & (allcart(j,i,jcon),j=1,3)
490 cwork(j,i)=cwork(j,i)+przes(j)
491 c(j,i)=c(j,i)+cwork(j,i)*boltz
492 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
499 c rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2
502 c write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms
507 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
509 c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
515 tempfac(1,i)=tempfac(1,i)+csq(j,i)
516 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
518 tempfac(1,i)=dsqrt(tempfac(1,i))
519 tempfac(2,i)=dsqrt(tempfac(2,i))
523 c------------------------------------------------------------------------------
524 subroutine fittest_coord(igr)
527 include 'sizesclu.dat'
528 include 'COMMON.IOUNITS'
529 include 'COMMON.CONTROL'
530 include 'COMMON.CLUSTER'
531 include 'COMMON.CHAIN'
532 include 'COMMON.INTERACT'
534 include 'COMMON.FFIELD'
535 include 'COMMON.TORCNSTR'
536 include 'COMMON.SAXS'
538 double precision przes(3),obrot(3,3)
539 double precision xx(3,maxres2),yy(3,maxres2)
540 integer i,ii,j,k,icon,jcon,jconmin,igr
541 double precision rms,rmsmin,cwork(3,maxres2)
542 double precision ehpb,Esaxs_constr,edihcnstr
549 c(j,i)=allcart(j,i,jcon)
552 call int_from_cart1(.false.)
556 if (nsaxs.gt.0) call e_saxs(Esaxs_constr)
558 if (ndih_constr.gt.0) call etor_constr(edihcnstr)
559 rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr
560 c write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb,
561 c & " Edihcnstr",edihcnstr
562 if (rms.lt.rmsmin) then
567 write (iout,*) "fittest conformation",jconmin," penalty",rmsmin
570 c(j,i)=allcart(j,i,jconmin)
575 c------------------------------------------------------------------------------
576 subroutine closest_coord(igr)
579 include 'sizesclu.dat'
580 include 'COMMON.IOUNITS'
581 include 'COMMON.CONTROL'
582 include 'COMMON.CLUSTER'
583 include 'COMMON.CHAIN'
584 include 'COMMON.INTERACT'
587 double precision przes(3),obrot(3,3)
588 integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin
589 double precision rms,rmsmin,cwork(3,maxres2)
590 double precision xx(3,maxres2),yy(3,maxres2)
591 double precision rmscalc
594 c write (iout,*) "CLOSEST_COORD: Average coords"
601 yy(j,i)=allcart(j,i,jcon)
604 rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin)
605 c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
606 if (non_conv) print *,non_conv,icon,jcon
607 if (rms.lt.rmsmin) then
612 c write (iout,*) "rmsmin",rmsmin," rms",rms
616 c(j,i)=allcart(j,i,jconmin)
621 c------------------------------------------------------------------------------
625 include 'sizesclu.dat'
626 include 'COMMON.IOUNITS'
627 include 'COMMON.CONTROL'
628 include 'COMMON.CLUSTER'
629 include 'COMMON.CHAIN'
630 include 'COMMON.INTERACT'
632 double precision przes(3)
633 integer i,ii,j,k,icon,jcon,jconmin,igr
637 przes(j)=przes(j)+c(j,i)
641 przes(j)=przes(j)/nres
645 c(j,i)=c(j,i)-przes(j)