rename
[unres4.git] / source / cluster / io_clust.f90
diff --git a/source/cluster/io_clust.f90 b/source/cluster/io_clust.f90
deleted file mode 100644 (file)
index db74bf6..0000000
+++ /dev/null
@@ -1,1824 +0,0 @@
-      module io_clust
-!-----------------------------------------------------------------------------
-      use clust_data
-      use io_units
-!      use names
-      use io_base !, only: ilen
-      use geometry_data, only: nres,c
-      use energy_data, only: nnt,nct,nss
-      use control_data, only: lside
-      implicit none
-!-----------------------------------------------------------------------------
-!
-!
-!-----------------------------------------------------------------------------
-      contains
-!-----------------------------------------------------------------------------
-! wrtclust.f
-!-----------------------------------------------------------------------------
-      SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
-
-      use hc_, only: ioffset
-      use control_data, only: lprint_cart,lprint_int,titel
-      use geometry, only: int_from_cart1,nres
-!      implicit real*8 (a-h,o-z)
-!      include 'DIMENSIONS'
-!      include 'sizesclu.dat'
-      integer,parameter :: num_in_line=5
-      LOGICAL :: PRINTANG(max_cut)
-      integer :: PRINTPDB(max_cut),printmol2(max_cut)
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.HEADER'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.VAR'
-!      include 'COMMON.CLUSTER'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.GEO'
-!      include 'COMMON.FREE'
-!      include 'COMMON.TEMPFAC'
-      real(kind=8) :: rmsave(maxgr)
-      CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
-      character(len=80) :: cfname
-      character(len=8) :: ctemper
-      DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
-           MUMM /'000'/
-!      external ilen
-      integer :: ncon,icut,ib
-      integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
-                 irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
-      real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
-                 ave_dim,amax_dim,emin1
-  
-
-      allocate(tempfac(2,nres))
-
-      do i=1,64
-        cfname(i:i)=" "
-      enddo
-!      print *,"calling WRTCLUST",ncon
-!      write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
-      rewind 80
-      call flush(iout)
-      temper=1.0d0/(beta_h(ib)*1.987d-3)
-      if (temper.lt.100.0d0) then
-        write(ctemper,'(f3.0)') temper
-        ctemper(3:3)=" "
-      else if (temper.lt.1000.0) then
-        write (ctemper,'(f4.0)') temper
-        ctemper(4:4)=" "
-      else
-        write (ctemper,'(f5.0)') temper
-        ctemper(5:5)=" "
-      endif
-
-      do i=1,ncon*(ncon-1)/2
-        read (80) diss(i)
-      enddo
-      close(80,status='delete')
-!
-!  PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
-!
-      ii1= index(intinname,'/')
-      ii2=ii1
-      ii1=ii1+1
-      do while (ii2.gt.0) 
-        ii1=ii1+ii2
-        ii2=index(intinname(ii1:),'/')
-      enddo 
-      ii = ii1+index(intinname(ii1:),'.')-1
-      if (ii.eq.0) then
-        ii=ilen(intinname)
-      else
-        ii=ii-1
-      endif
-      prefixp=intinname(ii1:ii)
-!d    print *,icut,printang(icut),printpdb(icut),printmol2(icut)
-!d    print *,'ecut=',ecut
-      WRITE (iout,100) NGR
-      DO 19 IGR=1,NGR
-      WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
-      NRECORD=LICZ(IGR)/num_in_line
-      IND1=1
-      DO 63 IRECORD=1,NRECORD
-      IND2=IND1+num_in_line-1
-      WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
-          totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
-      IND1=IND2+1
-   63 CONTINUE
-      WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
-         totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
-      IND1=1
-      ICON=list_conf(NCONF(IGR,1))
-!      WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
-! 12/8/93 Estimation of "diameters" of the subsequent families.
-      ave_dim=0.0
-      amax_dim=0.0
-!      write (iout,*) "ecut",ecut
-      do i=2,licz(igr)
-        ii=nconf(igr,i)
-        if (totfree(ii)-emin .gt. ecut) goto 10
-        do j=1,i-1
-          jj=nconf(igr,j)
-          if (jj.eq.1) exit
-          if (ii.lt.jj) then
-            ind=ioffset(ncon,ii,jj)
-          else
-            ind=ioffset(ncon,jj,ii)
-          endif
-!          write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
-!     &     " ind",ind
-          call flush(iout)
-          curr_dist=dabs(diss(ind)+0.0d0)
-!          write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
-!     &      list_conf(jj),curr_dist
-          if (curr_dist .gt. amax_dim) amax_dim=curr_dist
-          ave_dim=ave_dim+curr_dist**2
-        enddo
-      enddo   
-   10 if (licz(igr) .gt. 1) &
-       ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
-      write (iout,'(/A,F8.1,A,F8.1)') &
-       'Max. distance in the family:',amax_dim,&
-       '; average distance in the family:',ave_dim 
-      rmsave(igr)=0.0d0
-      qpart=0.0d0
-      do i=1,licz(igr)
-        icon=nconf(igr,i)
-        boltz=dexp(-totfree(icon))
-        rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
-        qpart=qpart+boltz
-      enddo
-      rmsave(igr)=rmsave(igr)/qpart
-      write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
-   19 CONTINUE
-      WRITE (iout,400)
-      WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
-!      print *,icut,printang(icut)
-      IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
-        emin=totfree_gr(1)
-!        print *,'emin',emin,' ngr',ngr
-        if (lprint_cart) then
-          cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-            //"K"//".x"
-        else
-          cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-            //"K"//".int"
-        endif
-        do igr=1,ngr
-          icon=nconf(igr,1)
-          if (totfree_gr(igr)-emin.le.ecut) then
-            if (lprint_cart) then
-              call cartout(igr,icon,totfree(icon)/beta_h(ib),&
-                totfree_gr(igr)/beta_h(ib),&
-                rmstb(icon),cfname)
-            else 
-!              print '(a)','calling briefout'
-              do i=1,2*nres
-                do j=1,3
-                  c(j,i)=allcart(j,i,icon)
-                enddo
-              enddo
-              call int_from_cart1(.false.)
-!el              call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
-!el                totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
-!el                jhpb_all(1,icon),cfname)
-              call briefout(igr,totfree(icon)/beta_h(ib),&
-                totfree_gr(igr))
-!              print '(a)','exit briefout'
-            endif
-          endif
-        enddo
-        close(igeom)
-      ENDIF
-      IF (PRINTPDB(ICUT).gt.0) THEN
-! Write out a number of conformations from each family in PDB format and
-! create InsightII command file for their displaying in different colors
-        cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-          //"K_"//'ave'//exten
-        write (iout,*) "cfname",cfname
-        OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
-        write (ipdb,'(a,f8.2)') &
-          "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
-        close (ipdb)
-        I=1
-        ICON=NCONF(1,1)
-        EMIN=totfree_gr(I)
-        emin1=totfree(icon)
-        DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
-!          write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
-          write (NUMM,'(bz,i4.4)') i
-          ncon_lim=min0(licz(i),printpdb(icut))
-          cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-            //"K_"//numm(:ilen(numm))//exten
-          OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
-          write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
-           i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
-! Write conformations of the family i to PDB files
-          ncon_out=1
-          do while (ncon_out.lt.printpdb(icut) .and. &
-           ncon_out.lt.licz(i).and. &
-           totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
-            ncon_out=ncon_out+1
-!            write (iout,*) i,ncon_out,nconf(i,ncon_out),
-!     &        totfree(nconf(i,ncon_out)),emin1,ecut
-          enddo
-          write (iout,*) "ncon_out",ncon_out
-          call flush(iout)
-          do j=1,nres
-            tempfac(1,j)=5.0d0
-            tempfac(2,j)=5.0d0
-          enddo
-          do j=1,ncon_out
-            icon=nconf(i,j)
-            do ii=1,2*nres
-              do k=1,3
-                c(k,ii)=allcart(k,ii,icon)
-              enddo
-            enddo
-            call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
-            write (ipdb,'("TER")')
-          enddo
-          close(ipdb)
-! Average structures and structures closest to average
-          cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-          //"K_"//'ave'//exten
-          OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
-           position="APPEND")
-          call ave_coord(i)
-          write (ipdb,'(a,i5)') "REMARK CLUSTER",i
-          call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
-          write (ipdb,'("TER")')
-          call closest_coord(i)
-          call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
-          write (ipdb,'("TER")')
-          close (ipdb)
-          I=I+1
-          ICON=NCONF(I,1)
-          emin1=totfree(icon)
-        ENDDO
-      ENDIF 
-      IF (printmol2(icut).gt.0) THEN
-! Write out a number of conformations from each family in PDB format and
-! create InsightII command file for their displaying in different colors
-        I=1
-        ICON=NCONF(1,1)
-        EMIN=ENERGY(ICON)
-        emin1=emin
-        DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
-          write (NUMM,'(bz,i4.4)') i
-          cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
-            //"K_"//numm(:ilen(numm))//extmol
-          OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
-          ncon_out=1
-          do while (ncon_out.lt.printmol2(icut) .and. &
-           ncon_out.lt.licz(i).and. &
-           totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
-            ncon_out=ncon_out+1
-          enddo
-          do j=1,ncon_out
-            icon=nconf(i,j)
-            do ii=1,2*nres
-              do k=1,3
-                c(k,ii)=allcart(k,ii,icon)
-              enddo
-            enddo
-            CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
-          enddo
-          CLOSE(imol2)
-          I=I+1
-          ICON=NCONF(I,1)
-          emin1=totfree(icon)
-        ENDDO
-      ENDIF 
-  100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
-  200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
-       ' CONTAINS ',I4,' CONFORMATION(S): ')
-! 300 FORMAT ( 8(I4,F6.1))
-  300 FORMAT (5(I4,1pe12.3))
-  400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
-  500 FORMAT (8(2I4,2X)) 
-  600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
-      RETURN
-      END SUBROUTINE WRTCLUST
-!------------------------------------------------------------------------------
-      subroutine ave_coord(igr)
-
-      use control_data, only:lside
-      use regularize_, only:fitsq,matvec
-!      implicit none
-!      include 'DIMENSIONS'
-!      include 'sizesclu.dat'
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.CLUSTER'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.INTERACT'
-!      include 'COMMON.VAR'
-!      include 'COMMON.TEMPFAC'
-!      include 'COMMON.IOUNITS'
-      logical :: non_conv
-      real(kind=8) :: przes(3),obrot(3,3)
-      real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
-      real(kind=8) :: eref
-      integer :: i,ii,j,k,icon,jcon,igr
-      real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
-!      write (iout,*) "AVE_COORD: igr",igr
-      jcon=nconf(igr,1)
-      eref=totfree(jcon)
-      boltz = dexp(-totfree(jcon)+eref)
-      qpart=boltz
-      do i=1,2*nres
-        do j=1,3
-          c(j,i)=allcart(j,i,jcon)*boltz
-          cref1(j,i)=allcart(j,i,jcon)
-          csq(j,i)=allcart(j,i,jcon)**2*boltz
-        enddo
-      enddo
-      DO K=2,LICZ(IGR)
-      jcon=nconf(igr,k)
-      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)=cref1(j,i)
-          enddo
-        enddo
-        do i=nnt,nct
-!          if (itype(i).ne.10) then
-            ii=ii+1
-            do j=1,3
-              xx(j,ii)=allcart(j,i+nres,jcon)
-              yy(j,ii)=cref1(j,i+nres)
-            enddo
-!          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
-            cwork(j,i)=allcart(j,i,jcon)
-          enddo
-        enddo
-        call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
-             ,non_conv)
-      endif
-!      write (iout,*) "rms",rms
-!      do i=1,3
-!        write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
-!      enddo
-      if (rms.lt.0.0) then
-        print *,'error, rms^2 = ',rms,icon,jcon
-        stop
-      endif
-      if (non_conv) print *,non_conv,icon,jcon
-      boltz=dexp(-totfree(jcon)+eref)
-      qpart = qpart + boltz
-      do i=1,2*nres
-        do j=1,3
-          xx(j,i)=allcart(j,i,jcon)
-        enddo
-      enddo
-      call matvec(cwork,obrot,xx,2*nres)
-      do i=1,2*nres
-!        write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
-!     &    (allcart(j,i,jcon),j=1,3)
-        do j=1,3
-          cwork(j,i)=cwork(j,i)+przes(j)
-          c(j,i)=c(j,i)+cwork(j,i)*boltz
-          csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz 
-        enddo
-      enddo
-      ENDDO ! K
-      do i=1,2*nres
-        do j=1,3
-          c(j,i)=c(j,i)/qpart
-          csq(j,i)=csq(j,i)/qpart-c(j,i)**2
-        enddo
-!        write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
-      enddo
-      do i=nnt,nct
-        tempfac(1,i)=0.0d0
-        tempfac(2,i)=0.0d0
-        do j=1,3
-          tempfac(1,i)=tempfac(1,i)+csq(j,i)
-          tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
-        enddo
-        tempfac(1,i)=dsqrt(tempfac(1,i))
-        tempfac(2,i)=dsqrt(tempfac(2,i))
-      enddo
-      return
-      end subroutine ave_coord
-!------------------------------------------------------------------------------
-      subroutine closest_coord(igr)
-
-      use regularize_, only:fitsq
-!      implicit none
-!      include 'DIMENSIONS'
-!      include 'sizesclu.dat'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.CLUSTER'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.INTERACT'
-!      include 'COMMON.VAR'
-      logical :: non_conv
-      real(kind=8) :: przes(3),obrot(3,3)
-      real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
-      integer :: i,ii,j,k,icon,jcon,jconmin,igr
-      real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
-      rmsmin=1.0d10
-      jconmin=nconf(igr,1)
-      DO K=1,LICZ(IGR)
-      jcon=nconf(igr,k)
-      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)=c(j,i)
-          enddo
-        enddo
-        do i=nnt,nct
-!          if (itype(i).ne.10) then
-            ii=ii+1
-            do j=1,3
-              xx(j,ii)=allcart(j,i+nres,jcon)
-              yy(j,ii)=c(j,i+nres)
-            enddo
-!          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
-            cwork(j,i)=allcart(j,i,jcon)
-          enddo
-        enddo
-        call fitsq(rms,cwork(1,nnt),c(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
-!      write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
-      if (non_conv) print *,non_conv,icon,jcon
-      if (rms.lt.rmsmin) then
-        rmsmin=rms
-        jconmin=jcon
-      endif
-      ENDDO ! K
-!      write (iout,*) "rmsmin",rmsmin," rms",rms
-      call flush(iout)
-      do i=1,2*nres
-        do j=1,3
-          c(j,i)=allcart(j,i,jconmin)
-        enddo
-      enddo
-      return
-      end subroutine closest_coord
-!-----------------------------------------------------------------------------
-! read_coords.F
-!-----------------------------------------------------------------------------
-      subroutine read_coords(ncon,*)
-
-      use energy_data, only: ihpb,jhpb,max_ene
-      use control_data, only: from_bx,from_cx
-      use control, only: tcpu
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-#ifdef MPI
-      use MPI_data
-      include "mpif.h"
-      integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
-!      include "COMMON.MPI"
-#else
-      use MPI_data, only: nprocs
-#endif
-!      include "COMMON.CONTROL"
-!      include "COMMON.CHAIN"
-!      include "COMMON.INTERACT"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.VAR"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.GEO"
-!      include "COMMON.CLUSTER"
-      character(len=3) :: liczba
-      integer :: ncon
-      integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
-        nn,nn1,inan
-      integer :: ixdrf,iret,itmp
-      real(kind=4) :: prec,reini,refree,rmsdev
-      integer :: nrec,nlines,iscor,lenrec,lenrec_in
-      real(kind=8) :: energ,t_acq !,tcpu
-!el      integer ilen,iroof
-!el      external ilen,iroof
-      real(kind=8) :: rjunk
-      integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
-      logical :: lerr
-      real(kind=8) :: energia(0:max_ene),etot
-      real(kind=4) :: csingle(3,2*nres+2)
-      integer :: Previous,Next
-      character(len=256) :: bprotfiles
-!      print *,"Processor",me," calls read_protein_data"
-#ifdef MPI
-      if (me.eq.master) then
-        Previous=MPI_PROC_NULL
-      else
-        Previous=me-1
-      endif
-      if (me.eq.nprocs-1) then
-        Next=MPI_PROC_NULL
-      else
-        Next=me+1
-      endif
-! Set the scratchfile names
-      write (liczba,'(bz,i3.3)') me
-
-      allocate(STATUS(MPI_STATUS_SIZE))
-#endif
-! 1/27/05 AL Change stored coordinates to single precision and don't store 
-!         energy components in the binary databases.
-      lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
-      lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
-#ifdef DEBUG
-      write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
-      write (iout,*) "lenrec_in",lenrec_in
-#endif
-      bprotfiles=scratchdir(:ilen(scratchdir))// &
-             "/"//prefix(:ilen(prefix))//liczba//".xbin"
-! EL
-! allocate cluster arrays
-      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
-      allocate(entfac(maxconf)) !(maxconf)
-      allocate(rmstb(maxconf)) !(maxconf)
-      allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
-      allocate(nss_all(maxstr_proc)) !(maxstr_proc)
-      allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
-      allocate(iscore(maxconf)) !(maxconf)
-
-
-#ifdef CHUJ
-      ICON=1
-  123 continue
-      if (from_cart .and. .not. from_bx .and. .not. from_cx) then
-        if (efree) then
-        read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
-          rmstb(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),rmstb(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),rmstb(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
-!old        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)
-        call read_angles(intin,*13)
-        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
-!
-! CALCULATE DISTANCES
-!
-   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
-      call flush(iout)
-      jj_old=1
-      open (icbase,file=bprotfiles,status="unknown",&
-         form="unformatted",access="direct",recl=lenrec)
-! Read conformations from binary DA files (one per batch) and write them to 
-! a binary DA scratchfile.
-      jj=0
-      jjj=0
-#ifdef MPI
-      write (liczba,'(bz,i3.3)') me
-      IF (ME.EQ.MASTER) THEN
-! Only the master reads the database; it'll send it to the other procs
-! through a ring.
-#endif
-        t_acq = tcpu()
-        icount=0
-
-        if (from_bx) then
-
-          open (intin,file=intinname,status="old",form="unformatted",&
-                  access="direct",recl=lenrec_in)
-
-        else if (from_cx) then
-#if (defined(AIX) && !defined(JUBL))
-          call xdrfopen_(ixdrf,intinname, "r", iret)
-#else
-          call xdrfopen(ixdrf,intinname, "r", iret)
-#endif
-          prec=10000.0
-          write (iout,*) "xdrfopen: iret",iret
-          if (iret.eq.0) then
-            write (iout,*) "Error: coordinate file ",&
-             intinname(:ilen(intinname))," does not exist."
-            call flush(iout)
-#ifdef MPI
-            call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
-            stop
-          endif
-        else
-          write (iout,*) "Error: coordinate format not specified"
-          call flush(iout)
-#ifdef MPI
-          call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
-#else
-          stop
-#endif
-        endif
-
-!#define DEBUG
-#ifdef DEBUG
-        write (iout,*) "Opening file ",intinname(:ilen(intinname))
-        write (iout,*) "lenrec",lenrec_in
-        call flush(iout)
-#endif
-!#undef DEBUG
-!        write (iout,*) "maxconf",maxconf
-        i=0
-        do while (.true.)
-           i=i+1
-!el           if (i.gt.maxconf) then
-!el             write (iout,*) "Error: too many conformations ",&
-!el              "(",maxconf,") maximum."
-!#ifdef MPI
-!el             call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
-!#endif
-!el             stop
-!el           endif
-!          write (iout,*) "i",i
-!          call flush(iout)
-          if (from_bx) then
-            read(intin,err=101,end=101) &
-             ((csingle(l,k),l=1,3),k=1,nres),&
-             ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
-             nss,(ihpb(k),jhpb(k),k=1,nss),&
-             energy(jj+1),&
-             entfac(jj+1),rmstb(jj+1),iscor
-             do j=1,2*nres
-               do k=1,3
-                 c(k,j)=csingle(k,j)
-               enddo
-             enddo
-          else
-#if (defined(AIX) && !defined(JUBL))
-            call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
-            if (iret.eq.0) goto 101
-            call xdrfint_(ixdrf, nss, iret)
-            if (iret.eq.0) goto 101
-            do j=1,nss
-              call xdrfint_(ixdrf, ihpb(j), iret)
-              if (iret.eq.0) goto 101
-              call xdrfint_(ixdrf, jhpb(j), iret)
-              if (iret.eq.0) goto 101
-            enddo
-            call xdrffloat_(ixdrf,reini,iret)
-            if (iret.eq.0) goto 101
-            call xdrffloat_(ixdrf,refree,iret)
-            if (iret.eq.0) goto 101
-            call xdrffloat_(ixdrf,rmsdev,iret)
-            if (iret.eq.0) goto 101
-            call xdrfint_(ixdrf,iscor,iret)
-            if (iret.eq.0) goto 101
-#else
-!            write (iout,*) "calling xdrf3dfcoord"
-            call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
-!            write (iout,*) "iret",iret
-!            call flush(iout)
-            if (iret.eq.0) goto 101
-            call xdrfint(ixdrf, nss, iret)
-!            write (iout,*) "iret",iret
-!            write (iout,*) "nss",nss
-            call flush(iout)
-            if (iret.eq.0) goto 101
-            do k=1,nss
-              call xdrfint(ixdrf, ihpb(k), iret)
-              if (iret.eq.0) goto 101
-              call xdrfint(ixdrf, jhpb(k), iret)
-              if (iret.eq.0) goto 101
-            enddo
-            call xdrffloat(ixdrf,reini,iret)
-            if (iret.eq.0) goto 101
-            call xdrffloat(ixdrf,refree,iret)
-            if (iret.eq.0) goto 101
-            call xdrffloat(ixdrf,rmsdev,iret)
-            if (iret.eq.0) goto 101
-            call xdrfint(ixdrf,iscor,iret)
-            if (iret.eq.0) goto 101
-#endif
-            energy(jj+1)=reini
-            entfac(jj+1)=refree
-            rmstb(jj+1)=rmsdev
-            do k=1,nres
-              do l=1,3
-                c(l,k)=csingle(l,k)
-              enddo
-            enddo
-            do k=nnt,nct
-              do l=1,3
-                c(l,nres+k)=csingle(l,nres+k-nnt+1)
-              enddo
-            enddo
-          endif
-#ifdef DEBUG
-          write (iout,'(5hREAD ,i5,3f15.4,i10)') &
-           jj+1,energy(jj+1),entfac(jj+1),&
-           rmstb(jj+1),iscor
-          write (iout,*) "Conformation",jjj+1,jj+1
-          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)
-          call flush(iout)
-#endif
-          call add_new_cconf(jjj,jj,jj_old,icount,Next)
-        enddo
-  101   continue
-        write (iout,*) i-1," conformations read from DA file ",&
-          intinname(:ilen(intinname))
-        write (iout,*) jj," conformations read so far"
-        if (from_bx) then
-          close(intin)
-        else
-#if (defined(AIX) && !defined(JUBL))
-          call xdrfclose_(ixdrf, iret)
-#else
-          call xdrfclose(ixdrf, iret)
-#endif
-        endif
-#ifdef MPI
-!#ifdef DEBUG  
-        write (iout,*) "jj_old",jj_old," jj",jj
-!#endif
-        call write_and_send_cconf(icount,jj_old,jj,Next)
-        call MPI_Send(0,1,MPI_INTEGER,Next,570,&
-                   MPI_COMM_WORLD,IERROR)
-        jj_old=jj+1
-#else
-        call write_and_send_cconf(icount,jj_old,jj,Next)
-#endif
-        t_acq = tcpu() - t_acq
-#ifdef MPI
-        write (iout,*) "Processor",me,&
-          " time for conformation read/send",t_acq
-      ELSE
-! A worker gets the confs from the master and sends them to its neighbor
-        t_acq = tcpu()
-        call receive_and_pass_cconf(icount,jj_old,jj,&
-          Previous,Next)
-        t_acq = tcpu() - t_acq
-      ENDIF
-#endif
-      ncon=jj
-!      close(icbase)
-      close(intin)
-
-      write(iout,*)"A total of",ncon," conformations read."
-
-      allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
-#ifdef MPI
-! Check if everyone has the same number of conformations
-      call MPI_Allgather(ncon,1,MPI_INTEGER,&
-        ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
-      lerr=.false.
-      do i=0,nprocs-1
-        if (i.ne.me) then
-            if (ncon.ne.ntot_all(i)) then
-              write (iout,*) "Number of conformations at processor",i,&
-               " differs from that at processor",me,&
-               ncon,ntot_all(i)
-              lerr = .true.
-            endif
-        endif
-      enddo
-      if (lerr) then
-        write (iout,*)
-        write (iout,*) "Number of conformations read by processors"
-        write (iout,*)
-        do i=0,nprocs-1
-          write (iout,'(8i10)') i,ntot_all(i)
-        enddo
-        write (iout,*) "Calculation terminated."
-        call flush(iout)
-        return 1
-      endif
-      return
-#endif
- 1111 write(iout,*) "Error opening coordinate file ",&
-       intinname(:ilen(intinname))
-      call flush(iout)
-      return 1
-      end subroutine read_coords
-!------------------------------------------------------------------------------
-      subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
-
-      use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
-      use energy_data, only: itel,itype,dsc,max_ene
-      use control_data, only: symetr
-      use geometry, only: int_from_cart1
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-!      include "COMMON.CLUSTER"
-!      include "COMMON.CONTROL"
-!      include "COMMON.CHAIN"
-!      include "COMMON.INTERACT"
-!      include "COMMON.LOCAL"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.NAMES"
-!      include "COMMON.VAR"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.GEO"
-      integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
-        nn,nn1,inan,Next,itj,chalen
-      real(kind=8) :: etot,energia(0:max_ene)
-      jjj=jjj+1
-      chalen=int((nct-nnt+2)/symetr)
-      call int_from_cart1(.false.)
-      do j=nnt+1,nct
-        if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
-         if (j.gt.2) then
-          if (itel(j).ne.0 .and. itel(j-1).ne.0) then
-          write (iout,*) "Conformation",jjj,jj+1
-          write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
-       chalen
-          write (iout,*) "The Cartesian geometry is:"
-          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
-          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
-          write (iout,*) "The internal geometry is:"
-          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
-          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+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,*) &
-            "This conformation WILL NOT be added to the database."
-          return
-          endif
-         endif
-        endif
-      enddo
-      do j=nnt,nct
-        itj=itype(j)
-        if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) &
-                                  .gt.2.0d0) then
-          write (iout,*) "Conformation",jjj,jj+1
-          write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
-          write (iout,*) "The Cartesian geometry is:"
-          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
-          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
-          write (iout,*) "The internal geometry is:"
-          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
-          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+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,*) &
-            "This conformation WILL NOT be added to the database."
-          return
-        endif
-      enddo
-      do j=3,nres
-        if (theta(j).le.0.0d0) then
-          write (iout,*) &
-            "Zero theta angle(s) in conformation",jjj,jj+1
-          write (iout,*) "The Cartesian geometry is:"
-          write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
-          write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
-          write (iout,*) "The internal geometry is:"
-          write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
-          write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+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,*) &
-            "This conformation WILL NOT be added to the database."
-          return
-        endif
-        if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
-      enddo
-      jj=jj+1
-#ifdef DEBUG
-      write (iout,*) "Conformation",jjj,jj
-      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,'(8f10.4)') (vbld(k),k=nnt+1,nct)
-      write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+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)') (vbld(k+nres),k=nnt,nct)
-      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)
-      write (iout,'(e15.5,16i5)') entfac(icount+1)
-!     &        iscore(icount+1,0)
-#endif
-      icount=icount+1
-      call store_cconf_from_file(jj,icount)
-      if (icount.eq.maxstr_proc) then
-#ifdef DEBUG
-        write (iout,* ) "jj_old",jj_old," jj",jj
-#endif
-        call write_and_send_cconf(icount,jj_old,jj,Next)
-        jj_old=jj+1
-        icount=0
-      endif
-      return
-      end subroutine add_new_cconf
-!------------------------------------------------------------------------------
-      subroutine store_cconf_from_file(jj,icount)
-   
-      use energy_data, only: ihpb,jhpb
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-!      include "COMMON.CLUSTER"
-!      include "COMMON.CHAIN"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.INTERACT"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.VAR"
-      integer :: i,j,jj,icount
-! Store the conformation that has been read in
-      do i=1,2*nres
-        do j=1,3
-          allcart(j,i,icount)=c(j,i)
-        enddo
-      enddo
-      nss_all(icount)=nss
-      do i=1,nss
-        ihpb_all(i,icount)=ihpb(i)
-        jhpb_all(i,icount)=jhpb(i)
-      enddo
-      return
-      end subroutine store_cconf_from_file
-!------------------------------------------------------------------------------
-      subroutine write_and_send_cconf(icount,jj_old,jj,Next)
-
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-#ifdef MPI
-      use MPI_data
-      include "mpif.h"
-      integer :: IERROR
-!      include "COMMON.MPI"
-#endif
-!      include "COMMON.CHAIN"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.INTERACT"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.CLUSTER"
-!      include "COMMON.VAR"
-      integer :: icount,jj_old,jj,Next
-! Write the structures to a scratch file
-#ifdef MPI
-! Master sends the portion of conformations that have been read in to the neighbor
-#ifdef DEBUG
-      write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
-      call flush(iout)
-#endif
-      call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
-          Next,571,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
-          Next,572,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
-          Next,573,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
-          Next,577,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
-          Next,579,MPI_COMM_WORLD,IERROR)
-      call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
-          MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
-#endif
-      call dawrite_ccoords(jj_old,jj,icbase)
-      return
-      end subroutine write_and_send_cconf
-!------------------------------------------------------------------------------
-#ifdef MPI
-      subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
-
-      use MPI_data
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-      include "mpif.h"
-      integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
-!      include "COMMON.MPI"
-!      include "COMMON.CHAIN"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.INTERACT"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.VAR"
-!      include "COMMON.GEO"
-!      include "COMMON.CLUSTER"
-      integer :: i,j,k,icount,jj_old,jj,Previous,Next
-      icount=1
-#ifdef DEBUG
-      write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
-      call flush(iout)
-#endif
-      do while (icount.gt.0) 
-      call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
-           STATUS,IERROR)
-      call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
-           IERROR)
-#ifdef DEBUG
-      write (iout,*) "Processor",me," icount",icount
-#endif
-      if (icount.eq.0) return
-      call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
-          Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
-        Next,571,MPI_COMM_WORLD,IERROR)
-      call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
-          Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
-        Next,572,MPI_COMM_WORLD,IERROR)
-      call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
-          Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
-        Next,573,MPI_COMM_WORLD,IERROR)
-      call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
-        Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
-        Next,577,MPI_COMM_WORLD,IERROR)
-      call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
-        Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
-        Next,579,MPI_COMM_WORLD,IERROR)
-      call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
-        MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
-      call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
-        MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
-      jj=jj_old+icount-1
-      call dawrite_ccoords(jj_old,jj,icbase)
-      jj_old=jj+1
-#ifdef DEBUG
-      write (iout,*) "Processor",me," received",icount," conformations"
-      do i=1,icount
-        write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
-        write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
-        write (iout,'(e15.5,16i5)') entfac(i)
-      enddo
-#endif
-      enddo
-      return
-      end subroutine receive_and_pass_cconf
-#endif
-!------------------------------------------------------------------------------
-      subroutine daread_ccoords(istart_conf,iend_conf)
-
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-#ifdef MPI
-      use MPI_data
-      include "mpif.h"
-!      include "COMMON.MPI"
-#endif
-!      include "COMMON.CHAIN"
-!      include "COMMON.CLUSTER"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.INTERACT"
-!      include "COMMON.VAR"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.GEO"
-      integer :: istart_conf,iend_conf
-      integer :: i,j,ij,ii,iii
-      integer :: len
-      character(len=16) :: form,acc
-      character(len=32) :: nam
-!
-! Read conformations off a DA scratchfile.
-!
-#ifdef DEBUG
-      write (iout,*) "DAREAD_COORDS"
-      write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
-      inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
-      write (iout,*) "len=",len," form=",form," acc=",acc
-      write (iout,*) "nam=",nam
-      call flush(iout)
-#endif
-      do ii=istart_conf,iend_conf
-        ij = ii - istart_conf + 1
-        iii=list_conf(ii)
-#ifdef DEBUG
-        write (iout,*) "Reading binary file, record",iii," ii",ii
-        call flush(iout)
-#endif
-        read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
-          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
-          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
-          entfac(ii),rmstb(ii)
-#ifdef DEBUG
-        write (iout,*) ii,iii,ij,entfac(ii)
-        write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
-        write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
-          i=nnt+nres,nct+nres)
-        write (iout,'(2e15.5)') entfac(ij)
-        write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
-          jhpb_all(i,ij),i=1,nss)
-        call flush(iout)
-#endif
-      enddo
-      return
-      end subroutine daread_ccoords
-!------------------------------------------------------------------------------
-      subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
-
-!      implicit none
-!      include "DIMENSIONS"
-!      include "sizesclu.dat"
-#ifdef MPI
-      use MPI_data
-      include "mpif.h"
-!      include "COMMON.MPI"
-#endif
-!      include "COMMON.CHAIN"
-!      include "COMMON.INTERACT"
-!      include "COMMON.IOUNITS"
-!      include "COMMON.VAR"
-!      include "COMMON.SBRIDGE"
-!      include "COMMON.GEO"
-!      include "COMMON.CLUSTER"
-      integer :: istart_conf,iend_conf
-      integer :: i,j,ii,ij,iii,unit_out
-      integer :: len
-      character(len=16) :: form,acc
-      character(len=32) :: nam
-!
-! Write conformations to a DA scratchfile.
-!
-#ifdef DEBUG
-      write (iout,*) "DAWRITE_COORDS"
-      write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
-      write (iout,*) "lenrec",lenrec
-      inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
-      write (iout,*) "len=",len," form=",form," acc=",acc
-      write (iout,*) "nam=",nam
-      call flush(iout)
-#endif
-      do ii=istart_conf,iend_conf
-        iii=list_conf(ii)
-        ij = ii - istart_conf + 1
-#ifdef DEBUG
-        write (iout,*) "Writing binary file, record",iii," ii",ii
-        call flush(iout)
-#endif
-        write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
-          ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
-          nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
-          entfac(ii),rmstb(ii)
-#ifdef DEBUG
-        write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
-        write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
-         nct+nres)
-        write (iout,'(2e15.5)') entfac(ij)
-        write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
-         nss_all(ij))
-        call flush(iout)
-#endif
-      enddo
-      return
-      end subroutine dawrite_ccoords
-!-----------------------------------------------------------------------------
-! readrtns.F
-!-----------------------------------------------------------------------------
-      subroutine read_control
-!
-! Read molecular data
-!
-      use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
-      use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
-                 iscode,symetr,punch_dist,print_dist,nstart,nend,&
-                 caonly,iopt,efree,lprint_cart,lprint_int
-!      implicit none
-!      include 'DIMENSIONS'
-!      include 'sizesclu.dat'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.TIME1'
-!      include 'COMMON.SBRIDGE'
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.CLUSTER'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.HEADER'
-!      include 'COMMON.FFIELD'
-!      include 'COMMON.FREE'
-!      include 'COMMON.INTERACT'
-      character(len=320) :: controlcard !,ucase
-!#ifdef MPL
-!      include 'COMMON.INFO'
-!#endif
-      integer :: i
-
-      read (INP,'(a80)') titel
-      call card_concat(controlcard,.true.)
-
-      call readi(controlcard,'NRES',nres,0)
-
-!      call alloc_clust_arrays
-      allocate(rcutoff(max_cut+1)) !(max_cut+1)
-      allocate(beta_h(maxT)) !(maxT)
-      allocate(mult(nres)) !(maxres)
-
-
-      call readi(controlcard,'RESCALE',rescale_mode,2)
-      call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
-      write (iout,*) "DISTCHAINMAX",distchainmax
-      call readi(controlcard,'PDBOUT',outpdb,0)
-      call readi(controlcard,'MOL2OUT',outmol2,0)
-      refstr=(index(controlcard,'REFSTR').gt.0)
-      write (iout,*) "REFSTR",refstr
-      pdbref=(index(controlcard,'PDBREF').gt.0)
-      iscode=index(controlcard,'ONE_LETTER')
-      tree=(index(controlcard,'MAKE_TREE').gt.0)
-      min_var=(index(controlcard,'MINVAR').gt.0)
-      plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
-      punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
-      call readi(controlcard,'NCUT',ncut,1)
-      call readi(controlcard,'SYM',symetr,1)
-      write (iout,*) 'sym', symetr
-      call readi(controlcard,'NSTART',nstart,0)
-      call readi(controlcard,'NEND',nend,0)
-      call reada(controlcard,'ECUT',ecut,10.0d0)
-      call reada(controlcard,'PROB',prob_limit,0.99d0)
-      write (iout,*) "Probability limit",prob_limit
-      lgrp=(index(controlcard,'LGRP').gt.0)
-      caonly=(index(controlcard,'CA_ONLY').gt.0)
-      print_dist=(index(controlcard,'PRINT_DIST').gt.0)
-      call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
-      call readi(controlcard,'IOPT',iopt,2)
-      lside = index(controlcard,"SIDE").gt.0
-      efree = index(controlcard,"EFREE").gt.0
-      call readi(controlcard,'NTEMP',nT,1)
-      write (iout,*) "nT",nT
-!el      call reada(controlcard,'TEMP0',temp0,300.0d0) !el
-      call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
-      write (iout,*) "nT",nT
-      write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
-      do i=1,nT
-        beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
-      enddo
-      write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
-      lprint_cart=index(controlcard,"PRINT_CART") .gt.0
-      lprint_int=index(controlcard,"PRINT_INT") .gt.0
-      if (min_var) iopt=1
-      return
-      end subroutine read_control
-!-----------------------------------------------------------------------------
-      subroutine molread
-!
-! Read molecular data.
-!
-      use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
-      use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
-!                 wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
-!                 wturn3,wturn4,wturn6,wvdwpp,weights
-      use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
-                 indpdb
-      use geometry, only: chainbuild,alloc_geo_arrays
-      use energy, only: alloc_ener_arrays
-      use control, only: rescode,setup_var,init_int_table
-      use conform_compar, only: contact
-!      implicit none
-!      include 'DIMENSIONS'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.GEO'
-!      include 'COMMON.VAR'
-!      include 'COMMON.INTERACT'
-!      include 'COMMON.LOCAL'
-!      include 'COMMON.NAMES'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.FFIELD'
-!      include 'COMMON.SBRIDGE'
-!      include 'COMMON.HEADER'
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.CONTACTS'
-!      include 'COMMON.TIME1'
-!#ifdef MPL
-!      include 'COMMON.INFO'
-!#endif
-      character(len=4) ::  sequence(nres) !(maxres)
-      character(len=800) :: weightcard
-!      integer rescode
-      real(kind=8) :: x(6*nres) !(maxvar)
-      integer  :: itype_pdb(nres) !(maxres)
-!      logical seq_comp
-      integer :: i,j,kkk
-!
-! Body
-!
-!el      allocate(weights(n_ene))
-      allocate(weights(max_ene))
-      call alloc_geo_arrays
-      call alloc_ener_arrays
-!-----------------------------
-      allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
-      allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
-      allocate(itype(nres+2)) !(maxres)
-      allocate(itel(nres+2))
-
-      do i=1,2*nres+2
-        do j=1,3
-          c(j,i)=0
-          dc(j,i)=0
-        enddo
-      enddo
-      do i=1,nres+2
-        itype(i)=0
-        itel(i)=0
-      enddo
-!--------------------------
-! Read weights of the subsequent energy terms.
-      call card_concat(weightcard,.true.)
-      call reada(weightcard,'WSC',wsc,1.0d0)
-      call reada(weightcard,'WLONG',wsc,wsc)
-      call reada(weightcard,'WSCP',wscp,1.0d0)
-      call reada(weightcard,'WELEC',welec,1.0D0)
-      call reada(weightcard,'WVDWPP',wvdwpp,welec)
-      call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
-      call reada(weightcard,'WCORR4',wcorr4,0.0D0)
-      call reada(weightcard,'WCORR5',wcorr5,0.0D0)
-      call reada(weightcard,'WCORR6',wcorr6,0.0D0)
-      call reada(weightcard,'WTURN3',wturn3,1.0D0)
-      call reada(weightcard,'WTURN4',wturn4,1.0D0)
-      call reada(weightcard,'WTURN6',wturn6,1.0D0)
-      call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
-      call reada(weightcard,'WSCCOR',wsccor,1.0D0)
-      call reada(weightcard,'WBOND',wbond,1.0D0)
-      call reada(weightcard,'WTOR',wtor,1.0D0)
-      call reada(weightcard,'WTORD',wtor_d,1.0D0)
-      call reada(weightcard,'WANG',wang,1.0D0)
-      call reada(weightcard,'WSCLOC',wscloc,1.0D0)
-      call reada(weightcard,'SCAL14',scal14,0.4D0)
-      call reada(weightcard,'SCALSCP',scalscp,1.0d0)
-      call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
-      call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
-      call reada(weightcard,'TEMP0',temp0,300.0d0)   !!! el
-      if (index(weightcard,'SOFT').gt.0) ipot=6
-! 12/1/95 Added weight for the multi-body term WCORR
-      call reada(weightcard,'WCORRH',wcorr,1.0D0)
-      if (wcorr4.gt.0.0d0) wcorr=wcorr4
-      weights(1)=wsc
-      weights(2)=wscp
-      weights(3)=welec
-      weights(4)=wcorr
-      weights(5)=wcorr5
-      weights(6)=wcorr6
-      weights(7)=wel_loc
-      weights(8)=wturn3
-      weights(9)=wturn4
-      weights(10)=wturn6
-      weights(11)=wang
-      weights(12)=wscloc
-      weights(13)=wtor
-      weights(14)=wtor_d
-      weights(15)=wstrain
-      weights(16)=wvdwpp
-      weights(17)=wbond
-      weights(18)=scal14
-!el      weights(19)=wsccor !!!!!!!!!!!!!!!!
-      weights(21)=wsccor
-      write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
-        wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
-        wturn4,wturn6,wsccor
-   10 format (/'Energy-term weights (unscaled):'// &
-       'WSCC=   ',f10.6,' (SC-SC)'/ &
-       'WSCP=   ',f10.6,' (SC-p)'/ &
-       'WELEC=  ',f10.6,' (p-p electr)'/ &
-       'WVDWPP= ',f10.6,' (p-p VDW)'/ &
-       'WBOND=  ',f10.6,' (stretching)'/ &
-       'WANG=   ',f10.6,' (bending)'/ &
-       'WSCLOC= ',f10.6,' (SC local)'/ &
-       'WTOR=   ',f10.6,' (torsional)'/ &
-       'WTORD=  ',f10.6,' (double torsional)'/ &
-       'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
-       'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
-       'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
-       'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
-       'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
-       'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
-       'WTURN4= ',f10.6,' (turns, 4th order)'/ &
-       'WTURN6= ',f10.6,' (turns, 6th order)'/ &
-       'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
-
-      if (wcorr4.gt.0.0d0) then
-        write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
-         'between contact pairs of peptide groups'
-        write (iout,'(2(a,f5.3/))') &
-        'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
-        'Range of quenching the correlation terms:',2*delt_corr
-      else if (wcorr.gt.0.0d0) then
-        write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
-         'between contact pairs of peptide groups'
-      endif
-      write (iout,'(a,f8.3)') &
-        'Scaling factor of 1,4 SC-p interactions:',scal14
-      write (iout,'(a,f8.3)') &
-        'General scaling factor of SC-p interactions:',scalscp
-      r0_corr=cutoff_corr-delt_corr
-      do i=1,20
-        aad(i,1)=scalscp*aad(i,1)
-        aad(i,2)=scalscp*aad(i,2)
-        bad(i,1)=scalscp*bad(i,1)
-        bad(i,2)=scalscp*bad(i,2)
-      enddo
-
-      call flush(iout)
-      print *,'indpdb=',indpdb,' pdbref=',pdbref
-
-! Read sequence if not taken from the pdb file.
-      if (iscode.gt.0) then
-        read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
-      else
-        read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
-      endif
-! Convert sequence to numeric code
-      do i=1,nres
-        itype(i)=rescode(i,sequence(i),iscode)
-      enddo
-      print *,nres
-      print '(20i4)',(itype(i),i=1,nres)
-
-      do i=1,nres
-#ifdef PROCOR
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
-#else
-        if (itype(i).eq.ntyp1) then
-#endif
-          itel(i)=0
-#ifdef PROCOR
-        else if (iabs(itype(i+1)).ne.20) then
-#else
-        else if (iabs(itype(i)).ne.20) then
-#endif
-          itel(i)=1
-        else
-          itel(i)=2
-        endif
-      enddo
-      write (iout,*) "ITEL"
-      do i=1,nres-1
-        write (iout,*) i,itype(i),itel(i)
-      enddo
-
-      print *,'Call Read_Bridge.'
-      call read_bridge
-      nnt=1
-      nct=nres
-      print *,'NNT=',NNT,' NCT=',NCT
-      if (itype(1).eq.ntyp1) nnt=2
-      if (itype(nres).eq.ntyp1) nct=nct-1
-      if (nstart.lt.nnt) nstart=nnt
-      if (nend.gt.nct .or. nend.eq.0) nend=nct
-      write (iout,*) "nstart",nstart," nend",nend
-      nres0=nres
-!      if (pdbref) then
-!        read(inp,'(a)') pdbfile
-!        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
-!        open(ipdbin,file=pdbfile,status='old',err=33)
-!        goto 34 
-!  33    write (iout,'(a)') 'Error opening PDB file.'
-!        stop
-!  34    continue
-!        print *,'Begin reading pdb data'
-!        call readpdb
-!        print *,'Finished reading pdb data'
-!        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
-!        do i=1,nres
-!          itype_pdb(i)=itype(i)
-!        enddo
-!        close (ipdbin)
-!        write (iout,'(a,i3)') 'nsup=',nsup
-!        nstart_seq=nnt
-!        if (nsup.le.(nct-nnt+1)) then
-!          do i=0,nct-nnt+1-nsup
-!            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
-!              nstart_seq=nnt+i
-!              goto 111
-!            endif
-!          enddo
-!          write (iout,'(a)') 
-!     &            'Error - sequences to be superposed do not match.'
-!          stop
-!        else
-!          do i=0,nsup-(nct-nnt+1)
-!            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
-!     &      then
-!              nstart_sup=nstart_sup+i
-!              nsup=nct-nnt+1
-!              goto 111
-!            endif
-!          enddo 
-!          write (iout,'(a)') 
-!     &            'Error - sequences to be superposed do not match.'
-!        endif
-!  111   continue
-!        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
-!     &                 ' nstart_seq=',nstart_seq
-!      endif
-write(iout,*)"przed ini_int_tab"
-      call init_int_table
-write(iout,*)"po ini_int_tab"
-write(iout,*)"przed setup var"
-      call setup_var
-write(iout,*)"po setup var"
-      write (iout,*) "molread: REFSTR",refstr
-      if (refstr) then
-        if (.not.pdbref) then
-          call read_angles(inp,*38)
-          goto 39
-   38     write (iout,'(a)') 'Error reading reference structure.'
-#ifdef MPL
-          call mp_stopall(Error_Msg)
-#else
-          stop 'Error reading reference structure'
-#endif
-   39     call chainbuild
-          nstart_sup=nnt
-          nstart_seq=nnt
-          nsup=nct-nnt+1
-          kkk=1
-          do i=1,2*nres
-            do j=1,3
-              cref(j,i,kkk)=c(j,i)
-            enddo
-          enddo
-        endif
-        call contact(.true.,ncont_ref,icont_ref)
-      endif
-      return
-      end subroutine molread
-!-----------------------------------------------------------------------------
-      subroutine openunits
-!      implicit none
-!      include 'DIMENSIONS'
-      use control_data, only: from_cx,from_bx,from_cart
-#ifdef MPI
-      use MPI_data
-      include "mpif.h"
-      character(len=3) :: liczba
-!      include "COMMON.MPI"
-#endif
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.CONTROL'
-      integer :: lenpre,lenpot !,ilen
-!      external ilen
-      character(len=16) :: cformat,cprint
-!      character(len=16) ucase
-      integer :: lenint,lenout
-      call getenv('INPUT',prefix)
-      call getenv('OUTPUT',prefout)
-      call getenv('INTIN',prefintin)
-      call getenv('COORD',cformat)
-      call getenv('PRINTCOOR',cprint)
-      call getenv('SCRATCHDIR',scratchdir)
-      from_bx=.true.
-      from_cx=.false.
-      if (index(ucase(cformat),'CX').gt.0) then
-        from_cx=.true.
-        from_bx=.false.
-      endif
-      from_cart=.true.
-      lenpre=ilen(prefix)
-      lenout=ilen(prefout)
-      lenint=ilen(prefintin)
-! Get the names and open the input files
-      open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
-#ifdef MPI
-      write (liczba,'(bz,i3.3)') me
-      outname=prefout(:lenout)//'_clust.out_'//liczba
-#else
-      outname=prefout(:lenout)//'_clust.out'
-#endif
-      if (from_bx) then
-        intinname=prefintin(:lenint)//'.bx'
-      else if (from_cx) then
-        intinname=prefintin(:lenint)//'.cx'
-      else
-        intinname=prefintin(:lenint)//'.int'
-      endif
-      rmsname=prefintin(:lenint)//'.rms'
-      open (jplot,file=prefout(:ilen(prefout))//'.tex',&
-             status='unknown')
-      open (jrms,file=rmsname,status='unknown')
-      open(iout,file=outname,status='unknown')
-! Get parameter filenames and open the parameter files.
-      call getenv('BONDPAR',bondname)
-      open (ibond,file=bondname,status='old')
-      call getenv('THETPAR',thetname)
-      open (ithep,file=thetname,status='old')
-      call getenv('ROTPAR',rotname)
-      open (irotam,file=rotname,status='old')
-      call getenv('TORPAR',torname)
-      open (itorp,file=torname,status='old')
-      call getenv('TORDPAR',tordname)
-      open (itordp,file=tordname,status='old')
-      call getenv('FOURIER',fouriername)
-      open (ifourier,file=fouriername,status='old')
-      call getenv('ELEPAR',elename)
-      open (ielep,file=elename,status='old')
-      call getenv('SIDEPAR',sidename)
-      open (isidep,file=sidename,status='old')
-      call getenv('SIDEP',sidepname)
-      open (isidep1,file=sidepname,status="old")
-      call getenv('SCCORPAR',sccorname)
-      open (isccor,file=sccorname,status="old")
-#ifndef OLDSCP
-!
-! 8/9/01 In the newest version SCp interaction constants are read from a file
-! Use -DOLDSCP to use hard-coded constants instead.
-!
-      call getenv('SCPPAR',scpname)
-      open (iscpp,file=scpname,status='old')
-#endif
-      return
-      end subroutine openunits
-!-----------------------------------------------------------------------------
-! geomout.F
-!-----------------------------------------------------------------------------
-      subroutine pdboutC(etot,rmsd,tytul)
-
-      use energy_data, only: ihpb,jhpb,itype
-!      implicit real*8 (a-h,o-z)
-!      include 'DIMENSIONS'
-!      include 'COMMON.CONTROL'
-!      include 'COMMON.CHAIN'
-!      include 'COMMON.INTERACT'
-!      include 'COMMON.NAMES'
-!      include 'COMMON.IOUNITS'
-!      include 'COMMON.HEADER'
-!      include 'COMMON.SBRIDGE'
-!      include 'COMMON.TEMPFAC'
-      character(len=50) :: tytul
-      character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
-                                               'G','H','I','J'/)
-      integer :: ica(nres)
-      real(kind=8) :: etot,rmsd
-      integer :: iatom,ichain,ires,i,j,iti
-
-      write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
-        ' ENERGY ',etot,' RMS ',rmsd
-      iatom=0
-      ichain=1
-      ires=0
-      do i=nnt,nct
-        iti=itype(i)
-        if (iti.eq.ntyp1) then
-          ichain=ichain+1
-          ires=0
-          write (ipdb,'(a)') 'TER'
-        else
-        ires=ires+1
-        iatom=iatom+1
-        ica(i)=iatom
-        write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
-           ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
-        if (iti.ne.10) then
-          iatom=iatom+1
-          write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
-            ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
-        endif
-        endif
-      enddo
-      write (ipdb,'(a)') 'TER'
-      do i=nnt,nct-1
-        if (itype(i).eq.ntyp1) cycle
-        if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
-          write (ipdb,30) ica(i),ica(i+1)
-        else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
-          write (ipdb,30) ica(i),ica(i+1),ica(i)+1
-        else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
-          write (ipdb,30) ica(i),ica(i)+1
-        endif
-      enddo
-      if (itype(nct).ne.10) then
-        write (ipdb,30) ica(nct),ica(nct)+1
-      endif
-      do i=1,nss
-        write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
-      enddo
-      write (ipdb,'(a6)') 'ENDMDL'
-  10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
-  20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
-  30  FORMAT ('CONECT',8I5)
-      return
-      end subroutine pdboutC
-!-----------------------------------------------------------------------------
-      subroutine cartout(igr,i,etot,free,rmsd,plik)
-!     implicit real*8 (a-h,o-z)
-!     include 'DIMENSIONS'
-!     include 'sizesclu.dat'
-!     include 'COMMON.IOUNITS'
-!     include 'COMMON.CHAIN'
-!     include 'COMMON.VAR'
-!     include 'COMMON.LOCAL'
-!     include 'COMMON.INTERACT'
-!     include 'COMMON.NAMES'
-!     include 'COMMON.GEO'
-!     include 'COMMON.CLUSTER'
-      integer :: igr,i,j,k
-      real(kind=8) :: etot,free,rmsd
-      character(len=80) :: plik
-      open (igeom,file=plik,position='append')
-      write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
-      write (igeom,'(i4,$)') &
-        nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
-      write (igeom,'(i10)') iscore(i)
-      write (igeom,'(8f10.5)') &
-        ((allcart(k,j,i),k=1,3),j=1,nres),&
-        ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
-      return
-      end subroutine cartout
-!------------------------------------------------------------------------------
-!      subroutine alloc_clust_arrays(n_conf)
-
-!      integer :: n_conf
-!COMMON.CLUSTER
-!      common /clu/
-!      allocate(diss(maxdist)) !(maxdist)
-!el      allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
-!      allocatable :: enetb !(max_ene,maxstr_proc)
-!el      allocate(entfac(maxconf)) !(maxconf)
-!      allocatable :: totfree_gr !(maxgr)
-!el      allocate(rcutoff(max_cut+1)) !(max_cut+1)
-!      common /clu1/
-!      allocatable :: licz,iass !(maxgr)
-!      allocatable :: nconf !(maxgr,maxingr)
-!      allocatable :: iass_tot !(maxgr,max_cut)
-!      allocatable :: list_conf !(maxconf)
-!      common /alles/
-!el      allocatable :: allcart !(3,maxres2,maxstr_proc)
-!el      allocate(rmstb(maxconf)) !(maxconf)
-!el      allocate(mult(nres)) !(maxres)
-!el      allocatable :: nss_all !(maxstr_proc)
-!el      allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
-!      allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
-!COMMON.TEMPFAC
-!      common /factemp/
-!      allocatable :: tempfac !(2,maxres)
-!COMMON.FREE
-!      common /free/
-!el      allocate(beta_h(maxT)) !(maxT)
-
-!      end subroutine alloc_clust_arrays
-!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
-      end module io_clust