subroutine read_control
C
C Read molecular data
C
      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'
      include "COMMON.SPLITELE"
      character*320 controlcard,ucase
#ifdef MPL
      include 'COMMON.INFO'
#endif
      integer i

      read (INP,'(a80)') titel
      call card_concat(controlcard)

      call readi(controlcard,'NRES',nres,0)
      call readi(controlcard,'RESCALE',rescale_mode,2)
      call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
      write (iout,*) "DISTCHAINMAX",distchainmax
C Reading the dimensions of box in x,y,z coordinates
      call reada(controlcard,'BOXX',boxxsize,100.0d0)
      call reada(controlcard,'BOXY',boxysize,100.0d0)
      call reada(controlcard,'BOXZ',boxzsize,100.0d0)
c Cutoff range for interactions
      call reada(controlcard,"R_CUT",r_cut,15.0d0)
      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
      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
      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
c--------------------------------------------------------------------------
      subroutine molread
C
C Read molecular data.
C
      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*4 sequence(maxres)
      character*800 weightcard
      integer rescode
      double precision x(maxvar)
      integer itype_pdb(maxres)
      logical seq_comp
      integer i,j,kkk
C
C Body
C
C Read weights of the subsequent energy terms.
      call card_concat(weightcard)
      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)
      if (index(weightcard,'SOFT').gt.0) ipot=6
C 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
      weights(19)=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

C 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
C 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
c      if (pdbref) then
c        read(inp,'(a)') pdbfile
c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
c        open(ipdbin,file=pdbfile,status='old',err=33)
c        goto 34 
c  33    write (iout,'(a)') 'Error opening PDB file.'
c        stop
c  34    continue
c        print *,'Begin reading pdb data'
c        call readpdb
c        print *,'Finished reading pdb data'
c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
c        do i=1,nres
c          itype_pdb(i)=itype(i)
c        enddo
c        close (ipdbin)
c        write (iout,'(a,i3)') 'nsup=',nsup
c        nstart_seq=nnt
c        if (nsup.le.(nct-nnt+1)) then
c          do i=0,nct-nnt+1-nsup
c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
c              nstart_seq=nnt+i
c              goto 111
c            endif
c          enddo
c          write (iout,'(a)') 
c     &            'Error - sequences to be superposed do not match.'
c          stop
c        else
c          do i=0,nsup-(nct-nnt+1)
c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
c     &      then
c              nstart_sup=nstart_sup+i
c              nsup=nct-nnt+1
c              goto 111
c            endif
c          enddo 
c          write (iout,'(a)') 
c     &            'Error - sequences to be superposed do not match.'
c        endif
c  111   continue
c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
c     &                 ' nstart_seq=',nstart_seq
c      endif
      call init_int_table
      call 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
c-----------------------------------------------------------------------------
      logical function seq_comp(itypea,itypeb,length)
      implicit none
      integer length,itypea(length),itypeb(length)
      integer i
      do i=1,length
        if (itypea(i).ne.itypeb(i)) then
          seq_comp=.false.
          return
        endif
      enddo
      seq_comp=.true.
      return
      end
c-----------------------------------------------------------------------------
      subroutine read_bridge
C Read information about disulfide bridges.
      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.TIME1'
#ifdef MPL
      include 'COMMON.INFO'
#endif
      integer i,j
C Read bridging residues.
      read (inp,*) ns,(iss(i),i=1,ns)
      print *,'ns=',ns
C Check whether the specified bridging residues are cystines.
      do i=1,ns
	if (itype(iss(i)).ne.1) then
	  write (iout,'(2a,i3,a)') 
     &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
     &   ' can form a disulfide bridge?!!!'
	  write (*,'(2a,i3,a)') 
     &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
     &   ' can form a disulfide bridge?!!!'
#ifdef MPL
	 call mp_stopall(error_msg)
#else
         stop
#endif
        endif
      enddo
C Read preformed bridges.
      if (ns.gt.0) then
      read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
      if (nss.gt.0) then
        nhpb=nss
C Check if the residues involved in bridges are in the specified list of
C 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 MPL
	      call mp_stopall(error_msg)
#else
              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)=dbr
          forcon(i)=fbr
        enddo
        do i=1,nss
          ihpb(i)=ihpb(i)+nres
          jhpb(i)=jhpb(i)+nres
        enddo
      endif
      endif
      return
      end
c----------------------------------------------------------------------------
      subroutine read_angles(kanal,*)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      integer i,kanal
      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
        theta(i)=deg2rad*theta(i)
        phi(i)=deg2rad*phi(i)
        alph(i)=deg2rad*alph(i)
        omeg(i)=deg2rad*omeg(i)
      enddo
      return
   10 return1
      end
c----------------------------------------------------------------------------
      subroutine reada(rekord,lancuch,wartosc,default)
      implicit none
      character*(*) rekord,lancuch
      double precision wartosc,default
      integer ilen,iread
      external ilen
      iread=index(rekord,lancuch)
      if (iread.eq.0) then
        wartosc=default 
        return
      endif   
      iread=iread+ilen(lancuch)+1
      read (rekord(iread:),*) wartosc
      return
      end
c----------------------------------------------------------------------------
      subroutine multreada(rekord,lancuch,tablica,dim,default)
      implicit none
      integer dim,i
      double precision tablica(dim),default
      character*(*) rekord,lancuch
      integer ilen,iread
      external ilen
      do i=1,dim
        tablica(i)=default 
      enddo
      iread=index(rekord,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
c----------------------------------------------------------------------------
      subroutine readi(rekord,lancuch,wartosc,default)
      implicit none
      character*(*) rekord,lancuch
      integer wartosc,default
      integer ilen,iread
      external ilen
      iread=index(rekord,lancuch)
      if (iread.eq.0) then
        wartosc=default 
        return
      endif   
      iread=iread+ilen(lancuch)+1
      read (rekord(iread:),*) wartosc
      return
      end
c----------------------------------------------------------------------------
      subroutine card_concat(card)
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      character*(*) card
      character*80 karta,ucase
      external ilen
      read (inp,'(a)') karta
      karta=ucase(karta)
      card=' '
      do while (karta(80:80).eq.'&')
        card=card(:ilen(card)+1)//karta(:79)
        read (inp,'(a)') karta
        karta=ucase(karta)
      enddo
      card=card(:ilen(card)+1)//karta
      return
      end
c----------------------------------------------------------------------------
      subroutine openunits
      implicit none
      include 'DIMENSIONS'    
#ifdef MPI
      include "mpif.h"
      character*3 liczba
      include "COMMON.MPI"
#endif
      include 'COMMON.IOUNITS'
      include 'COMMON.CONTROL'
      integer lenpre,lenpot,ilen
      external ilen
      character*16 cformat,cprint
      character*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)
C 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')
C 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
C
C 8/9/01 In the newest version SCp interaction constants are read from a file
C Use -DOLDSCP to use hard-coded constants instead.
C
      call getenv('SCPPAR',scpname)
      open (iscpp,file=scpname,status='old')
#endif
      return
      end