hc.f
initialize.f
intcor.f
- main_clust.f
+ main_clust.F
matmult.f
misc.f
noyes.f
timing.F
track.F
wrtclust.f
+ cxread.F
)
set(UNRES_CLUSTER_PP_SRC
#=========================================
# System specific flags
#=========================================
+set(CPPFLAGS "NEWUNRES")
if(${CMAKE_SYSTEM_NAME} MATCHES "Linux")
- set(CPPFLAGS "${CPPFLAGS} LINUX")
+ set(CPPFLAGS "${CPPFLAGS} -DLINUX")
endif(${CMAKE_SYSTEM_NAME} MATCHES "Linux")
#========================================
# Setting binary name
#========================================
-set(UNRES_CLUSTER_BIN "cluster_unres_${Fortran_COMPILER_NAME}_${UNRES_MD_FF}.exe")
+set(UNRES_CLUSTER_BIN "cluster_unres_${Fortran_COMPILER_NAME}.exe")
#=========================================
# Set full unres CLUSTER sources
set_target_properties(UNRES_CLUSTER_BIN PROPERTIES OUTPUT_NAME ${UNRES_CLUSTER_BIN})
set_property(TARGET UNRES_CLUSTER_BIN PROPERTY RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin )
+#=========================================
+# Link libraries
+#=========================================
+# link libxdrf.a
+target_link_libraries( UNRES_CLUSTER_BIN xdrf )
+
#=========================================
# Install Path
#=========================================
install(TARGETS UNRES_CLUSTER_BIN DESTINATION ${CMAKE_INSTALL_PREFIX}/cluster)
+
+#=========================================
+# TESTS
+#=========================================
+
+FILE(WRITE ${CMAKE_CURRENT_BINARY_DIR}/scripts/cluster_int.sh
+"#!/bin/sh
+export INPUT=$1
+export INTIN=1l2y_csa_GB000
+export OUTPUT=1l2y_csa_GB000
+#-----------------------------------------------------------------------------
+CLUSTER_BIN=${CMAKE_BINARY_DIR}/bin/${UNRES_CLUSTER_BIN}
+#-----------------------------------------------------------------------------
+echo CTEST_FULL_OUTPUT
+$CLUSTER_BIN
+./cluster_check.sh $1
+")
+
+#
+# File permissions workaround
+#
+FILE( COPY ${CMAKE_CURRENT_BINARY_DIR}/scripts/cluster_int.sh
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
+ FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE
+)
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/cluster_check.sh
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
+ FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE
+)
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1l2y_clust_int.inp
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1l2y_csa_GB000.int
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
+FILE(COPY ${CMAKE_SOURCE_DIR}/ctest/1L2Y.pdb
+ DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+
+
+add_test(NAME CLUSTER_INT COMMAND sh ${CMAKE_CURRENT_BINARY_DIR}/cluster_int.sh 1l2y_clust_int )
- double precision betaT
- integer iscode,indpdb,outpdb,outmol2,iopt
+ double precision betaT,ts,te
+ integer iscode,indpdb,outpdb,outmol2,iopt,is,ie,isampl
logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
- & lprint_cart,lprint_int,from_cart,efree
- common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2,
- & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
- & from_cart,efree,iopt
+ & lprint_cart,lprint_int,from_cart,from_cx,efree
+ common /cntrl/ betaT,ts,te,is,ie,isampl,iscode,indpdb,refstr,
+ & pdbref,outpdb,outmol2,punch_dist,print_dist,caonly,lside,
+ & lprint_cart,lprint_int,from_cart,from_cx,efree,iopt
C I/O units used by the program
C-----------------------------------------------------------------------
integer inp,iout,igeom,intin,ipdb,imol2,ipdbin,jplot,jstatin,
- & jstatout
+ & jstatout,nfiles
common /iounits/ inp,iout,igeom,intin,ipdb,imol2,ipdbin,jplot,
- & jrms,jstatin,jstatout
+ & jrms,jstatin,jstatout,nfiles
character*256 outname,intname,intinname,pdbname,mol2name,prefix,
- & prefout,prefintin,statname,rmsname,statinname,statoutname
+ & prefout,prefintin,statname,rmsname,statinname,statoutname,
+ & cxfiles(maxfiles)
common /fnames/ outname,intname,pdbname,mol2name,intinname,prefix,
- & prefout,prefintin,statinname,statoutname
+ & prefout,prefintin,statinname,statoutname,cxfiles
C-----------------------------------------------------------------------
C INP - main input file
C IOUT - list file
parameter (mxch=1)
C Maximum number of generated conformations
parameter (mxio=1000)
+C Maximum number of coordinate files
+ parameter (maxfiles=50)
implicit real*8 (a-h,o-z)
include 'COMMON.GEO'
IF (DABS(X).LT.1.0D0) GOTO 1
- ARCOS=0.5D0*(PI+DSIGN(X,1.0D0)*PI)
+ ARCOS=0.5D0*(PI-DSIGN(X,1.0D0)*PI)
RETURN
1 ARCOS=DACOS(X)
RETURN
--- /dev/null
+ subroutine cxread(icon,*)
+ include 'DIMENSIONS'
+ include 'sizesclu.dat'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.NAMES'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.HEADER'
+ include 'COMMON.SBRIDGE'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.CLUSTER'
+ character*64 nazwa
+ real*4 rtime,rpotE,ruconst,rt_bath,rprop(20)
+ double precision time
+ integer iret,itmp
+ real xoord(3,maxres2+2),prec
+ double precision cm(3)
+ integer nstep
+ integer ilen
+ external ilen
+ integer icon
+
+c print *,"is",is," ie",ie," isampl",isampl
+ print *,nazwa
+ nstep=0
+ icon=0
+ nprop=0
+ nprop_prev=0
+ do i=1,20
+ rprop(i)=0.0d0
+ enddo
+
+ DO IFILE = 1, NFILES
+
+ print *,"CXREAD: opening file ",
+ & cxfiles(ifile)(:ilen(cxfiles(ifile)))
+ write (iout,*) "CXREAD: opening file ",
+ & cxfiles(ifile)(:ilen(cxfiles(ifile)))
+#if (defined(AIX) && !defined(JUBL))
+ call xdrfopen_(ixdrf,cxfiles(ifile), "r", iret)
+#else
+ call xdrfopen(ixdrf,cxfiles(ifile), "r", iret)
+#endif
+ if (iret.eq.0) cycle
+
+ print *,"CXREAD: reading file ",
+ & cxfiles(ifile)(:ilen(cxfiles(ifile)))
+ write(iout,*) "CXREAD: reading file ",
+ & cxfiles(ifile)(:ilen(cxfiles(ifile)))
+
+ do while (iret.gt.0)
+
+#if (defined(AIX) && !defined(JUBL))
+ call xdrffloat_(ixdrf, rtime, iret)
+ call xdrffloat_(ixdrf, rpotE, iret)
+#ifdef DEBUG
+ write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
+#endif
+ call flush(iout)
+ call xdrffloat_(ixdrf, ruconst, iret)
+#ifdef NEWUNRES
+ call xdrffloat(ixdrf, ruconst_back, iret)
+c print *,"uconst_back",ruconst_back
+#endif
+ call xdrffloat_(ixdrf, rt_bath, iret)
+ call xdrfint_(ixdrf, nss, iret)
+#ifdef DEBUG
+ write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
+#endif
+ do j=1,nss
+ call xdrfint_(ixdrf, ihpb(j), iret)
+ call xdrfint_(ixdrf, jhpb(j), iret)
+ enddo
+ call xdrfint_(ixdrf, nprop, iret)
+ do i=1,nprop
+ call xdrffloat_(ixdrf, rprop(i), iret)
+ enddo
+#else
+ call xdrffloat(ixdrf, rtime, iret)
+ call xdrffloat(ixdrf, rpotE, iret)
+#ifdef DEBUG
+ write (iout,*) "rtime",rtime," rpotE",rpotE," iret",iret
+#endif
+ call flush(iout)
+ call xdrffloat(ixdrf, ruconst, iret)
+#ifdef NEWUNRES
+ call xdrffloat(ixdrf, ruconst_back, iret)
+c print *,"uconst_back",ruconst_back
+#endif
+ call xdrffloat(ixdrf, rt_bath, iret)
+ call xdrfint(ixdrf, nss, iret)
+#ifdef DEBUG
+ write (iout,*) "ruconst",ruconst," rt_bath",rt_bath," nss",nss
+#endif
+ do j=1,nss
+ call xdrfint(ixdrf, ihpb(j), iret)
+ call xdrfint(ixdrf, jhpb(j), iret)
+ enddo
+ call xdrfint(ixdrf, nprop, iret)
+c write (iout,*) "nprop",nprop
+ if (it.gt.0 .and. nprop.ne.nprop_prev) then
+ write (iout,*) "Warning previous nprop",nprop_prev,
+ & " current",nprop
+ nprop=nprop_prev
+ else
+ nprop_prev=nprop
+ endif
+ do i=1,nprop
+ call xdrffloat(ixdrf, rprop(i), iret)
+ enddo
+#endif
+ if (iret.eq.0) exit
+#ifdef DEBUG
+ write (iout,*) rtime,rpotE,rt_bath,nss,
+ & (ihpb(j),jhpb(j),j=1,nss),(rprop(j),j=1,nprop)
+ write (iout,*) "nprop",nprop
+ call flush(iout)
+#endif
+ prec=10000.0
+
+ itmp=0
+#if (defined(AIX) && !defined(JUBL))
+ call xdrf3dfcoord_(ixdrf, xoord, itmp, prec, iret)
+#else
+ call xdrf3dfcoord(ixdrf, xoord, itmp, prec, iret)
+#endif
+#ifdef DEBUG
+ write (iout,'(10f8.3)') ((xoord(j,i),j=1,3),i=1,itmp)
+#endif
+ if (iret.eq.0) exit
+ if (itmp .ne. nres + nct - nnt + 1) then
+ write (iout,*) "Error: inconsistent sizes",itmp,nres+nct-nnt+1
+ call flush(iout)
+ exit
+ endif
+
+ time=rtime
+
+ do i=1,3
+ cm(i)=0.0d0
+ enddo
+
+ do i=1,nres
+ do j=1,3
+ c(j,i)=xoord(j,i)
+ cm(j)=cm(j)+c(j,i)
+ enddo
+ enddo
+ do i=1,nct-nnt+1
+ do j=1,3
+ c(j,i+nres+nnt-1)=xoord(j,i+nres)
+ enddo
+ enddo
+
+ do i=1,3
+ cm(i)=cm(i)/nres
+ enddo
+
+ do i=1,nres
+ do j=1,3
+ c(j,i)=c(j,i)-cm(j)
+ enddo
+ enddo
+ do i=1,nct-nnt+1
+ do j=1,3
+ c(j,i+nres+nnt-1)=c(j,i+nres+nnt-1)-cm(j)
+ enddo
+ enddo
+
+ nstep=nstep+1
+
+ if (nstep.gt.ie .or. rtime.gt.te) return
+
+ if((nstep.ge.is.or.rtime.ge.ts) .and. mod(nstep,isampl).eq.0)then
+
+ icon=icon+1
+#ifdef DEBUG
+ write (iout,*) "conformation, record",nstep,icon
+ write (iout,*) "pote",rpotE," time",rtime
+c write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss",nss
+c write (iout,*) "length",nres*4+(nct-nnt+1)*4+4+2*nss*4
+ call flush(iout)
+#endif
+ energy(icon)=rpotE
+ totfree(icon)=rpotE
+ rmstab(icon)=rmsdev
+ nss_all(icon)=nss
+ do k=1,nss
+ ihpb_all(k,icon)=ihpb(k)
+ jhpb_all(k,icon)=jhpb(k)
+ enddo
+ iscore(icon)=icon
+ do k=1,2*nres
+ do l=1,3
+ allcart(l,k,icon)=c(l,k)
+ enddo
+ enddo
+
+#ifdef DEBUG
+ call int_from_cart(.true.,.false.)
+ write (iout,*) "Storing conformation, record",icon
+ write (iout,*) "Cartesian coordinates"
+ write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
+ write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
+ write (iout,*) "Internal coordinates"
+ write (iout,'(8f10.4)') (dist(k-1,k),k=nnt+1,nct)
+ write (iout,'(8f10.4)') (dist(k,k+nres),k=nnt,nct)
+ write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
+ write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
+ write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
+ write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
+ write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
+c write (iout,'(8f10.5)') (rprop(j),j=1,nQ)
+ write (iout,'(16i5)') iscor
+ call flush(iout)
+#endif
+ endif
+
+ 112 continue
+
+ enddo
+
+#if (defined(AIX) && !defined(JUBL))
+ call xdrfclose_(ixdrf, iret)
+#else
+ call xdrfclose(ixdrf, iret)
+#endif
+ write (iout,*) nstep," conformations read so far file",
+ & cxfiles(ifile)(:ilen(cxfiles(ifile)))
+ call flush(iout)
+
+ ENDDO ! IFILE
+
+ return
+ end
subroutine fitsq(rms,x,y,nn,t,b,non_conv)
- implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
include 'COMMON.IOUNITS'
c x and y are the vectors of coordinates (dimensioned (3,n)) of the two
c structures to be superimposed. nn is 3*n, where n is the number of
C---------------------------------------------------------------C
SUBROUTINE HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,
X CRITVAL,HEIGHT)
+ include 'DIMENSIONS'
include 'sizesclu.dat'
include 'COMMON.IOUNITS'
integer ICLASS(maxconf,maxconf-1)
C C
C-------------------------------------------------C
SUBROUTINE HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
+ include 'DIMENSIONS'
include 'COMMON.IOUNITS'
CHARACTER*80 LINE
INTEGER IORDER(LEV),HEIGHT(LEV)
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.CHAIN'
- if (i1.lt.1 .or. i2.lt.1 .or. i3.lt.1 .or. i4.lt.1) then
+ if (i1.lt.1 .or. i2.lt.1 .or. i3.lt.1) then
alpha=0.0
return
endif
--- /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)
+ integer len
+
+ 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
+ if (from_cx) then
+
+ call cxread(ncon,*101)
+ goto 111
+ 101 write (iout,*) "Error in CX file"
+ stop
+ 111 continue
+
+ else
+
+ 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
+ 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
+
+ endif
+
+#ifdef DEBUG
+ do icon=1,ncon
+ write (iout,*) "Conformation",icon
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,
+ & (allcart(j,i,icon),j=1,3),(allcart(j,i+nres,icon),j=1,3)
+ enddo
+ enddo
+#endif
+ 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()
+C
+C CALCULATE DISTANCES
+C
+ DO I=1,NCON-1
+ if (mod(i,100).eq.0) print *,'Calculating RMS i=',i
+ if (from_cart .or. from_cx) 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)
+c print *,"HC"
+ LEV = N-1
+ CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT)
+c print *,"HCASS"
+c print *,"lev",lev
+c print *,"call HCDEN"
+ CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL)
+c print *,"HCDEN"
+
+ 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)
+c do i=1,nres
+c phi(i)=phiall(i,jcon)
+c theta(i)=thetall(i,jcon)
+c alph(i)=alphall(i,jcon)
+c omeg(i)=omall(i,jcon)
+c enddo
+c 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)=allcart(j,i,icon)
+ 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)=allcart(j,i+nres,icon)
+ enddo
+c endif
+ enddo
+ call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
+ else
+ do i=nnt,nct
+c do j=1,3
+c c(j,i)=allcart(j,i,jcon)
+c 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,allcart(1,nnt,jcon),allcart(1,nnt,icon),
+ & 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
+++ /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')
*
* Max. number of conformations in the data set.
*
- PARAMETER (MAXCONF=2500)
+ PARAMETER (MAXCONF=5000)
*
* Max. number of "distances" between conformations.
*
C=========================================================================
C
subroutine dajczas(rntime,hrtime,mintime,sectime)
+ include 'DIMENSIONS'
include 'COMMON.IOUNITS'
real*8 rntime,hrtime,mintime,sectime
hrtime=rntime/3600.0D0
C
C PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
C
- ii1= index(intinname,'/')
+ ii1= index(intname,'/')
ii2=ii1
ii1=ii1+1
do while (ii2.gt.0)
ii1=ii1+ii2
- ii2=index(intinname(ii1:),'/')
+ ii2=index(intname(ii1:),'/')
enddo
- ii = ii1+index(intinname(ii1:),'.')-1
+ ii = ii1+index(intname(ii1:),'.')-1
if (ii.eq.0) then
- ii=ilen(intinname)
+ ii=ilen(intname)
else
ii=ii-1
endif
- prefixp=intinname(ii1:ii)
+ prefixp=intname(ii1:ii)
cd print *,icut,printang(icut),printpdb(icut),printmol2(icut)
cd print *,'ecut=',ecut
WRITE (iout,100) NGR
write (iout,'(/A,F8.1,A,F8.1)')
& 'Max. distance in the family:',amax_dim,
& '; average distance in the family:',ave_dim
- if (refstr .or. pdbref) write (iout,'(a,i5,f8.3)')
+ if (refstr .or. pdbref) then
+ write (iout,'(a,i5,2f8.3)')
& "RMSD of the lowest-energy conformation #",nconf(igr,1),
- & rmsnat(nconf(igr,1))
+ & rmsnat(nconf(igr,1)),rmstab(nconf(igr,1))
+ rmsave=0.0d0
+ do i=1,licz(igr)
+ rmsave=rmsave+rmsnat(nconf(igr,i))
+ enddo
+ rmsave=rmsave/licz(igr)
+ write (iout,'(a,f8.3)') "Average RMSD in the family",rmsave
+ endif
19 CONTINUE
WRITE (iout,400)
WRITE (iout,500) (I,IASS(I),I=1,NCON)
else
c Produce only a single PDB file for the leading member of the family
write (iout,*) 'Writing pdb file: icon=',icon
+ if (from_cart .or. from_cx) then
+
+ do ii=1,2*nres
+ do j=1,3
+ c(j,ii)=allcart(j,ii,icon)
+ enddo
+ enddo
+
+ else
+
do ii=1,nres
phi(ii)=phiall(ii,icon)
theta(ii)=thetall(ii,icon)
omeg(ii)=omall(ii,icon)
enddo
call chainbuild
+
+ endif
+
cfname=prefixp(:ilen(prefixp))//numm(:ilen(numm))//exten
OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
c print *,'Calling pdbout'
projects / unres.git / commitdiff