subroutine setsys
        implicit none
        include 'pimc.par'
        integer il,nb,l,js,jsp,i,is,mslice,k,ik,mstep,mmav
        integer lm,jm,jk,mmm,mmmm,n,m,nj,mp,inf,js1,is1,it,kp
        integer mb,nxtalp,ncell(mdim),nnd,nsew,nshadd,nvadd
        integer iaddks,kk,ks,mult,kcomp(mdim)
        integer nvacan
        real*8 xin(mparts),rmin,rmax,epsinf,eta,epsbarinv,aicutk
        real*8 wsites(mdim,mparts),displ,rronf,ronf,pll(mdim)
        real*8 sang,cutkew,omegaau,tempau,rsau,chargeau,erfc
        real*8 rrho,rrs,rksq,rtmin,rtmax,gamma,ellmin
        external sang
        logical isex
        include 'pimc.cm'
        
c    ** Input run parameters **	
        write (*,*) 'name of the run'
        read (*,'(a14)') filen
        ln=index(filen,' ')-1
c does this file exist?
        inquire(file=filen(1:ln)//'.in',exist=ifex)
        if(ifex) then
         open(1,file=filen(1:ln)//'.in',status='old',form='formatted')
         open(2,file=filen(1:ln)//'.out',STATUS= 'UNKNOWN',
     &               form='formatted')
         open(4,file=filen(1:ln)//'.dat',STATUS= 'UNKNOWN',
     &               form='formatted')
         open(10,file=filen(1:ln)//'.eff',STATUS= 'UNKNOWN',
     &         form='formatted')
        else 
         stop ' file does not exist '
        endif
        read(1,'(a80)') title
        write(*,'(a80)') title
        write(2,'(a80)') title

        read(1,'(a80)') title
        read(1,*) ndim,amass,charge
        write(*,'(a80)') title
        write(*,'(i5,2f15.6)') ndim,amass,charge
        write(2,'(a80)') title
        write(2,'(i5,2f15.6)') ndim,amass,charge

        read(1,'(a80)') title
        read(1,*) temp,bond,zp
        write(*,'(a80)') title
        write(*,'(3f15.6)') temp,bond,zp
        write(2,'(a80)') title
        write(2,'(3f15.6)') temp,bond,zp

        read(1,'(a80)') title
        read(1,*) np,nslices
        write(*,'(a80)') title
        write(*,'(6i8,f10.5)') 
     &       np,nslices
        write(2,'(a80)') title
        write(2,'(6i8,f10.5)') 
     &       np,nslices

        read(1,'(a80)') title
        write(*,'(a80)') title
        write(2,'(a80)') title

        read(1,'(a80)') title
        read(1,*) nblock,nstep,nstart,nmovers,nanal
        write(*,'(a80)') title
        write(*,'(5i8)') nblock,nstep,nstart,nmovers,nanal
        write(2,'(a80)') title
        write(2,'(5i8)') nblock,nstep,nstart,nmovers,nanal

        read(1,'(a80)') title
        read(1,*) nalg,ntrial,hop,gammadsp
        write(*,'(a80)') title
        write(*,'(2i8,2f10.5)') nalg,ntrial,hop,gammadsp
        write(2,'(a80)') title
        write(2,'(2i8,2f10.5)') nalg,ntrial,hop,gammadsp

        read(1,'(a80)') title
        read(1,*) iord,epsph,nspill
        write(*,'(a80)') title
        write(*,'(i7,g18.6,i8)') iord,epsph,nspill
        write(2,'(a80)') title
        write(2,'(i5,g18.6,i8)') iord,epsph,nspill

        read(1,'(a80)') title
        read(1,*) idbg!,itabu
        write(*,'(a80)') title
        write(*,'(i7,2g12.6,2i5)') idbg!,itabu
        write(2,'(a80)') title
        write(2,'(i5,2g12.6,2i5)') idbg!,itabu

        if (np.gt.mparts) stop'np.gt.mparts'
        if (nslices.gt.mslices) stop 'nslices.gt.mslices'
        if (ndim.ne.mdim) stop 'ndim.ne.mdim'
        if(mod(nslices,2).ne.0)
     &   pause'odd number of time slices not allowed'
        if (iord.gt.iomax) stop 'iord.gt.iomax'
        if (np.gt.1.and.ncumtr.gt.0) 
     &      stop 'np>1 & ncumtr>0 not implemented'
        if(np.ge.1) then
         open(20,file=filen(1:ln)//'.gr',STATUS= 'UNKNOWN',
     &                 form='formatted')
         if(ifewald.ne.0)
     &    open(21,file=filen(1:ln)//'.sk',STATUS= 'UNKNOWN',
     &    form='formatted')
        endif
        if (nspill.ge.0.and.nspill.le.nstep) then
         open(31,file=filen(1:ln)//'.pc',form='unformatted',
     &           status='new')
        endif
        pi=4.d0*datan(1.d0)
        spi=dsqrt(pi)
        sq2=dsqrt(2.d0)
        sq2pi=sq2*spi
c atomic units: m_electr=1, a_0=hbar2/e2/me=1 (Bohr radius)
c 
c change unit system from atomic (hydrogen) to polaronic units
c 
        beta = 1.d0 / temp
        tau=beta/nslices
c alamb=1/(2 )    
        alamb=0.5d0
        sigma=dsqrt(alamb*tau)
        sigma2=dsqrt(2.0d0*alamb*tau)
        cke=0.25d0/(alamb*tau)

        costact=0.d0
        call gauleg(-0.5d0,0.5d0,xnode,wg,iord)
        io2=(iord+1)/2
        io2nl=iord*(iord-1)/2-(iord+1)/2+iord+1
        iosg=(((iord+1)/2)**2-(iord+1)/2)/2+(iord+1)/2
        do i=1,iord
         up2(i)=(.5d0+xnode(i))
         um2(i)=(.5d0-xnode(i))
        enddo
        write (*,*)
        write (*,*) 'Gaussian integration of order ',iord
        write (2,*)
        write (2,*) 'Gaussian integration of order ',iord
c set the beginning of the random sequence
        inquire(file='rmaseed.dat',exist=ifex)
        if (ifex) then
         call rmaset(1)
        else
         call rmaset(0)
        endif

c    ** Initial position **
        nin=0
        nout=nblock
        inquire(file=filen(1:ln)//'.rs',exist=ifex)
        if(ifex) then
         open(3,file=filen(1:ln)//'.rs',STATUS= 'UNKNOWN',
     &              form='formatted')
         read(3,*) nnd,(ell(l),l=1,nnd),rnn
         read(3,*) mp,mslice,(((x(l,nj,js),l=1,nnd),nj=1,mp)
     &                       ,js=1,mslice)
     $                       ,((proton(l,nj),l=1,nnd),nj=1,2)
         write (*,*) 'positions read in qid.rs'
         write (2,*) 'positions read in qid.rs'
         if (mp.ne.np) stop'mp.ne.np'
         if (mslice.ne.nslices) stop'mslice.ne.nslices'
         if (nnd.ne.ndim) stop'nnd.ne.ndim'
         if (nstart.gt.0) then
          read(3,*) mb,mstep,mmav,(av(i),anorm(i),i=1,mmav)
          if (mstep.ne.nstep) stop 'mstep.ne.nstep'
          if (mmav.gt.mav) stop 'mmav.ne.mav'
          nin=mb
          nout=nout+mb
          write(*,*) 'nstart>0 : averages read in qid.rs'
          write(*,*) 'nin, nout :',nin,nout
          write(2,*) 'nstart>0 : averages read in qid.rs'
          write(2,*) 'nin, nout :',nin,nout,anorm(1)
         endif
        else
         write (*,*) ' rs file not found '
         write (*,*) 'start from random configuration (sites)'
         do l=1,ndim
          proton(l,1)=0.d0
          proton(l,2)=0.d0
          if(l.eq.1) then
           proton(1,1)=-0.5*bond
           proton(1,2)=+0.5*bond
          endif
          do n=1,np
           do js=1,nslices
            call rmar(ronf)
            x(l,n,js) = proton(l,n)
     &                + (2.0d0*ronf-1.0d0)*sigma
           enddo
          enddo
         enddo
c cutoff radial tables at 1/2 the smallest box dimension
        endif
        if (nspill.ge.0.and.nspill.le.nstep) 
     &  write(31) np,nslices,ndim,nstep*nblock/nspill
c sets rptable for tabulation of the potential 
        rtmin=0.05d0
        cutr=3.d0
        rptable(1)=rtmin
        rptable(2)=cutr
        rptable(3)=(lptable-1.)/(cutr-rtmin)
        rptable(4)=1.-rtmin*rptable(3)
        rptable(5)=1.d0/rptable(3)
        rptable(6)=rptable(5)**2/6.d0
        dpi=real(lptable)/cutr
        dr=1.d0/dpi
        level=1

        if (nstart.eq.-1) then
         inquire(file=filen(1:ln)//'.pc',exist=ifex)
         if (ifex) then
          open(31,file=filen(1:ln)//'.pc',form='unformatted',
     &         status='old')
          read(31) mp,mslice,mmm,mmmm,(pll(l),l=1,ndim)
          if (mp.ne.np) stop'mp.ne.np in pc file'
          if (mslice.ne.nslices) stop'mslice.ne.nslices in pc file'
          if (mmm.ne.ndim) stop'mmm.ne.ndim in pc file'
          if (mmmm.lt.nstep*nblock) then
           write(*,*) 'number of confs in pc file:',mmmm
           stop 'mmmm.lt.nstep*nblock in pc file'
          endif
          do l=1,ndim
           if (ell(l).ne.pll(l)) stop'ell(l).ne.pll(l)'
          enddo
          read(1,'(a80)') title
          read(1,'(a80)') title
          read(1,*) epsph,level
          write (*,*) 'imaginary time cut off has been changed'
          write(*,'(a80)') title
          write(*,'(g15.6,i5)') epsph,level
          write (2,*) 'imaginary time cut off has been changed'
          write(2,'(a80)') title
          write(2,'(g15.6,i5)') epsph,level
          close(1)
         else
          stop 'pc file not found'
         endif
        endif
        return
        end