subroutine read_ref_structure(*)
C
C Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
C angles.
C
      implicit none
      include 'DIMENSIONS'
      include 'sizesclu.dat'
      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'
      character*4 sequence(maxres)
      integer rescode
      double precision x(maxvar)
      integer itype_pdb(maxres)
      logical seq_comp
      integer i,j,k,nres_pdb,iaux
      double precision ddsc,dist
      integer ilen
      external ilen
C
      nres0=nres
      write (iout,*) "pdbref",pdbref
      if (pdbref) then
        read(inp,'(a)') pdbfile
        write (iout,'(2a,1h.)') 'PDB data will be read from file ',
     &    pdbfile(:ilen(pdbfile))
        call flush(iout)
        open(ipdbin,file=pdbfile,status='old',err=33)
        goto 34 
  33    write (iout,'(a)') 'Error opening PDB file.'
        return1
  34    continue
        do i=1,nres
          itype_pdb(i)=itype(i)
        enddo
        call readpdb(.true.)
        do i=1,nres
          iaux=itype_pdb(i)
          itype_pdb(i)=itype(i)
          itype(i)=iaux
        enddo
        close (ipdbin)
        nres_pdb=nres
        nres=nres0
        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
              do j=nnt+nsup-1,nnt,-1
                do k=1,3
                  cref(k,nres+j+i,1)=cref(k,nres_pdb+j,1)
                enddo
              enddo
              do j=nnt+nsup-1,nnt,-1
                do k=1,3
                  cref(k,j+i,1)=cref(k,j,1)
                enddo
                phi_ref(j+i)=phi_ref(j)
                theta_ref(j+i)=theta_ref(j)
                alph_ref(j+i)=alph_ref(j)
                omeg_ref(j+i)=omeg_ref(j)
              enddo
              nstart_seq=nnt+i
              nstart_sup=nnt+i
              goto 111
            endif
          enddo
#define DEBUG
#ifdef DEBUG
          do j=nnt,nct
            write (iout,'(i5,3f10.5,5x,3f10.5)') 
     &         j,(cref(k,j,1),k=1,3),(cref(k,j+nres,1),k=1,3)
          enddo
          call flush(iout)
#endif
#undef DEBUG
          write (iout,'(a)') 
     &            'Error - sequences to be superposed do not match.'
          return1
        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,'(a,i5)') 
     &   'Experimental structure begins at residue',nstart_seq
#define DEBUG
#ifdef DEBUG
        do j=nnt,nct
          write (iout,'(i5,3f10.5,5x,3f10.5)') 
     &        j,(cref(k,j,1),k=1,3),(cref(k,j+nres,1),k=1,3)
        enddo
        call flush(iout)
#endif
#undef DEBUG
      else
        call read_angles(inp,*38)
        goto 39
   38   write (iout,'(a)') 'Error reading reference structure.'
        return1
   39   call chainbuild     
        nstart_sup=nnt
        nstart_seq=nnt
        nsup=nct-nnt+1
        do i=1,2*nres
          do j=1,3
            cref(j,i,1)=c(j,i)
          enddo
        enddo
      endif
      nend_sup=nstart_sup+nsup-1
      do i=1,2*nres
        do j=1,3
          c(j,i)=cref(j,i,1)
        enddo
      enddo
      do i=1,nres
        do j=1,3
          dc(j,nres+i)=cref(j,nres+i,1)-cref(j,i,1)
        enddo
        if (itype(i).ne.10) then
          ddsc = dist(i,nres+i)
          do j=1,3
            dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
          enddo
        else
          do j=1,3
            dc_norm(j,nres+i)=0.0d0
          enddo
        endif
c        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
c         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
c         dc_norm(3,nres+i)**2
        do j=1,3
          dc(j,i)=c(j,i+1)-c(j,i)
        enddo
        ddsc = dist(i,i+1)
        do j=1,3
          dc_norm(j,i)=dc(j,i)/ddsc
        enddo
      enddo
      write (iout,'(a,i3,a,i3,a,i3,a)')
     &    'Number of residues to be superposed:',nsup,
     &    ' (from residue',nstart_sup,' to residue',
     &    nend_sup,').'
      call flush(iout)
      return
      end