From 550f832d956213047c82b1c17397e45ded1534bf Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sat, 14 Mar 2020 19:07:59 +0100 Subject: [PATCH] cluster unres --- source/cluster/unres/src/CMakeLists.txt | 56 +++- source/cluster/unres/src/COMMON.CONTROL | 12 +- source/cluster/unres/src/COMMON.IOUNITS | 9 +- source/cluster/unres/src/DIMENSIONS | 2 + source/cluster/unres/src/arcos.f | 2 +- source/cluster/unres/src/cxread.F | 237 ++++++++++++++ source/cluster/unres/src/fitsq.f | 2 +- source/cluster/unres/src/hc.f | 2 + source/cluster/unres/src/intcor.f | 2 +- source/cluster/unres/src/main_clust.F | 525 +++++++++++++++++++++++++++++++ source/cluster/unres/src/main_clust.f | 498 ----------------------------- source/cluster/unres/src/readrtns.F | 81 ++++- source/cluster/unres/src/sizesclu.dat | 2 +- source/cluster/unres/src/timing.F | 1 + source/cluster/unres/src/wrtclust.f | 35 ++- 15 files changed, 934 insertions(+), 532 deletions(-) create mode 100644 source/cluster/unres/src/cxread.F create mode 100644 source/cluster/unres/src/main_clust.F delete mode 100644 source/cluster/unres/src/main_clust.f diff --git a/source/cluster/unres/src/CMakeLists.txt b/source/cluster/unres/src/CMakeLists.txt index 26a6b4d..c26fe29 100644 --- a/source/cluster/unres/src/CMakeLists.txt +++ b/source/cluster/unres/src/CMakeLists.txt @@ -18,7 +18,7 @@ set(UNRES_CLUSTER_UNRES_SRC0 hc.f initialize.f intcor.f - main_clust.f + main_clust.F matmult.f misc.f noyes.f @@ -30,6 +30,7 @@ set(UNRES_CLUSTER_UNRES_SRC0 timing.F track.F wrtclust.f + cxread.F ) set(UNRES_CLUSTER_PP_SRC @@ -56,8 +57,9 @@ set_property(SOURCE ${UNRES_CLUSTER_SRC0} PROPERTY COMPILE_FLAGS ${FFLAGS0} ) #========================================= # 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") @@ -71,7 +73,7 @@ set_property(SOURCE ${UNRES_CLUSTER_PP_SRC} PROPERTY COMPILE_DEFINITIONS ${CPPFL #======================================== # 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 @@ -85,9 +87,57 @@ add_executable(UNRES_CLUSTER_BIN ${UNRES_CLUSTER_SRCS} ) 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 ) diff --git a/source/cluster/unres/src/COMMON.CONTROL b/source/cluster/unres/src/COMMON.CONTROL index 1fa1444..a9bf4e7 100644 --- a/source/cluster/unres/src/COMMON.CONTROL +++ b/source/cluster/unres/src/COMMON.CONTROL @@ -1,7 +1,7 @@ - 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 diff --git a/source/cluster/unres/src/COMMON.IOUNITS b/source/cluster/unres/src/COMMON.IOUNITS index b80aa37..334014f 100644 --- a/source/cluster/unres/src/COMMON.IOUNITS +++ b/source/cluster/unres/src/COMMON.IOUNITS @@ -2,13 +2,14 @@ C----------------------------------------------------------------------- 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 diff --git a/source/cluster/unres/src/DIMENSIONS b/source/cluster/unres/src/DIMENSIONS index 44a47a7..37cdb6e 100644 --- a/source/cluster/unres/src/DIMENSIONS +++ b/source/cluster/unres/src/DIMENSIONS @@ -71,3 +71,5 @@ C Maximum number of chains parameter (mxch=1) C Maximum number of generated conformations parameter (mxio=1000) +C Maximum number of coordinate files + parameter (maxfiles=50) diff --git a/source/cluster/unres/src/arcos.f b/source/cluster/unres/src/arcos.f index 052a1e4..698f704 100644 --- a/source/cluster/unres/src/arcos.f +++ b/source/cluster/unres/src/arcos.f @@ -2,7 +2,7 @@ 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 diff --git a/source/cluster/unres/src/cxread.F b/source/cluster/unres/src/cxread.F new file mode 100644 index 0000000..5d4b511 --- /dev/null +++ b/source/cluster/unres/src/cxread.F @@ -0,0 +1,237 @@ + 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 diff --git a/source/cluster/unres/src/fitsq.f b/source/cluster/unres/src/fitsq.f index 17d92ee..fbcf223 100644 --- a/source/cluster/unres/src/fitsq.f +++ b/source/cluster/unres/src/fitsq.f @@ -1,5 +1,5 @@ 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 diff --git a/source/cluster/unres/src/hc.f b/source/cluster/unres/src/hc.f index 2e03f24..da10b8f 100644 --- a/source/cluster/unres/src/hc.f +++ b/source/cluster/unres/src/hc.f @@ -282,6 +282,7 @@ C statement was not necessary. C 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) @@ -410,6 +411,7 @@ C F. Murtagh, ESA/ESO/STECF, Garching, Feb. 1986.C C C C-------------------------------------------------C SUBROUTINE HCDEN(LEV,IORDER,HEIGHT,CRITVAL) + include 'DIMENSIONS' include 'COMMON.IOUNITS' CHARACTER*80 LINE INTEGER IORDER(LEV),HEIGHT(LEV) diff --git a/source/cluster/unres/src/intcor.f b/source/cluster/unres/src/intcor.f index 7d0b849..045cce1 100644 --- a/source/cluster/unres/src/intcor.f +++ b/source/cluster/unres/src/intcor.f @@ -8,7 +8,7 @@ c 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 diff --git a/source/cluster/unres/src/main_clust.F b/source/cluster/unres/src/main_clust.F new file mode 100644 index 0000000..3838eaf --- /dev/null +++ b/source/cluster/unres/src/main_clust.F @@ -0,0 +1,525 @@ +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 diff --git a/source/cluster/unres/src/main_clust.f b/source/cluster/unres/src/main_clust.f deleted file mode 100644 index e3dd395..0000000 --- a/source/cluster/unres/src/main_clust.f +++ /dev/null @@ -1,498 +0,0 @@ -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 diff --git a/source/cluster/unres/src/readrtns.F b/source/cluster/unres/src/readrtns.F index 9d55005..9838839 100644 --- a/source/cluster/unres/src/readrtns.F +++ b/source/cluster/unres/src/readrtns.F @@ -48,10 +48,22 @@ C Set up the time limit (caution! The time must be input in minutes!) 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 @@ -353,6 +365,37 @@ c---------------------------------------------------------------------------- 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' @@ -385,26 +428,44 @@ c---------------------------------------------------------------------------- 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 @@ -412,7 +473,7 @@ C Get the names and open the input files 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') diff --git a/source/cluster/unres/src/sizesclu.dat b/source/cluster/unres/src/sizesclu.dat index c58a5f2..cb6c733 100644 --- a/source/cluster/unres/src/sizesclu.dat +++ b/source/cluster/unres/src/sizesclu.dat @@ -4,7 +4,7 @@ * * Max. number of conformations in the data set. * - PARAMETER (MAXCONF=2500) + PARAMETER (MAXCONF=5000) * * Max. number of "distances" between conformations. * diff --git a/source/cluster/unres/src/timing.F b/source/cluster/unres/src/timing.F index ecab9b4..9bca353 100644 --- a/source/cluster/unres/src/timing.F +++ b/source/cluster/unres/src/timing.F @@ -170,6 +170,7 @@ C Next definitions for Linux C========================================================================= C subroutine dajczas(rntime,hrtime,mintime,sectime) + include 'DIMENSIONS' include 'COMMON.IOUNITS' real*8 rntime,hrtime,mintime,sectime hrtime=rntime/3600.0D0 diff --git a/source/cluster/unres/src/wrtclust.f b/source/cluster/unres/src/wrtclust.f index ca43876..3f057b2 100644 --- a/source/cluster/unres/src/wrtclust.f +++ b/source/cluster/unres/src/wrtclust.f @@ -25,20 +25,20 @@ c ICANT(I,J)=((NCON+NCON-J)*(J-1))/2+I-J 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 @@ -94,9 +94,17 @@ cd print '(3i4,f12.4)',ind,ii,jj,curr_dist 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) @@ -237,6 +245,16 @@ c Write conformations of the family i to PDB files 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) @@ -244,6 +262,9 @@ c Produce only a single PDB file for the leading member of the family 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' -- 1.7.9.5