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 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      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.CONTACTS1'
      include 'COMMON.PEPTCONT'
      include 'COMMON.TIME1'
      include 'COMMON.COMPAR'
      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,kkk
      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))
        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)
        do kkk=1,nperm
        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,kkk)=cref(k,nres_pdb+j,kkk)
                enddo
              enddo
              do j=nnt+nsup-1,nnt,-1
                do k=1,3
                  cref(k,j+i,kkk)=cref(k,j,kkk)
                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
#ifdef DEBUG
              do j=nnt,nct
                write (iout,'(i5,3f10.5,5x,3f10.5)') 
     &            j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
              enddo
#endif
              nstart_seq=nnt+i
              nstart_sup=nnt+i
              goto 111
            endif
          enddo
          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
        enddo
  111   continue
        write (iout,'(a,i5)') 
     &   'Experimental structure begins at residue',nstart_seq
      else
        call read_angles(inp,*38)
        goto 39
   38   write (iout,'(a)') 'Error reading reference structure.'
        return1
   39   call chainbuild 
        kkk=1    
        nstart_sup=nnt
        nstart_seq=nnt
        nsup=nct-nnt+1
        do i=1,2*nres
          do j=1,3
            cref(j,i,kkk)=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,kkk)
        enddo
      enddo
      do i=1,nres
        do j=1,3
          dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
        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
c      print *,"Calling contact"
      call contact(.true.,ncont_ref,icont_ref(1,1),
     &  nstart_sup,nend_sup)
c      print *,"Calling elecont"
      call elecont(.true.,ncont_pept_ref,
     &   icont_pept_ref(1,1),
     &   nstart_sup,nend_sup)
       write (iout,'(a,i3,a,i3,a,i3,a)')
     &    'Number of residues to be superposed:',nsup,
     &    ' (from residue',nstart_sup,' to residue',
     &    nend_sup,').'
      return
      end