Added main_clust.F
authorAdam Liwo <adam@piasek4.chem.univ.gda.pl>
Wed, 30 Sep 2015 10:02:39 +0000 (12:02 +0200)
committerAdam Liwo <adam@piasek4.chem.univ.gda.pl>
Wed, 30 Sep 2015 10:02:39 +0000 (12:02 +0200)
source/cluster/unres/src/cxread.F [new file with mode: 0644]
source/cluster/unres/src/main_clust.F [new file with mode: 0644]

diff --git a/source/cluster/unres/src/cxread.F b/source/cluster/unres/src/cxread.F
new file mode 100644 (file)
index 0000000..5d4b511
--- /dev/null
@@ -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/main_clust.F b/source/cluster/unres/src/main_clust.F
new file mode 100644 (file)
index 0000000..3838eaf
--- /dev/null
@@ -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