+++ /dev/null
-C
-C Program to cluster united-residue MCM results.
-C
- include 'DIMENSIONS'
- include 'sizesclu.dat'
- include 'COMMON.TIME1'
- include 'COMMON.INTERACT'
- include 'COMMON.NAMES'
- include 'COMMON.GEO'
- include 'COMMON.HEADER'
- include 'COMMON.CONTROL'
- include 'COMMON.CONTACTS'
- include 'COMMON.CHAIN'
- include 'COMMON.VAR'
- include 'COMMON.CLUSTER'
- include 'COMMON.IOUNITS'
- logical printang(max_cut)
- integer printpdb(max_cut)
- integer printmol2(max_cut)
- character*240 lineh
- REAL CRIT(maxconf),MEMBR(maxconf)
- REAL CRITVAL(maxconf-1)
- REAL DISS(maxdist)
- INTEGER IA(maxconf),IB(maxconf)
- INTEGER ICLASS(maxconf,maxconf-1),HVALS(maxconf-1)
- INTEGER IORDER(maxconf-1),HEIGHT(maxconf-1)
- DIMENSION NN(maxconf),DISNN(maxconf)
- LOGICAL FLAG(maxconf)
-
- double precision varia(maxvar)
- double precision hrtime,mintime,sectime
- logical eof
- double precision eee(n_ene,maxconf)
-
- call initialize
- call openunits
- call read_control
- call molread
- do i=1,nres
- phi(i)=0.0D0
- theta(i)=0.0D0
- alph(i)=0.0D0
- omeg(i)=0.0D0
- enddo
-
- print *,'MAIN: nnt=',nnt,' nct=',nct
-
- DO I=1,NCUT
- PRINTANG(I)=.FALSE.
- PRINTPDB(I)=0
- printmol2(i)=0
- IF (RCUTOFF(I).LT.0.0) THEN
- RCUTOFF(I)=ABS(RCUTOFF(I))
- PRINTANG(I)=.TRUE.
- PRINTPDB(I)=outpdb
- printmol2(i)=outmol2
- ENDIF
- ENDDO
- PRINT *,'Number of cutoffs:',NCUT
- PRINT *,'Cutoff values:'
- DO ICUT=1,NCUT
- PRINT '(8HRCUTOFF(,I2,2H)=,F8.1,2i2)',ICUT,RCUTOFF(ICUT),
- & printpdb(icut),printmol2(icut)
- ENDDO
- DO I=1,NRES-3
- MULT(I)=1
- ENDDO
- ICON=1
- 123 continue
- if (from_cart) then
- if (efree) then
- read (intin,*,end=13,err=11) energy(icon),totfree(icon),
- & rmstab(icon),
- & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
- & i=1,nss_all(icon)),iscore(icon)
- else
- read (intin,*,end=13,err=11) energy(icon),rmstab(icon),
- & nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),
- & i=1,nss_all(icon)),iscore(icon)
- endif
- read (intin,'(8f10.5)',end=13,err=10)
- & ((allcart(j,i,icon),j=1,3),i=1,nres),
- & ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
- print *,icon,energy(icon),nss_all(icon),rmstab(icon)
- else
- read(intin,'(a80)',end=13,err=12) lineh
- read(lineh(:5),*,err=8) ic
- if (efree) then
- read(lineh(6:),*,err=8) energy(icon)
- else
- read(lineh(6:),*,err=8) energy(icon)
- endif
- goto 9
- 8 ic=1
- print *,'error, assuming e=1d10',lineh
- energy(icon)=1d10
- nss=0
- 9 continue
-cold read(lineh(18:),*,end=13,err=11) nss_all(icon)
- ii = index(lineh(15:)," ")+15
- read(lineh(ii:),*,end=13,err=11) nss_all(icon)
- IF (NSS_all(icon).LT.9) THEN
- read (lineh(20:),*,end=102)
- & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),
- & iscore(icon)
- ELSE
- read (lineh(20:),*,end=102)
- & (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
- read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),
- & I=9,NSS_all(icon)),iscore(icon)
- ENDIF
-
- 102 continue
-
- PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
- totfree(icon)=energy(icon)
- call read_angles(intin,*13)
- call chainbuild
- do ii=1,2*nres
- do jj=1,3
- allcart(jj,ii,icon)=c(jj,ii)
- enddo
- enddo
- do i=1,nres
- phiall(i,icon)=phi(i)
- thetall(i,icon)=theta(i)
- alphall(i,icon)=alph(i)
- omall(i,icon)=omeg(i)
- enddo
- endif
- ICON=ICON+1
- GOTO 123
-C
-C CALCULATE DISTANCES
-C
- 10 print *,'something wrong with angles'
- goto 13
- 11 print *,'something wrong with NSS',nss
- goto 13
- 12 print *,'something wrong with header'
-
- 13 NCON=ICON-1
- if (lgrp) then
- i=0
- do while (i.lt.ncon)
-c read (jstatin,*,err=1111,end=1111) idum,(eee(j,i),j=1,n_ene)
- read (jstatin,'(a)') lineh
- if (index(lineh,'#').gt.0) goto 1123
- i=i+1
- read (lineh,*) idum,(eee(j,i),j=1,n_ene)
- print *,i,idum,(eee(j,i),j=1,n_ene)
- energy(i)=eee(15,i)
- 1123 continue
- enddo
- endif
- goto 1112
- 1111 print *,'Error in the statin file; statout will not be produced.'
- write (iout,*)
- & 'Error in the statin file; statout will not be produced.'
- lgrp=.false.
- 1112 continue
- print *,'ncon',ncon
- DO I=1,NCON
- ICC(I)=I
- ENDDO
- WRITE (iout,'(A80)') TITEL
- t1=tcpu()
- DO I=1,NCON-1
- if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
- if (from_cart) then
- do ii=1,2*nres
- do jj=1,3
- c(jj,ii)=allcart(jj,ii,i)
- enddo
- enddo
- call int_from_cart(.true.,.false.)
- do ii=1,nres
- phiall(ii,i)=phi(ii)
- thetall(ii,i)=theta(ii)
- alphall(ii,i)=alph(ii)
- omall(ii,i)=omeg(ii)
- enddo
- else
- do ii=1,nres
- phi(ii)=phiall(ii,i)
- theta(ii)=thetall(ii,i)
- alph(ii)=alphall(ii,i)
- omeg(ii)=omall(ii,i)
- enddo
- call chainbuild
- endif
- do ii=1,nres
- do jj=1,3
- cref(jj,ii)=c(jj,ii)
- enddo
- enddo
- DO J=I+1,NCON
- IND=IOFFSET(NCON,I,J)
- ATTALUMS(IND)=DIFCONF(I,J)
-c write (iout,'(2i4,f10.5)') i,j,ATTALUMS(IND)
- DISS(IND)=ATTALUMS(IND)
- ENDDO
- ENDDO
- t2=tcpu()
- WRITE (iout,'(/a,1pe14.5,a/)')
- & 'Time for distance calculation:',T2-T1,' sec.'
- t1=tcpu()
- PRINT '(a)','End of distance computation'
-
- if (punch_dist) then
- do i=1,ncon-1
- do j=i+1,ncon
- IND=IOFFSET(NCON,I,J)
- write (jrms,'(2i5,2f10.5)') i,j,attalums(IND),
- & energy(j)-energy(i)
- enddo
- enddo
- endif
-C
-C Print out the RMS deviation matrix.
-C
- if (print_dist) CALL DISTOUT(NCON)
-C
-C call hierarchical clustering HC from F. Murtagh
-C
- N=NCON
- LEN = (N*(N-1))/2
- write(iout,*) "-------------------------------------------"
- write(iout,*) "HIERARCHICAL CLUSTERING using"
- if (iopt.eq.1) then
- write(iout,*) "WARD'S MINIMUM VARIANCE METHOD"
- elseif (iopt.eq.2) then
- write(iout,*) "SINGLE LINK METHOD"
- elseif (iopt.eq.3) then
- write(iout,*) "COMPLETE LINK METHOD"
- elseif (iopt.eq.4) then
- write(iout,*) "AVERAGE LINK (OR GROUP AVERAGE) METHOD"
- elseif (iopt.eq.5) then
- write(iout,*) "MCQUITTY'S METHOD"
- elseif (iopt.eq.6) then
- write(iout,*) "MEDIAN (GOWER'S) METHOD"
- elseif (iopt.eq.7) then
- write(iout,*) "CENTROID METHOD"
- else
- write(iout,*) "IOPT=",iopt," IS INVALID, use 1-7"
- write(*,*) "IOPT=",iopt," IS INVALID, use 1-7"
- stop
- endif
- write(iout,*)
- write(iout,*) "hc.f by F. Murtagh, ESA/ESO/STECF, Garching"
- write(iout,*) "February 1986"
- write(iout,*) "References:"
- write(iout,*) "1. Multidimensional clustering algorithms"
- write(iout,*) " Fionn Murtagh"
- write(iout,*) " Vienna : Physica-Verlag, 1985."
- write(iout,*) "2. Multivariate data analysis"
- write(iout,*) " Fionn Murtagh and Andre Heck"
- write(iout,*) " Kluwer Academic Publishers, 1987"
- write(iout,*) "-------------------------------------------"
- write(iout,*)
-
- CALL HC(N,M,LEN,IOPT,IA,IB,CRIT,MEMBR,NN,DISNN,FLAG,DISS)
- LEV = N-1
- CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
- CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
-
- icut=1
- i=1
- NGR=i+1
- do j=1,n
- licz(iclass(j,i))=licz(iclass(j,i))+1
- nconf(iclass(j,i),licz(iclass(j,i)))=j
- enddo
- do i=1,lev-1
-
- idum=lev-i
- DO L=1,LEV
- IF (HEIGHT(L).EQ.IDUM) GOTO 190
- ENDDO
- 190 IDUM=L
-cd print *,i+1,CRITVAL(IDUM)
- IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN
- WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut)
- write (iout,'(a,f8.2)') 'Maximum distance found:',
- & CRITVAL(IDUM)
- CALL SRTCLUST(ICUT,ncon)
- CALL TRACK(ICUT)
- CALL WRTCLUST(ncon,ICUT,PRINTANG,PRINTPDB,PRINTMOL2)
- icut=icut+1
- if (icut.gt.ncut) goto 191
- ENDIF
- NGR=i+1
- do l=1,maxgr
- licz(l)=0
- enddo
- do j=1,n
- licz(iclass(j,i))=licz(iclass(j,i))+1
- nconf(iclass(j,i),licz(iclass(j,i)))=j
-cd print *,j,iclass(j,i),
-cd & licz(iclass(j,i)),nconf(iclass(j,i),licz(iclass(j,i)))
- enddo
- enddo
- 191 continue
-C
- if (plot_tree) then
- CALL WRITRACK
- CALL PLOTREE
- endif
-C
- t2=tcpu()
- WRITE (iout,'(/a,1pe14.5,a/)')
- & 'Total time for clustering:',T2-T1,' sec.'
- if (lgrp) then
- write (jstatout,'(a5,18a12,10f4.1)')"# ",
- & "EVDW SC-SC","EVDW2 SC-p","EES p-p",
- & "EBE bending","ESC SCloc","ETORS ",
- & "ECORR4 ","ECORR5 ","ECORR6 ",
- & "EELLO ","ETURN3 ","ETURN4 ","ETURN6 "
- & ,"ETOT total","RMSD","nat.contact","nnt.contact",
- & "cont.order",(rcutoff(i),i=1,ncut)
- do i=1,ncon
- write(jstatout,'(i5,18f12.4,10i6)') i,(eee(j,i),j=1,n_ene),
- & (iass_tot(i,icut),icut=1,ncut)
- enddo
- endif
-C
- stop '********** Program terminated normally.'
- end
-c---------------------------------------------------------------------------
- double precision function difconf(icon,jcon)
- include 'DIMENSIONS'
- include 'sizesclu.dat'
- 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)
- double precision xx(3,maxres2),yy(3,maxres2)
- do i=1,nres
- phi(i)=phiall(i,jcon)
- theta(i)=thetall(i,jcon)
- alph(i)=alphall(i,jcon)
- omeg(i)=omall(i,jcon)
- enddo
- call chainbuild
-c do i=1,nres
-c print '(i4,2(3f10.5,5x))',i,(cref(j,i),j=1,3),(c(j,i),j=1,3)
-c enddo
- if (lside) then
- ii=0
- do i=nnt,nct
- ii=ii+1
- do j=1,3
- xx(j,ii)=allcart(j,i,jcon)
- yy(j,ii)=cref(j,i)
- enddo
- enddo
- do i=nnt,nct
-c if (itype(i).ne.10) then
- ii=ii+1
- do j=1,3
- xx(j,ii)=allcart(j,i+nres,jcon)
- yy(j,ii)=cref(j,i+nres)
- enddo
-c endif
- enddo
- call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
- else
- do i=nnt,nct
- do j=1,3
- c(j,i)=allcart(j,i,jcon)
- enddo
-c write (2,'(i5,3f10.5,5x,3f10.5)') i,
-c & (c(j,i),j=1,3),(cref(j,i),j=1,3)
- enddo
- call fitsq(rms,c(1,nnt),cref(1,nnt),nct-nnt+1,przes,obrot
- & ,non_conv)
- endif
- if (rms.lt.0.0) then
- print *,error,'rms^2 = ',rms,icon,jcon
- stop
- endif
- if (non_conv) print *,non_conv,icon,jcon
- difconf=dsqrt(rms)
- RETURN
- END
-C------------------------------------------------------------------------------
- double precision function rmsnat(jcon)
- include 'DIMENSIONS'
- include 'sizesclu.dat'
- 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)
- double precision xx(3,maxres2),yy(3,maxres2)
- if (lside) then
- ii=0
- do i=nnt,nct
- ii=ii+1
- do j=1,3
- xx(j,ii)=allcart(j,i,jcon)
- yy(j,ii)=cref_pdb(j,i)
- enddo
- enddo
- do i=nnt,nct
-c if (itype(i).ne.10) then
- ii=ii+1
- do j=1,3
- xx(j,ii)=allcart(j,i+nres,jcon)
- yy(j,ii)=cref_pdb(j,i+nres)
- enddo
-c endif
- enddo
- call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
- else
- do i=nnt,nct
- do j=1,3
- c(j,i)=allcart(j,i,jcon)
- enddo
- enddo
- call fitsq(rms,c(1,nnt),cref_pdb(1,nnt),nct-nnt+1,przes,obrot
- & ,non_conv)
- endif
- if (rms.lt.0.0) then
- print *,error,'rms^2 = ',rms,icon,jcon
- stop
- endif
- if (non_conv) print *,non_conv,icon,jcon
- rmsnat=dsqrt(rms)
- RETURN
- END
-C------------------------------------------------------------------------------
- subroutine distout(ncon)
- include 'DIMENSIONS'
- include 'sizesclu.dat'
- parameter (ncol=10)
- include 'COMMON.IOUNITS'
- include 'COMMON.CLUSTER'
- dimension b(ncol)
- write (iout,'(a)') 'The distance matrix'
- do 1 i=1,ncon,ncol
- nlim=min0(i+ncol-1,ncon)
- write (iout,1000) (k,k=i,nlim)
- write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
- 1000 format (/8x,10(i4,3x))
- 1020 format (/1x,80(1h-)/)
- do 2 j=i,ncon
- jlim=min0(j,nlim)
- if (jlim.eq.j) then
- b(jlim-i+1)=0.0d0
- jlim1=jlim-1
- else
- jlim1=jlim
- endif
- do 3 k=i,jlim1
- if (j.lt.k) then
- IND=IOFFSET(NCON,j,k)
- else
- IND=IOFFSET(NCON,k,j)
- endif
- 3 b(k-i+1)=attalums(IND)
- write (iout,1010) j,(b(k),k=1,jlim-i+1)
- 2 continue
- 1 continue
- 1010 format (i5,3x,10(f6.2,1x))
- return
- end
-C-----------------------------------------------------------------------
- block data
- include 'DIMENSIONS'
- include 'COMMON.LOCAL'
- data dsc / 1.237, ! CYS (type 1)
- & 2.142, ! MET (type 2)
- & 2.299, ! PHE (type 3)
- & 1.776, ! ILE (type 4)
- & 1.939, ! LEU (type 5)
- & 1.410, ! VAL (type 6)
- & 2.605, ! TRP (type 7)
- & 2.484, ! TYR (type 8)
- & 0.743, ! ALA (type 9)
- & 0.000, ! GLY (type 10)
- & 1.393, ! THR (type 11)
- & 1.150, ! SER (type 12)
- & 2.240, ! GLN (type 13)
- & 1.684, ! ASN (type 14)
- & 2.254, ! GLU (type 15)
- & 1.709, ! ASP (type 16)
- & 2.113, ! HIS (type 17)
- & 3.020, ! ARG (type 18)
- & 2.541, ! LYS (type 19)
- & 1.345, ! PRO (type 20)
- & 0.000 /! D (type 21)
- end
lgrp=(index(controlcard,'LGRP').gt.0)
caonly=(index(controlcard,'CA_ONLY').gt.0)
print_dist=(index(controlcard,'PRINT_DIST').gt.0)
+ print *,"print_dist",print_dist
call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
call readi(controlcard,'IOPT',iopt,2)
+ print *,"iopt",iopt
lside = index(controlcard,"SIDE").gt.0
efree = index(controlcard,"EFREE").gt.0
+ print *,controlcard
+ call readi(controlcard,'ISTART',is,1)
+ print *,"is",is
+ call readi(controlcard,'IEND',ie,10000000)
+ print *,"ie",ie
+ call readi(controlcard,'ISAMPL',isampl,1)
+ print *,"isampl",isampl
+ call reada(controlcard,'TS',ts,1.0d10)
+ call reada(controlcard,'TE',te,1.0d10)
+ print *,"is",is," ie",ie
if (min_var) iopt=1
return
end
return
end
c----------------------------------------------------------------------------
+ subroutine split_string(rekord,tablica,dim,nsub)
+ implicit none
+ integer dim,nsub,i,ii,ll,kk
+ character*(*) tablica(dim)
+ character*(*) rekord
+ integer ilen
+ external ilen
+ do i=1,dim
+ tablica(i)=" "
+ enddo
+ ii=1
+ ll = ilen(rekord)
+ nsub=0
+ do i=1,dim
+C Find the start of term name
+ kk = 0
+ do while (ii.le.ll .and. rekord(ii:ii).eq." ")
+ ii = ii+1
+ enddo
+C Parse the name into TABLICA(i) until blank found
+ do while (ii.le.ll .and. rekord(ii:ii).ne." ")
+ kk = kk+1
+ tablica(i)(kk:kk)=rekord(ii:ii)
+ ii = ii+1
+ enddo
+ if (kk.gt.0) nsub=nsub+1
+ if (ii.gt.ll) return
+ enddo
+ return
+ end
+c----------------------------------------------------------------------------
subroutine card_concat(card)
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
call getenv('PDB',cfrom_pdb)
call getenv('PRINTCOOR',cprint)
from_cart = index(ucase(cfrom_pdb),'CART').gt.0
- lprint_cart = index(ucase(cprint),'PRINT_CART').gt.0
- lprint_int = index(ucase(cprint),'PRINT_INT').gt.0
- if (from_cart .and. .not.lprint_int) then
- lprint_cart=.true.
- lprint_int=.false.
- endif
- if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.true.
+ from_cx = index(ucase(cfrom_pdb),'CX').gt.0
+ lprint_cart = index(ucase(cprint),'CART').gt.0
+ lprint_int = index(ucase(cprint),'INT').gt.0
+c if (.not.lprint_cart .and. .not.lprint_int) lprint_int=.true.
lenpre=ilen(prefix)
lenout=ilen(prefout)
lenint=ilen(prefintin)
C Get the names and open the input files
open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
outname=prefout(:lenout)//'_clust.out'
- if (lprint_cart) then
- intname=prefintin(:lenint)//'_clust'//'.x'
+ if (from_cart) then
intinname=prefintin(:lenint)//'.x'
+ else if (from_cx) then
+ write (iout,*) "cx files: ",prefintin(:ilen(prefintin))
+ call split_string(prefintin,cxfiles(1),maxfiles,nfiles)
+ write (iout,*) "nfiles",nfiles
+ write (iout,*) "Split cxfiles"
+ do i=1,nfiles
+ cxfiles(i)=cxfiles(i)(:ilen(cxfiles(i)))//'.cx'
+ write (iout,*) cxfiles(i)(:ilen(cxfiles(i)))
+ enddo
else
- intname=prefintin(:lenint)//'_clust'//'.int'
intinname=prefintin(:lenint)//'.int'
endif
+ if (lprint_cart) then
+ if (.not. from_cx) then
+ intname=prefintin(:lenint)//'_clust'//'.x'
+ else
+ intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.x'
+ endif
+ else if (lprint_int) then
+ if (.not. from_cx) then
+ intname=prefintin(:lenint)//'_clust'//'.int'
+ else
+ intname=cxfiles(1)(:ilen(cxfiles(1)))//'_clust'//'.int'
+ endif
+ endif
rmsname=prefintin(:lenint)//'.rms'
statinname=prefintin(:lenint)//'.stat'
print *,statinname
open (jplot,file=prefout(:ilen(prefout))//'.tex',
& status='unknown')
print *,'unit',jplot,' opened'
- open (intin,file=intinname,status='old')
+ if (.not. from_cx) open (intin,file=intinname,status='old')
print *,'unit',intin,' opened'
open (jrms,file=rmsname,status='unknown')
open (jstatin,file=statinname,status='unknown')