SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib) implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'sizesclu.dat' parameter (num_in_line=5) LOGICAL PRINTANG(max_cut) integer PRINTPDB(max_cut),printmol2(max_cut) include 'COMMON.CONTROL' include 'COMMON.HEADER' include 'COMMON.CHAIN' include 'COMMON.VAR' include 'COMMON.CLUSTER' include 'COMMON.IOUNITS' include 'COMMON.GEO' include 'COMMON.FREE' include 'COMMON.TEMPFAC' include 'COMMON.FFIELD' include 'COMMON.SBRIDGE' include 'COMMON.TORCNSTR' include 'COMMON.SAXS' CHARACTER*64 prefixp,NUMM,MUMM,EXTEN,extmol character*120 cfname character*8 ctemper DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,MUMM /'000'/ external ilen logical viol_nmr integer ib,list_peak_viol(maxdim_cont) double precision Esaxs_all(maxgr),Pcalc_all(maxsaxs,maxgr) do i=1,64 cfname(i:i)=" " enddo c print *,"calling WRTCLUST",ncon c write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut) rewind 80 call flush(iout) temper=1.0d0/(beta_h(ib)*1.987d-3) if (temper.lt.100.0d0) then write(ctemper,'(f3.0)') temper ctemper(3:3)=" " else if (temper.lt.1000.0) then write (ctemper,'(f4.0)') temper ctemper(4:4)=" " else write (ctemper,'(f5.0)') temper ctemper(5:5)=" " endif do i=1,ncon*(ncon-1)/2 read (80) diss(i) enddo close(80,status='delete') C C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS C ii1= index(intinname,'/') ii2=ii1 ii1=ii1+1 do while (ii2.gt.0) ii1=ii1+ii2 ii2=index(intinname(ii1:),'/') enddo ii = ii1+index(intinname(ii1:),'.')-1 if (ii.eq.0) then ii=ilen(intinname) else ii=ii-1 endif prefixp=intinname(ii1:ii) cd print *,icut,printang(icut),printpdb(icut),printmol2(icut) cd print *,'ecut=',ecut WRITE (iout,100) NGR DO 19 IGR=1,NGR WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR) NRECORD=LICZ(IGR)/num_in_line IND1=1 DO 63 IRECORD=1,NRECORD IND2=IND1+num_in_line-1 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)), & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2) IND1=IND2+1 63 CONTINUE WRITE (iout,300) (list_conf(NCONF(IGR,ICO)), & totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR)) IND1=1 ICON=list_conf(NCONF(IGR,1)) c WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3) C 12/8/93 Estimation of "diameters" of the subsequent families. ave_dim=0.0 amax_dim=0.0 c write (iout,*) "ecut",ecut emin=totfree(nconf(igr,1)) do i=2,licz(igr) ii=nconf(igr,i) if (totfree(ii)-emin .gt. ecut) goto 10 do j=1,i-1 jj=nconf(igr,j) if (jj.eq.1) exit if (ii.lt.jj) then ind=ioffset(ncon,ii,jj) else ind=ioffset(ncon,jj,ii) endif c write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj, c & " ind",ind," diss",diss(ind) c call flush(iout) curr_dist=dabs(diss(ind)+0.0d0) c write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii), c & list_conf(jj),curr_dist if (curr_dist .gt. amax_dim) amax_dim=curr_dist ave_dim=ave_dim+curr_dist**2 enddo enddo 10 if (licz(igr) .gt. 1) & ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2)) write (iout,'(/A,F8.1,A,F8.1)') & 'Max. distance in the family:',amax_dim, & '; average distance in the family:',ave_dim rmsave(igr)=0.0d0 gdt_ts_ave(igr)=0.0d0 gdt_ha_ave(igr)=0.0d0 tmscore_ave(igr)=0.0d0 qpart=0.0d0 e1=totfree(nconf(igr,1)) do i=1,licz(igr) icon=nconf(igr,i) boltz=dexp(-(totfree(icon)-e1)) rmsave(igr)=rmsave(igr)+boltz*rmstb(icon) gdt_ts_ave(igr)=gdt_ts_ave(igr)+boltz*gdt_ts_tb(icon) gdt_ha_ave(igr)=gdt_ha_ave(igr)+boltz*gdt_ha_tb(icon) tmscore_ave(igr)=tmscore_ave(igr)+boltz*tmscore_tb(icon) qpart=qpart+boltz c write (iout,'(2i5,10f10.5)') i,icon,boltz,rmstb(icon), c & gdt_ts_tb(icon),gdt_ha_tb(icon),tmscore_tb(icon) enddo c write (iout,*) "qpart",qpart rmsave(igr)=rmsave(igr)/qpart gdt_ts_ave(igr)=gdt_ts_ave(igr)/qpart gdt_ha_ave(igr)=gdt_ha_ave(igr)/qpart tmscore_ave(igr)=tmscore_ave(igr)/qpart write (iout,'(a,f5.2,a,3(a,f7.4))') "Cluster averages: RMSD", & rmsave(igr)," A, ", & "TMscore",tmscore_ave(igr), & ", GDT_TS",gdt_ts_ave(igr),", GDT_HA", & gdt_ha_ave(igr) 19 CONTINUE WRITE (iout,400) WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON) c print *,icut,printang(icut) IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then emin=totfree_gr(1) c print *,'emin',emin,' ngr',ngr if (lprint_cart) then cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K"//".x" else cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K"//".int" endif do igr=1,ngr icon=nconf(igr,1) if (totfree_gr(igr)-emin.le.ecut) then if (lprint_cart) then call cartout(igr,icon,totfree(icon)/beta_h(ib), & totfree_gr(igr)/beta_h(ib), & rmstb(icon),cfname) else c print '(a)','calling briefout' do i=1,2*nres do j=1,3 c(j,i)=allcart(j,i,icon) enddo enddo call int_from_cart1(.false.) call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib), & totfree_gr(igr),nss_all(icon),ihpb_all(1,icon), & jhpb_all(1,icon),cfname) c print '(a)','exit briefout' endif endif enddo close(igeom) ENDIF IF (PRINTPDB(ICUT).gt.0) THEN c Write out a number of conformations from each family in PDB format and c create InsightII command file for their displaying in different colors cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//'ave'//exten write (iout,*) "cfname",cfname OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED') write (ipdb,'(a,f8.2)') & "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper close (ipdb) I=1 ICON=NCONF(1,1) EMIN=totfree_gr(I) emin1=totfree(icon) DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT) c write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut write (NUMM,'(bz,i4.4)') i ncon_lim=min0(licz(i),printpdb(icut)) cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//numm(:ilen(numm))//exten OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED') write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5, & " AVE RMSD",0pf5.2)') & i,totfree_gr(i)/beta_h(ib),rmsave(i) c Write conformations of the family i to PDB files ncon_out=1 do while (ncon_out.lt.printpdb(icut) .and. & ncon_out.lt.licz(i).and. & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT) ncon_out=ncon_out+1 c write (iout,*) i,ncon_out,nconf(i,ncon_out), c & totfree(nconf(i,ncon_out)),emin1,ecut enddo c write (iout,*) "ncon_out",ncon_out call flush(iout) do j=1,nres tempfac(1,j)=5.0d0 tempfac(2,j)=5.0d0 enddo do j=1,ncon_out icon=nconf(i,j) do ii=1,2*nres do k=1,3 c(k,ii)=allcart(k,ii,icon) enddo enddo call center call pdbout(totfree(icon)/beta_h(ib),rmstb(icon),titel) write (ipdb,'("TER")') enddo close(ipdb) c Average structures and structures closest to average cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//'ave'//exten OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED', & position="APPEND") call ave_coord(i) write (ipdb,'(a,i5)') "REMARK CLUSTER",i call center call pdbout(totfree_gr(i)/beta_h(ib),rmsave(i),titel) write (ipdb,'("TER")') if (print_fittest.and.(nsaxs.gt.0 .or. nhpb.gt.0 & .or.npeak.gt.0)) then call fittest_coord(i) else call closest_coord(i) endif c write (iout,*) "Calling rmsnat" rms_closest(i) = rmsnat(i) write (iout,*) "Cluster",i call TMscore_sub(rmsd,gdt_ts_closest(i),gdt_ha_closest(i), & tmscore_closest(i),cfname,.true.) c write (iout,*) "WRTCLUST: nsaxs",nsaxs," i",i if (nsaxs.gt.0 .and. saxs_mode.eq.0) then call e_saxs(Esaxs_constr) Cnorm=0.0d0 do j=1,nsaxs-1 Cnorm=Cnorm+(distsaxs(j+1)-distsaxs(j))* & (Pcalc(j+1)+Pcalc(j))/2 enddo do j=1,nsaxs Pcalc_all(j,i)=Pcalc(j)/Cnorm enddo c write (iout,*) "Pcalc" c write (iout,'(f6.2,f10.5)') (distsaxs(j),Pcalc(j),j=1,nsaxs) Esaxs_all(i)=Esaxs_constr write (iout,*) "Esaxs",Esaxs_constr endif nviolxlink=0 if (link_start.gt.0) then do j=link_start,link_end if (irestr_type(j).eq.10 .or. irestr_type(j).eq. 11) then dxlink=dist(ihpb(j),jhpb(j)) if (dxlink.le.25.0d0) then write (iout,'(a,i2,2i5,f8.2)') "XLINK-", & irestr_type(j),ihpb(j),jhpb(j), & dxlink else nviolxlink=nviolxlink+1 write (iout,'(a,i2,2i5,f8.2,2h *)') "XLINK-", & irestr_type(j),ihpb(j),jhpb(j), & dxlink endif endif enddo if (nviolxlink.gt.0) & write (iout,*) nviolxlink," crosslink violations." c write (iout,*) "Family",i," rmsd",rmsd,"gdt_ts", c & gdt_ts_closest(i)," gdt_ha",gdt_ha_closest(i), c & "tmscore",tmscore_closest(i) endif c Determine # violated NMR restraints if (link_end_peak.gt.0) then nviolpeak=0 write (NUMM,'(bz,i4.4)') i cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//NUMM(:ilen(NUMM))//'.nmr' open(jrms,file=cfname) do j=link_start_peak,link_end_peak viol_nmr=.true. do ip=ipeak(1,j),ipeak(2,j) ii=ihpb_peak(ip) jj=jhpb_peak(ip) dd=dist(ii,jj) c iip=ip-ipeak(1,j)+1 C iii and jjj point to the residues for which the distance is assigned. if (ii.gt.nres) then iii=ii-nres jjj=jj-nres iiib=1 else iii=ii jjj=jj iiib=0 endif if (dd.lt.dhpb1_peak(ip)) then viol_nmr=.false. c write (iout,*) j,iii,jjj,iiib write (jrms,'(4i6)') j,iii,jjj,iiib endif enddo if (viol_nmr) then nviolpeak=nviolpeak+1 list_peak_viol(nviolpeak)=j endif enddo if (nviolpeak.gt.0) then write (iout,'(a,i5,2h (f8.4,2h%))') & "Number of violated NMR restraints:", & nviolpeak,100*(nviolpeak+0.)/npeak write (iout,'(a)')"List of violated restraints:" write (iout,'(16i5)') (list_peak_viol(j),j=1,nviolpeak) endif close(jrms) endif if (.not.raw_psipred .and. idihconstr_end.gt.0) then cfname=prefixp(:ilen(prefixp))//"_T" & //ctemper(:ilen(ctemper)) & //"K_"//NUMM(:ilen(NUMM))//'.angle' open(jrms,file=cfname) call int_from_cart1(.false.) nangviol=0 do j=idihconstr_start,idihconstr_end itori=idih_constr(j) phii=phi(itori) difi=pinorm(phii-phi0(j)) if (difi.gt.drange(j) .or. difi.lt.-drange(j)) & nangviol=nangviol+1 write (jrms,'(i5,3f10.3)') itori,phii*rad2deg, & phi0(j)*rad2deg,rad2deg*drange(j) enddo write (iout,'(a,i5)')"Number of angle-restraint violations:" & ,nangviol close(jrms) endif call center call pdbout(totfree_gr(i)/beta_h(ib),rms_closest(i),titel) write (ipdb,'("TER")') close (ipdb) I=I+1 ICON=NCONF(I,1) emin1=totfree(icon) ENDDO ngr_print=i-1 if (nsaxs.gt.0 .and. saxs_mode.eq.0) then cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//'ave'//'.dist' OPEN(99,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED') write (99,'(5h# ,10f10.5)') & (Esaxs_all(i)*wsaxs,i=1,ngr_print) do j=1,nsaxs write (99,'(f6.2,10f10.5)') distsaxs(j), & (Pcalc_all(j,i),i=1,ngr_print) enddo close(99) endif ENDIF IF (printmol2(icut).gt.0) THEN c Write out a number of conformations from each family in PDB format and c create InsightII command file for their displaying in different colors I=1 ICON=NCONF(1,1) EMIN=ENERGY(ICON) emin1=emin DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT) write (NUMM,'(bz,i4.4)') i cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper)) & //"K_"//numm(:ilen(numm))//extmol OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED') ncon_out=1 do while (ncon_out.lt.printmol2(icut) .and. & ncon_out.lt.licz(i).and. & totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT) ncon_out=ncon_out+1 enddo do j=1,ncon_out icon=nconf(i,j) do ii=1,2*nres do k=1,3 c(k,ii)=allcart(k,ii,icon) enddo enddo CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm) enddo CLOSE(imol2) I=I+1 ICON=NCONF(I,1) emin1=totfree(icon) ENDDO ENDIF call WRITE_STATS(ICUT,NCON,IB) 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS') 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5, & ' CONTAINS ',I4,' CONFORMATION(S): ') c 300 FORMAT ( 8(I4,F6.1)) 300 FORMAT (5(I4,1pe12.3)) 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:') 500 FORMAT (8(2I4,2X)) 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6) RETURN END c------------------------------------------------------------------------------ subroutine ave_coord(igr) implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.VAR' include 'COMMON.TEMPFAC' include 'COMMON.IOUNITS' logical non_conv double precision przes(3),obrot(3,3) double precision xx(3,maxres2),csq(3,maxres2) double precision eref double precision rmscalc c double precision rmscheck integer i,ii,j,k,icon,jcon,igr,ipermmin double precision rms,boltz,qpart,cwork(3,maxres2),cref1(3,maxres2) c write (iout,*) "AVE_COORD: igr",igr jcon=nconf(igr,1) eref=totfree(jcon) boltz = dexp(-totfree(jcon)+eref) qpart=boltz do i=1,2*nres do j=1,3 c(j,i)=allcart(j,i,jcon)*boltz cref1(j,i)=allcart(j,i,jcon) csq(j,i)=allcart(j,i,jcon)**2*boltz enddo enddo DO K=2,LICZ(IGR) jcon=nconf(igr,k) c write (iout,*) "k",k," jcon",jcon do i=1,2*nres do j=1,3 cwork(j,i)=allcart(j,i,jcon) enddo enddo rms=rmscalc(cwork(1,1),cref1(1,1),przes,obrot,ipermmin) c write (iout,*) "rms",rms," ipermmin",ipermmin c do i=1,3 c write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i), c & (obrot(i,j),j=1,3) c enddo c if (rms.lt.0.0) then c print *,'error, rms^2 = ',rms,icon,jcon c stop c endif c if (non_conv) print *,non_conv,icon,jcon boltz=dexp(-totfree(jcon)+eref) qpart = qpart + boltz do i=1,2*nres do j=1,3 xx(j,i)=allcart(j,i,jcon) enddo enddo call matvec(cwork,obrot,xx,2*nres) do i=1,2*nres c write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3), c & (allcart(j,i,jcon),j=1,3) do j=1,3 cwork(j,i)=cwork(j,i)+przes(j) c(j,i)=c(j,i)+cwork(j,i)*boltz csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz enddo enddo c rms check c rmscheck=0.0d0 c do i=nnt,nct c do j=1,3 c rmscheck=rmscheck+(cwork(j,i)-cref1(j,i))**2 c enddo c enddo c write (iout,*) "rmscheck",dsqrt(rmscheck/(nct-nnt+1)),rms ENDDO ! K do i=1,2*nres do j=1,3 c(j,i)=c(j,i)/qpart csq(j,i)=csq(j,i)/qpart-c(j,i)**2 enddo c write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3) enddo do i=nnt,nct tempfac(1,i)=0.0d0 tempfac(2,i)=0.0d0 do j=1,3 tempfac(1,i)=tempfac(1,i)+csq(j,i) tempfac(2,i)=tempfac(2,i)+csq(j,i+nres) enddo tempfac(1,i)=dsqrt(tempfac(1,i)) tempfac(2,i)=dsqrt(tempfac(2,i)) enddo return end c------------------------------------------------------------------------------ subroutine fittest_coord(igr) implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.VAR' include 'COMMON.FFIELD' include 'COMMON.TORCNSTR' include 'COMMON.SAXS' logical non_conv double precision przes(3),obrot(3,3) double precision xx(3,maxres2),yy(3,maxres2) integer i,ii,j,k,icon,jcon,jconmin,igr double precision rms,rmsmin,cwork(3,maxres2) double precision ehpb,Esaxs_constr,edihcnstr rmsmin=1.0d10 jconmin=nconf(igr,1) DO K=1,LICZ(IGR) jcon=nconf(igr,k) do i=1,2*nres do j=1,3 c(j,i)=allcart(j,i,jcon) enddo enddo call int_from_cart1(.false.) esaxs_constr=0 ehpb=0 edihcnstr=0 if (nsaxs.gt.0) call e_saxs(Esaxs_constr) call edis(ehpb) if (ndih_constr.gt.0) call etor_constr(edihcnstr) rms=wsaxs*esaxs_constr+wstrain*ehpb+edihcnstr c write (iout,*) "Esaxs_constr",esaxs_constr," Ehpb",ehpb, c & " Edihcnstr",edihcnstr if (rms.lt.rmsmin) then jconmin=nconf(igr,k) rmsmin=rms endif ENDDO ! K write (iout,*) "fittest conformation",jconmin," penalty",rmsmin do i=1,2*nres do j=1,3 c(j,i)=allcart(j,i,jconmin) enddo enddo return end c------------------------------------------------------------------------------ subroutine closest_coord(igr) implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.VAR' logical non_conv double precision przes(3),obrot(3,3) integer i,ii,j,k,icon,jcon,jconmin,igr,ipermmin double precision rms,rmsmin,cwork(3,maxres2) double precision xx(3,maxres2),yy(3,maxres2) double precision rmscalc rmsmin=1.0d10 jconmin=nconf(igr,1) DO K=1,LICZ(IGR) jcon=nconf(igr,k) do i=1,2*nres do j=1,3 xx(j,i)=c(j,i) yy(j,i)=allcart(j,i,jcon) enddo enddo rms=rmscalc(xx(1,1),yy(1,1),przes,obrot,ipermmin) c write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin if (non_conv) print *,non_conv,icon,jcon if (rms.lt.rmsmin) then rmsmin=rms jconmin=jcon endif ENDDO ! K c write (iout,*) "rmsmin",rmsmin," rms",rms c call flush(iout) do i=1,2*nres do j=1,3 c(j,i)=allcart(j,i,jconmin) enddo enddo return end c------------------------------------------------------------------------------ subroutine center implicit none include 'DIMENSIONS' include 'sizesclu.dat' include 'COMMON.IOUNITS' include 'COMMON.CONTROL' include 'COMMON.CLUSTER' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.VAR' double precision przes(3) integer i,ii,j,k,icon,jcon,jconmin,igr przes=0.0d0 do j=1,3 do i=1,nres przes(j)=przes(j)+c(j,i) enddo enddo do j=1,3 przes(j)=przes(j)/nres enddo do i=1,2*nres do j=1,3 c(j,i)=c(j,i)-przes(j) enddo enddo return end