rename
[unres4.git] / source / unres / io_base.F90
diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90
new file mode 100644 (file)
index 0000000..f86b4dd
--- /dev/null
@@ -0,0 +1,1326 @@
+      module io_base
+!-----------------------------------------------------------------------
+      use names
+      use io_units
+      implicit none
+!-----------------------------------------------------------------------------
+! Max. number of AA residues
+      integer,parameter :: maxres=6000!1200
+! Appr. max. number of interaction sites
+      integer,parameter :: maxres2=2*maxres
+!      parameter (maxres6=6*maxres)
+!      parameter (mmaxres2=(maxres2*(maxres2+1)/2))
+!-----------------------------------------------------------------------------
+! Max. number of S-S bridges
+!      integer,parameter :: maxss=20
+!-----------------------------------------------------------------------------
+! Max. number of derivatives of virtual-bond and side-chain vectors in theta
+! or phi.
+!      integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+      contains
+!-----------------------------------------------------------------------------
+! readrtns_CSA.F
+!-----------------------------------------------------------------------------
+      subroutine read_bridge
+! Read information about disulfide bridges.
+      use geometry_data, only: nres
+      use energy_data
+      use control_data, only:out1file
+      use MPI_data
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+!      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.DBASE'
+!      include 'COMMON.THREAD'
+!      include 'COMMON.TIME1'
+!      include 'COMMON.SETUP'
+!el local variables
+      integer :: i,j,ierror
+
+! Read bridging residues.
+      read (inp,*) ns
+      if (ns.gt.0) then
+        allocate(iss(ns))
+        read (inp,*) (iss(i),i=1,ns)
+      endif
+
+!      print *,'ns=',ns
+      if(me.eq.king.or..not.out1file) &
+        write (iout,*) 'ns=',ns
+      if (ns.gt.0) &
+        write(iout,*) ' iss:',(iss(i),i=1,ns)
+! Check whether the specified bridging residues are cystines.
+      do i=1,ns
+       if (itype(iss(i)).ne.1) then
+         if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
+         'Do you REALLY think that the residue ',&
+          restyp(itype(iss(i))),i,&
+         ' can form a disulfide bridge?!!!'
+         write (*,'(2a,i3,a)') &
+         'Do you REALLY think that the residue ',&
+          restyp(itype(iss(i))),i,&
+         ' can form a disulfide bridge?!!!'
+#ifdef MPI
+        call MPI_Finalize(MPI_COMM_WORLD,ierror)
+         stop
+#endif
+        endif
+      enddo
+! Read preformed bridges.
+      if (ns.gt.0) then
+        read (inp,*) nss
+      if (nss.gt.0) then
+        if(.not.allocated(ihpb)) allocate(ihpb(nss))
+        if(.not.allocated(jhpb)) allocate(jhpb(nss))
+        read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
+
+        if(.not.allocated(dhpb)) allocate(dhpb(nss))
+        if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
+        if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
+        if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
+! el Initialize the bridge array
+        do i=1,nss
+         dhpb(i)=0.0D0
+        enddo
+      endif
+!--------------------
+      if(fg_rank.eq.0) &
+       write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+      if (nss.gt.0) then
+       write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
+        nhpb=nss
+! Check if the residues involved in bridges are in the specified list of
+! bridging residues.
+        do i=1,nss
+          do j=1,i-1
+           if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
+            .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
+             write (iout,'(a,i3,a)') 'Disulfide pair',i,&
+            ' contains residues present in other pairs.'
+             write (*,'(a,i3,a)') 'Disulfide pair',i,&
+            ' contains residues present in other pairs.'
+#ifdef MPI
+             call MPI_Finalize(MPI_COMM_WORLD,ierror)
+              stop 
+#endif
+           endif
+          enddo
+         do j=1,ns
+           if (ihpb(i).eq.iss(j)) goto 10
+          enddo
+          write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+   10     continue
+         do j=1,ns
+           if (jhpb(i).eq.iss(j)) goto 20
+          enddo
+          write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
+   20     continue
+          dhpb(i)=d0cm
+          forcon(i)=akcm
+        enddo
+        do i=1,nss
+          ihpb(i)=ihpb(i)+nres
+          jhpb(i)=jhpb(i)+nres
+        enddo
+      endif
+      endif
+!      write(iout,*) "end read_bridge"
+      return
+      end subroutine read_bridge
+!-----------------------------------------------------------------------------
+      subroutine read_x(kanal,*)
+
+     use geometry_data
+     use energy_data
+     use geometry, only:int_from_cart1
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.GEO'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.INTERACT'
+! Read coordinates from input
+!
+!el local variables
+      integer :: l,k,j,i,kanal
+
+      read(kanal,'(8f10.5)',end=10,err=10) &
+        ((c(l,k),l=1,3),k=1,nres),&
+        ((c(l,k+nres),l=1,3),k=nnt,nct)
+      do j=1,3
+        c(j,nres+1)=c(j,1)
+        c(j,2*nres)=c(j,nres)
+      enddo
+      call int_from_cart1(.false.)
+      do i=1,nres-1
+        do j=1,3
+          dc(j,i)=c(j,i+1)-c(j,i)
+          dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+        enddo
+      enddo
+      do i=nnt,nct
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+          do j=1,3
+            dc(j,i+nres)=c(j,i+nres)-c(j,i)
+            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+          enddo
+        endif
+      enddo
+
+      return
+   10 return 1
+      end subroutine read_x
+!-----------------------------------------------------------------------------
+      subroutine read_threadbase
+
+      use geometry_data
+      use energy_data
+      use compare_data
+!      implicit real*8 (a-h,o-z)
+!      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.DBASE'
+!      include 'COMMON.THREAD'
+!      include 'COMMON.TIME1'
+
+!el local variables
+      integer :: k,j,i
+
+! Read pattern database for threading.
+      read (icbase,*) nseq
+      allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
+      allocate(nres_base(3,nseq)) !(3,maxseq)
+      allocate(str_nam(nseq)) !(maxseq)
+      do i=1,nseq
+        read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
+         nres_base(2,i),nres_base(3,i)
+        read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
+         nres_base(1,i))
+!       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
+!    &   nres_base(2,i),nres_base(3,i)
+!       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
+!    &   nres_base(1,i))
+      enddo
+      close (icbase)
+      if (weidis.eq.0.0D0) weidis=0.1D0
+      do i=nnt,nct
+        do j=i+2,nct
+          nhpb=nhpb+1
+          ihpb(nhpb)=i
+          jhpb(nhpb)=j
+          forcon(nhpb)=weidis
+        enddo
+      enddo 
+      read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
+      write (iout,'(a,i5)') 'nexcl: ',nexcl
+      write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
+      return
+      end subroutine read_threadbase
+!-----------------------------------------------------------------------------
+#ifdef WHAM_RUN
+!el      subroutine read_angles(kanal,iscor,energ,iprot,*)
+      subroutine read_angles(kanal,*)
+
+      use geometry_data
+      use energy_data
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'DIMENSIONS.ZSCOPT'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.GEO'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.IOUNITS'
+      character(len=80) :: lineh
+      integer :: iscor,iprot,ic
+      real(kind=8) :: energ
+#else
+      subroutine read_angles(kanal,*)
+
+      use geometry_data
+   !  use energy
+   !  use control
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.GEO'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CONTROL'
+#endif
+! Read angles from input 
+!
+!el local variables
+      integer :: i,kanal
+#ifdef WHAM_RUN
+      read(kanal,'(a80)',end=10,err=10) lineh
+      read(lineh(:5),*,err=8) ic
+      read(lineh(6:),*,err=8) energ
+      goto 9
+    8 ic=1
+      print *,'error, assuming e=1d10',lineh
+      energ=1d10
+      nss=0
+    9 continue
+      read(lineh(18:),*,end=10,err=10) nss
+      IF (NSS.LT.9) THEN
+        read (lineh(20:),*,end=10,err=10) &
+        (IHPB(I),JHPB(I),I=1,NSS),iscor
+      ELSE
+        read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
+        read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
+          I=9,NSS),iscor
+      ENDIF
+!      print *,"energy",energ," iscor",iscor
+#endif
+       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
+       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
+       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
+       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
+
+       do i=1,nres
+! 9/7/01 avoid 180 deg valence angle
+        if (theta(i).gt.179.99d0) theta(i)=179.99d0
+!
+        theta(i)=deg2rad*theta(i)
+        phi(i)=deg2rad*phi(i)
+        alph(i)=deg2rad*alph(i)
+        omeg(i)=deg2rad*omeg(i)
+       enddo
+      return
+   10 return 1
+      end subroutine read_angles
+!-----------------------------------------------------------------------------
+      subroutine reada(rekord,lancuch,wartosc,default)
+
+!      implicit none
+      character*(*) :: rekord,lancuch
+      real(kind=8) :: wartosc,default
+      integer :: iread !,ilen
+!el      external ilen
+      iread=index(rekord,lancuch)
+      if (iread.eq.0) then
+        wartosc=default 
+        return
+      endif   
+      iread=iread+ilen(lancuch)+1
+      read (rekord(iread:),*,err=10,end=10) wartosc
+      return
+  10  wartosc=default
+      return
+      end subroutine reada
+!-----------------------------------------------------------------------------
+      subroutine readi(rekord,lancuch,wartosc,default)
+
+!      implicit none
+      character*(*) :: rekord,lancuch
+      integer :: wartosc,default
+      integer :: iread !,ilen
+!el      external ilen
+      iread=index(rekord,lancuch)
+      if (iread.eq.0) then
+        wartosc=default 
+        return
+      endif   
+      iread=iread+ilen(lancuch)+1
+      read (rekord(iread:),*,err=10,end=10) wartosc
+      return
+  10  wartosc=default
+      return
+      end subroutine readi
+!-----------------------------------------------------------------------------
+      subroutine multreadi(rekord,lancuch,tablica,dim,default)
+
+!      implicit none
+      integer :: dim,i
+      integer :: tablica(dim),default
+      character*(*) :: rekord,lancuch
+      character(len=80) :: aux
+      integer :: iread !,ilen
+!el      external ilen
+      do i=1,dim
+        tablica(i)=default
+      enddo
+      iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+      if (iread.eq.0) return
+      iread=iread+ilen(lancuch)+1
+      read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+   10 return
+      end subroutine multreadi
+!-----------------------------------------------------------------------------
+      subroutine multreada(rekord,lancuch,tablica,dim,default)
+
+!      implicit none
+      integer :: dim,i
+      real(kind=8) :: tablica(dim),default
+      character*(*) :: rekord,lancuch
+      character(len=80) :: aux
+      integer :: iread !,ilen
+!el      external ilen
+      do i=1,dim
+        tablica(i)=default
+      enddo
+      iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+      if (iread.eq.0) return
+      iread=iread+ilen(lancuch)+1
+      read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+   10 return
+      end subroutine multreada
+!-----------------------------------------------------------------------------
+      subroutine card_concat(card,to_upper)
+
+! dla UNRESA to_upper jest zawsze .true.
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.IOUNITS'
+      character(*) :: card
+      character(len=80) :: karta       !,ucase
+      logical :: to_upper
+!el      external ilen
+      read (inp,'(a)') karta
+      if (to_upper) karta=ucase(karta)
+      card=' '
+      do while (karta(80:80).eq.'&')
+        card=card(:ilen(card)+1)//karta(:79)
+        read (inp,'(a)') karta
+        if (to_upper) karta=ucase(karta)
+      enddo
+      card=card(:ilen(card)+1)//karta
+      return
+      end subroutine card_concat
+!-----------------------------------------------------------------------------
+      subroutine read_dist_constr
+      use MPI_data
+   !  use control
+      use geometry, only: dist
+      use geometry_data
+      use control_data
+      use energy_data
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+!      include 'COMMON.SETUP'
+!      include 'COMMON.CONTROL'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.SBRIDGE'
+      integer,dimension(2,100) :: ifrag_,ipair_
+      real(kind=8),dimension(100) :: wfrag_,wpair_
+      character(len=640) :: controlcard
+
+!el local variables
+      integer :: i,k,j,ddjk,ii,jj,itemp
+      integer :: nfrag_,npair_,ndist_
+      real(kind=8) :: dist_cut
+
+!      write (iout,*) "Calling read_dist_constr"
+!      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+!      call flush(iout)
+      call card_concat(controlcard,.true.)
+      call readi(controlcard,"NFRAG",nfrag_,0)
+      call readi(controlcard,"NPAIR",npair_,0)
+      call readi(controlcard,"NDIST",ndist_,0)
+      call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+      call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+      call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+      call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+      call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+!      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+!      write (iout,*) "IFRAG"
+!      do i=1,nfrag_
+!        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+!      enddo
+!      write (iout,*) "IPAIR"
+!      do i=1,npair_
+!        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+!      enddo
+      if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
+      if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
+      if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
+      if(.not.allocated(forcon)) allocate(forcon(maxdim))
+
+      call flush(iout)
+      do i=1,nfrag_
+        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+        if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
+          ifrag_(2,i)=nstart_sup+nsup-1
+!        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+        call flush(iout)
+        if (wfrag_(i).gt.0.0d0) then
+        do j=ifrag_(1,i),ifrag_(2,i)-1
+          do k=j+1,ifrag_(2,i)
+!            write (iout,*) "j",j," k",k
+            ddjk=dist(j,k)
+            if (constr_dist.eq.1) then
+              nhpb=nhpb+1
+              ihpb(nhpb)=j
+              jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+              forcon(nhpb)=wfrag_(i) 
+            else if (constr_dist.eq.2) then
+              if (ddjk.le.dist_cut) then
+                nhpb=nhpb+1
+                ihpb(nhpb)=j
+                jhpb(nhpb)=k
+                dhpb(nhpb)=ddjk
+                forcon(nhpb)=wfrag_(i) 
+              endif
+            else
+              nhpb=nhpb+1
+              ihpb(nhpb)=j
+              jhpb(nhpb)=k
+              dhpb(nhpb)=ddjk
+              forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+            endif
+#ifdef MPI
+            if (.not.out1file .or. me.eq.king) &
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+             nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+            write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
+             nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+          enddo
+        enddo
+        endif
+      enddo
+      do i=1,npair_
+        if (wpair_(i).gt.0.0d0) then
+        ii = ipair_(1,i)
+        jj = ipair_(2,i)
+        if (ii.gt.jj) then
+          itemp=ii
+          ii=jj
+          jj=itemp
+        endif
+        do j=ifrag_(1,ii),ifrag_(2,ii)
+          do k=ifrag_(1,jj),ifrag_(2,jj)
+            nhpb=nhpb+1
+            ihpb(nhpb)=j
+            jhpb(nhpb)=k
+            forcon(nhpb)=wpair_(i)
+            dhpb(nhpb)=dist(j,k)
+#ifdef MPI
+            if (.not.out1file .or. me.eq.king) &
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+             nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+            write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+             nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+          enddo
+        enddo
+        endif
+      enddo 
+      do i=1,ndist_
+        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+        if (forcon(nhpb+1).gt.0.0d0) then
+          nhpb=nhpb+1
+          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+#ifdef MPI
+          if (.not.out1file .or. me.eq.king) &
+          write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#else
+          write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
+           nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+#endif
+        endif
+      enddo
+      call flush(iout)
+      return
+      end subroutine read_dist_constr
+!-----------------------------------------------------------------------------
+#ifdef WINIFL
+      subroutine flush(iu)
+      return
+      end subroutine flush
+#endif
+#ifdef AIX
+      subroutine flush(iu)
+      call flush_(iu)
+      return
+      end subroutine flush
+#endif
+!-----------------------------------------------------------------------------
+      subroutine copy_to_tmp(source)
+
+!      include "DIMENSIONS"
+!      include "COMMON.IOUNITS"
+      character*(*) :: source
+      character(len=256) :: tmpfile
+!      integer ilen
+!el      external ilen
+      logical :: ex
+      tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
+      inquire(file=tmpfile,exist=ex)
+      if (ex) then
+        write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
+         " to temporary directory..."
+        write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
+        call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
+      endif
+      return
+      end subroutine copy_to_tmp
+!-----------------------------------------------------------------------------
+      subroutine move_from_tmp(source)
+
+!      include "DIMENSIONS"
+!      include "COMMON.IOUNITS"
+      character*(*) :: source
+!      integer ilen
+!el      external ilen
+      write (*,*) "Moving ",source(:ilen(source)),&
+       " from temporary directory to working directory"
+      write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
+      call system("/bin/mv "//source(:ilen(source))//" "//curdir)
+      return
+      end subroutine move_from_tmp
+!-----------------------------------------------------------------------------
+! misc.f
+!-----------------------------------------------------------------------------
+! $Date: 1994/10/12 17:24:21 $
+! $Revision: 2.5 $
+
+      logical function find_arg(ipos,line,errflag)
+
+      integer, parameter :: maxlen = 80
+      character(len=80) :: line
+      character(len=1) :: empty=' ',equal='='
+      logical :: errflag
+      integer :: ipos
+! This function returns .TRUE., if an argument follows keyword keywd; if so
+! IPOS will point to the first non-blank character of the argument. Returns
+! .FALSE., if no argument follows the keyword; in this case IPOS points
+! to the first non-blank character of the next keyword.
+
+      do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
+        ipos=ipos+1
+      enddo 
+      errflag=.false.
+      if (line(ipos:ipos).eq.equal) then
+         find_arg=.true.
+         ipos=ipos+1
+         do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
+           ipos=ipos+1
+         enddo
+         if (ipos.gt.maxlen) errflag=.true.
+      else
+         find_arg=.false.
+      endif
+
+      return
+      end function find_arg
+!-----------------------------------------------------------------------------
+      logical function find_group(iunit,jout,key1)
+
+      character*(*) :: key1
+      character(len=80) :: karta       !,ucase
+      integer :: iunit,jout
+      integer :: ll !,ilen
+!EL      external ilen
+!EL      logical lcom
+      rewind (iunit)
+      karta=' '
+      ll=ilen(key1)
+      do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta)) 
+        read (iunit,'(a)',end=10) karta
+      enddo
+      write (jout,'(2a)') '> ',karta(1:78)
+      find_group=.true.
+      return
+   10 find_group=.false.
+      return
+      end function find_group
+!-----------------------------------------------------------------------------
+      logical function iblnk(charc)
+      character(len=1) :: charc
+      integer :: n
+      n = ichar(charc)
+      iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
+      return
+      end function iblnk
+!-----------------------------------------------------------------------------
+      integer function ilen(string)
+      character*(*) ::  string
+!EL      logical :: iblnk
+
+      ilen = len(string)
+1     if ( ilen .gt. 0 ) then
+         if ( iblnk( string(ilen:ilen) ) ) then
+            ilen = ilen - 1
+            goto 1
+         endif
+      endif
+      return
+      end function ilen
+!-----------------------------------------------------------------------------
+      integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
+      integer :: nkey,i,ikey,narg
+      character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
+!      character(len=16) :: ucase
+
+      do i=1,narg
+        if (ucase(keywd).eq.keywdset(i,ikey)) then
+! Match found
+          in_keywd_set=i
+          return
+        endif
+      enddo
+! No match to the allowed set of keywords if this point is reached. 
+      in_keywd_set=0
+      return
+      end function in_keywd_set
+!-----------------------------------------------------------------------------
+      character function lcase(string)
+      integer :: i, k, idiff
+      character*(*) :: string
+      character(len=1) :: c
+      character(len=40) :: chtmp
+!
+      i = len(lcase)
+      k = len(string)
+      if (i .lt. k) then
+         k = i
+         if (string(k+1:) .ne. ' ') then
+            chtmp = string
+         endif
+      endif
+      idiff = ichar('a') - ichar('A')
+      lcase = string
+      do 99 i = 1, k
+         c = string(i:i)
+         if (lge(c,'A') .and. lle(c,'Z')) then
+            lcase(i:i) = char(ichar(c) + idiff)
+         endif
+   99 continue
+      return
+      end function lcase
+!-----------------------------------------------------------------------------
+      logical function lcom(ipos,karta)
+      character(len=80) :: karta
+      character :: koment(2) = (/'!','#'/)
+      integer :: ipos,i
+
+      lcom=.false.
+      do i=1,2
+        if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
+      enddo
+      return
+      end function lcom
+!-----------------------------------------------------------------------------
+      logical function lower_case(ch)
+      character*(*) :: ch
+      lower_case=(ch.ge.'a' .and. ch.le.'z')
+      return
+      end function lower_case
+!-----------------------------------------------------------------------------
+      subroutine mykey(line,keywd,ipos,blankline,errflag)
+
+! This subroutine seeks a non-empty substring keywd in the string LINE.
+! The substring begins with the first character different from blank and
+! "=" encountered right to the pointer IPOS (inclusively) and terminates
+! at the character left to the first blank or "=". When the subroutine is 
+! exited, the pointer IPOS is moved to the position of the terminator in LINE. 
+! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
+! only separators or the maximum length of the data line (80) has been reached.
+! The logical variable ERRFLAG is set at .TRUE. if the string 
+! consists only from a "=".
+      integer, parameter :: maxlen=80
+      character(len=1) :: empty=' ',equal='=',comma=','
+      character*(*) :: keywd
+      character(len=80) :: line
+      logical :: blankline,errflag     !EL,lcom
+      integer :: ipos,istart,iend
+      errflag=.false.
+      do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
+        ipos=ipos+1
+      enddo
+      if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
+! At this point the rest of the input line turned out to contain only blanks
+! or to be commented out.
+        blankline=.true.
+        return
+      endif
+      blankline=.false.
+      istart=ipos
+! Checks whether the current char is a separator.
+      do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
+       .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
+        ipos=ipos+1
+      enddo
+      iend=ipos-1
+! Error flag set to .true., if the length of the keyword was found less than 1.
+      if (iend.lt.istart) then
+        errflag=.true.
+        return
+      endif
+      keywd=line(istart:iend)
+      return
+      end subroutine  mykey
+!-----------------------------------------------------------------------------
+      subroutine numstr(inum,numm)
+      character(len=10) :: huj='0123456789'
+      character*(*) :: numm
+      integer :: inumm,inum,inum1,inum2
+      inumm=inum
+      inum1=inumm/10
+      inum2=inumm-10*inum1
+      inumm=inum1
+      numm(3:3)=huj(inum2+1:inum2+1)
+      inum1=inumm/10
+      inum2=inumm-10*inum1
+      inumm=inum1
+      numm(2:2)=huj(inum2+1:inum2+1)
+      inum1=inumm/10
+      inum2=inumm-10*inum1 
+      inumm=inum1
+      numm(1:1)=huj(inum2+1:inum2+1)
+      return
+      end subroutine numstr
+!-----------------------------------------------------------------------------
+      function ucase(string)
+      integer :: i, k, idiff
+      character(*) :: string
+      character(len=len(string)) :: ucase
+      character(len=1) :: c
+      character(len=40) :: chtmp
+!
+      i = len(ucase)
+      k = len(string)
+      if (i .lt. k) then
+         k = i
+         if (string(k+1:) .ne. ' ') then
+            chtmp = string
+         endif
+      endif
+      idiff = ichar('a') - ichar('A')
+      ucase = string
+      do 99 i = 1, k
+         c = string(i:i)
+         if (lge(c,'a') .and. lle(c,'z')) then
+            ucase(i:i) = char(ichar(c) - idiff)
+         endif
+   99 continue
+      return
+      end function ucase
+!-----------------------------------------------------------------------------
+! geomout.F
+!-----------------------------------------------------------------------------
+      subroutine pdbout(etot,tytul,iunit)
+
+      use geometry_data, only: c,nres
+      use energy_data
+   !  use control
+      use compare_data
+      use MD_data
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.HEADER'
+!      include 'COMMON.SBRIDGE'
+!      include 'COMMON.DISTFIT'
+!      include 'COMMON.MD'
+!el      character(len=50) :: tytul
+      character*(*) :: tytul
+      character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
+      integer,dimension(nres) :: ica   !(maxres)
+
+!el  local variables
+      integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
+      real(kind=8) :: etot
+      integer :: nres2
+      nres2=2*nres
+
+      if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
+
+      write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
+!model      write (iunit,'(a5,i6)') 'MODEL',1
+      if (nhfrag.gt.0) then
+       do j=1,nhfrag
+        iti=itype(hfrag(1,j))
+        itj=itype(hfrag(2,j))
+        if (j.lt.10) then
+           write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
+                 'HELIX',j,'H',j,&
+                 restyp(iti),hfrag(1,j)-1,&
+                 restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+        else
+             write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
+                 'HELIX',j,'H',j,&
+                 restyp(iti),hfrag(1,j)-1,&
+                 restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
+        endif
+       enddo
+      endif
+
+      if (nbfrag.gt.0) then
+
+       do j=1,nbfrag
+
+        iti=itype(bfrag(1,j))
+        itj=itype(bfrag(2,j)-1)
+
+        write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
+                 'SHEET',1,'B',j,2,&
+                 restyp(iti),bfrag(1,j)-1,&
+                 restyp(itj),bfrag(2,j)-2,0
+
+        if (bfrag(3,j).gt.bfrag(4,j)) then
+
+         itk=itype(bfrag(3,j))
+         itl=itype(bfrag(4,j)+1)
+
+         write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
+                 'SHEET',2,'B',j,2,&
+                 restyp(itl),bfrag(4,j),&
+                 restyp(itk),bfrag(3,j)-1,-1,&
+                 "N",restyp(itk),bfrag(3,j)-1,&
+                 "O",restyp(iti),bfrag(1,j)-1
+
+        else
+
+         itk=itype(bfrag(3,j))
+         itl=itype(bfrag(4,j)-1)
+
+
+        write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
+                 'SHEET',2,'B',j,2,&
+                 restyp(itk),bfrag(3,j)-1,&
+                 restyp(itl),bfrag(4,j)-2,1,&
+                 "N",restyp(itk),bfrag(3,j)-1,&
+                 "O",restyp(iti),bfrag(1,j)-1
+
+
+
+        endif
+         
+       enddo
+      endif 
+
+      if (nss.gt.0) then
+        do i=1,nss
+         if (dyn_ss) then
+          write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
+               'SSBOND',i,'CYS',idssb(i)-nnt+1,&
+                          'CYS',jdssb(i)-nnt+1
+         else
+          write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
+               'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
+                          'CYS',jhpb(i)-nnt+1-nres
+         endif
+        enddo
+      endif
+      
+      iatom=0
+      ichain=1
+      ires=0
+      do i=nnt,nct
+        iti=itype(i)
+        if (iti.eq.ntyp1) then
+          ichain=ichain+1
+          ires=0
+          write (iunit,'(a)') 'TER'
+        else
+        ires=ires+1
+        iatom=iatom+1
+        ica(i)=iatom
+        write (iunit,10) iatom,restyp(iti),chainid(ichain),&
+           ires,(c(j,i),j=1,3),vtot(i)
+        if (iti.ne.10) then
+          iatom=iatom+1
+          write (iunit,20) iatom,restyp(iti),chainid(ichain),&
+            ires,(c(j,nres+i),j=1,3),&
+            vtot(i+nres)
+        endif
+        endif
+      enddo
+      write (iunit,'(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 (iunit,30) ica(i),ica(i+1)
+        else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+          write (iunit,30) ica(i),ica(i+1),ica(i)+1
+        else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+          write (iunit,30) ica(i),ica(i)+1
+        endif
+      enddo
+      if (itype(nct).ne.10) then
+        write (iunit,30) ica(nct),ica(nct)+1
+      endif
+      do i=1,nss
+       if (dyn_ss) then
+        write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
+       else
+        write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+       endif
+      enddo
+      write (iunit,'(a6)') 'ENDMDL'     
+  10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+  20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+  30  FORMAT ('CONECT',8I5)
+      return
+      end subroutine pdbout
+!-----------------------------------------------------------------------------
+      subroutine MOL2out(etot,tytul)
+! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2 
+! format.
+      use geometry_data, only: c
+      use energy_data
+   !  use control
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.HEADER'
+!      include 'COMMON.SBRIDGE'
+      character(len=32) :: tytul,fd
+      character(len=3) :: zahl
+      character(len=6) :: res_num,pom  !,ucase
+
+!el  local variables
+      integer :: i,j
+      real(kind=8) :: etot
+
+#ifdef AIX
+      call fdate_(fd)
+#elif (defined CRAY)
+      call date(fd)
+#else
+      call fdate(fd)
+#endif
+      write (imol2,'(a)') '#'
+      write (imol2,'(a)') &
+       '#         Creating user name:           unres'
+      write (imol2,'(2a)') '#         Creation time:                ',&
+       fd
+      write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
+      write (imol2,'(a)') tytul
+      write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
+      write (imol2,'(a)') 'SMALL'
+      write (imol2,'(a)') 'USER_CHARGES'
+      write (imol2,'(a)') '\@<TRIPOS>ATOM' 
+      do i=nnt,nct
+        write (zahl,'(i3)') i
+        pom=ucase(restyp(itype(i)))
+        res_num = pom(:3)//zahl(2:)
+        write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
+      enddo
+      write (imol2,'(a)') '\@<TRIPOS>BOND'
+      do i=nnt,nct-1
+        write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
+      enddo
+      do i=1,nss
+        write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
+      enddo
+      write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
+      do i=nnt,nct
+        write (zahl,'(i3)') i
+        pom = ucase(restyp(itype(i)))
+        res_num = pom(:3)//zahl(2:)
+        write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
+      enddo
+  10  FORMAT (I7,' CA      ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
+  30  FORMAT (I7,1x,A,I14,' RESIDUE',I13,' ****  ****')
+      return
+      end subroutine MOL2out
+!-----------------------------------------------------------------------------
+      subroutine intout
+
+      use geometry_data
+      use energy_data, only: itype
+   !  use control
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.VAR'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.GEO'
+!      include 'COMMON.TORSION'
+!el  local variables
+      integer :: i,iti
+
+      write (iout,'(/a)') 'Geometry of the virtual chain.'
+      write (iout,'(7a)') '  Res  ','         d','     Theta',&
+       '       Phi','       Dsc','     Alpha','      Omega'
+      do i=1,nres
+       iti=itype(i)
+        write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
+           rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
+           rad2deg*omeg(i)
+      enddo
+      return
+      end subroutine intout
+!-----------------------------------------------------------------------------
+#ifdef CLUSTER
+      subroutine briefout(it,ener,free)!,plik)
+#else
+      subroutine briefout(it,ener)
+#endif
+      use geometry_data
+      use energy_data
+   !  use control
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.VAR'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.GEO'
+!      include 'COMMON.SBRIDGE'
+!     print '(a,i5)',intname,igeom
+!el  local variables
+      integer :: i,it
+      real(kind=8) :: ener,free
+!     character(len=80) :: plik
+!      integer :: iii
+
+#if defined(AIX) || defined(PGI)
+      open (igeom,file=intname,position='append')
+#else
+      open (igeom,file=intname,access='append')
+#endif
+#ifdef WHAM_RUN
+!      iii=igeom
+      igeom=iout
+#endif
+      IF (NSS.LE.9) THEN
+#ifdef CLUSTER
+        WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+      ELSE
+        WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
+#else
+        WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
+      ELSE
+        WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
+#endif
+        WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
+      ENDIF
+!     IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+      WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
+      WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
+!     if (nvar.gt.nphi+ntheta) then
+        write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
+        write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
+!     endif
+      close(igeom)
+  180 format (I5,F12.3,I2,9(1X,2I3))
+  190 format (3X,11(1X,2I3))
+  200 format (8F10.4)
+      return
+      end subroutine briefout
+!-----------------------------------------------------------------------------
+#ifdef WINIFL
+      subroutine fdate(fd)
+      character(len=32) :: fd
+      write(fd,'(32x)')
+      return
+      end subroutine fdate
+#endif
+!-----------------------------------------------------------------------------
+#ifdef WHAM_RUN
+      real(kind=8) function gyrate(jcon)
+#else
+      real(kind=8) function gyrate()
+#endif
+
+      use geometry_data, only: c
+   !  use geometry
+      use energy_data
+!      implicit real*8 (a-h,o-z)
+!      include 'DIMENSIONS'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.CHAIN'
+      real(kind=8) :: rg
+      real(kind=8),dimension(3) :: cen
+!el  local variables
+      integer :: i,j,jcon
+
+      do j=1,3
+       cen(j)=0.0d0
+      enddo
+
+      do i=nnt,nct
+          do j=1,3
+            cen(j)=cen(j)+c(j,i)
+          enddo
+      enddo
+      do j=1,3
+            cen(j)=cen(j)/dble(nct-nnt+1)
+      enddo
+      rg = 0.0d0
+      do i = nnt, nct
+        do j=1,3
+         rg = rg + (c(j,i)-cen(j))**2 
+        enddo
+      end do
+#ifdef WHAM_RUN
+      gyrate = dsqrt(rg/dble(nct-nnt+1))
+#else
+      gyrate = sqrt(rg/dble(nct-nnt+1))
+#endif
+      return
+      end function gyrate
+#ifdef WHAM_RUN
+!-----------------------------------------------------------------------------
+! readrtns.F WHAM
+      subroutine reads(rekord,lancuch,wartosc,default)
+!      implicit none
+      character*(*) :: rekord,lancuch,wartosc,default
+      character(len=80) :: aux
+      integer :: lenlan,lenrec,iread,ireade
+!el      external ilen
+!el      logical iblnk
+!el      external iblnk
+      lenlan=ilen(lancuch)
+      lenrec=ilen(rekord)
+      iread=index(rekord,lancuch(:lenlan)//"=")
+!      print *,"rekord",rekord," lancuch",lancuch
+!      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
+      if (iread.eq.0) then
+        wartosc=default
+        return
+      endif
+      iread=iread+lenlan+1
+!      print *,"iread",iread
+!      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+      do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
+        iread=iread+1
+!      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+      enddo
+!      print *,"iread",iread
+      if (iread.gt.lenrec) then
+         wartosc=default
+        return
+      endif
+      ireade=iread+1
+!      print *,"ireade",ireade
+      do while (ireade.lt.lenrec .and. &
+         .not.iblnk(rekord(ireade:ireade)))
+        ireade=ireade+1
+      enddo
+      wartosc=rekord(iread:ireade)
+      return
+      end subroutine reads
+#endif
+!-----------------------------------------------------------------------------
+! permut.F
+!-----------------------------------------------------------------------------
+      subroutine permut(isym)
+
+      use geometry_data, only: tabperm
+!      implicit real*8 (a-h,o-z) 
+!      include 'DIMENSIONS'
+!      include 'COMMON.LOCAL'
+!      include 'COMMON.VAR'
+!      include 'COMMON.CHAIN'
+!      include 'COMMON.INTERACT'
+!      include 'COMMON.IOUNITS'
+!      include 'COMMON.GEO'
+!      include 'COMMON.NAMES'
+!      include 'COMMON.CONTROL'
+
+      integer :: n,isym
+!      logical nextp
+!el      external nextp
+      integer,dimension(isym) :: a
+!      parameter(n=symetr)
+!el local variables
+      integer :: kkk,i
+
+      n=isym
+      if (n.eq.1) then
+        tabperm(1,1)=1
+        return
+      endif
+      kkk=0
+      do i=1,n
+      a(i)=i
+      enddo
+   10 print *,(a(i),i=1,n)
+      kkk=kkk+1
+      do i=1,n
+      tabperm(kkk,i)=a(i)
+!      write (iout,*) "tututu", kkk
+      enddo
+      if(nextp(n,a)) go to 10
+      return
+      end subroutine permut
+!-----------------------------------------------------------------------------
+      logical function nextp(n,a)
+
+      integer :: n,i,j,k,t
+!      logical :: nextp
+      integer,dimension(n) :: a
+      i=n-1
+   10 if(a(i).lt.a(i+1)) go to 20
+      i=i-1
+      if(i.eq.0) go to 20
+      go to 10
+   20 j=i+1
+      k=n
+   30 t=a(j)
+      a(j)=a(k)
+      a(k)=t
+      j=j+1
+      k=k-1
+      if(j.lt.k) go to 30
+      j=i
+      if(j.ne.0) go to 40
+      nextp=.false.
+      return
+   40 j=j+1
+      if(a(j).lt.a(i)) go to 40
+      t=a(i)
+      a(i)=a(j)
+      a(j)=t
+      nextp=.true.
+      return
+      end function nextp
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+      end module io_base