subroutine read_ref_structure(iprot,*)
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 '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'
      include 'COMMON.COMPAR'
      include 'COMMON.VMCPAR'
      character*4 sequence(maxres)
      integer rescode
      double precision x(maxvar)
      integer itype_pdb(maxres)
      logical seq_comp
      integer i,j,k,ib,iref,iprot,nres_pdb
      double precision ddsc,dist
      integer ilen
      external ilen
C
      write (2,*) "PDBREF ",pdbref
      nres0=nres
      DO IB=1,NBETA(1,iprot)
      if (pdbref) then
        iref=0
        write (iout,'(2a,1h.)') 'PDB data will be read from file ',
     &    pdbfile(ib,iprot)(:ilen(pdbfile(ib,iprot)))
        open(ipdbin,file=pdbfile(ib,iprot),status='old',err=33)
        goto 34 
  33    write (iout,'(a)') 'Error opening PDB file.'
        return1
  34    continue
        nres_pdb=nres
        DO IREF=1,MAXREF
        write (2,*) "Reading reference structure",iref
        call readpdb(print_refstr,iprot)
        write (2,*) "After readpdb"
        call flush(iout)
        if (nres_pdb.eq.0) then
          exit
        endif
C Copy the coordinates to reference coordinates
        do i=1,2*nres
          do j=1,3
            cref(j,i,iref,ib,iprot)=c(j,i)
          enddo
        enddo
        do i=1,nres
          itype_pdb(i)=itype(i)
        enddo
        nres_pdb=nres
        nres=nres0
        nstart_seq=nnt
        if (nsup(iprot).le.(nct-nnt+1)) then
          do i=0,nct-nnt+1-nsup(iprot)
            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup(iprot)),
     &        nsup(iprot))) then
              do j=nnt+nsup(iprot)-1,nnt,-1
                do k=1,3
                  cref(k,j+i,iref,ib,iprot)=cref(k,j,iref,ib,iprot)
                  cref(k,nres+j+i,iref,ib,iprot)=
     &                             cref(k,nres_pdb+j,iref,ib,iprot)
                enddo
c                phi_ref(j+i,iprot)=phi_ref(j,iprot)
c                theta_ref(j+i,iprot)=theta_ref(j,iprot)
c                alph_ref(j+i,iprot)=alph_ref(j,iprot)
c                omeg_ref(j+i,iprot)=omeg_ref(j,iprot)
              enddo
#define DEBUG
#ifdef DEBUG
              do j=nnt,nct
                write (iout,'(i5,3f10.5,5x,3f10.5)') 
     &            j,(cref(k,j,iref,ib,iprot),k=1,3),
     &              (cref(k,j+nres,iref,ib,iprot),k=1,3)
              enddo
#endif
#undef DEBUG
              nstart_seq=nnt+i
              nstart_sup(iprot)=nnt+i
              goto 111
            endif
          enddo
          write (iout,'(a)') 
     &            'Error - sequences to be superposed do not match.'
          return1
        else
          do i=0,nsup(iprot)-(nct-nnt+1)
            if (seq_comp(itype(nnt),itype_pdb(nstart_sup(iprot)+i),
     &        nct-nnt+1)) 
     &      then
              nstart_sup(iprot)=nstart_sup(iprot)+i
              nsup(iprot)=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
        ENDDO ! IREF
        nref(ib,iprot)=iref-1
        write (iout,*) "NREF",ib,iprot,nref(ib,iprot)
        call flush(iout)
      else
        DO IREF=1,MAXREF
        call read_angles(inp,*38)
        goto 39
   38   nref(ib,iprot)=iref-1
        exit
   39   call chainbuild     
        nstart_sup(iprot)=nnt
        nstart_seq=nnt
        nsup(iprot)=nct-nnt+1
        do i=1,2*nres
          do j=1,3
            cref(j,i,iref,ib,iprot)=c(j,i)
          enddo
        enddo
        ENDDO ! IREF
        nref(ib,iprot)=iref-1
      endif
      nend_sup(iprot)=nstart_sup(iprot)+nsup(iprot)-1
c      do i=1,2*nres
c        do j=1,3
c         c(j,i)=cref(j,i,iprot)
c        enddo
c      enddo
      do i=1,nres
c        do j=1,3
c          dc(j,nres+i)=cref(j,nres+i,iprot)-cref(j,i,iprot)
c        enddo
c        if (itype(i).ne.10) then
c          ddsc = dist(i,nres+i)
c          do j=1,3
c            dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
c          enddo
c        else
c          do j=1,3
c            dc_norm(j,nres+i)=0.0d0
c          enddo
c        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"
c      call contact(print_contact,ncont_ref(iprot),icont_ref(1,1,iprot),
c     &  nstart_sup(iprot),nend_sup(iprot))
c      print *,"Calling elecont"
c      call elecont(print_contact,ncont_pept_ref(iprot),
c     &   icont_pept_ref(1,1,iprot),
c     &   nstart_sup(iprot),nend_sup(iprot))
c       write (iout,'(a,i3,a,i3,a,i3,a)')
c     &    'Number of residues to be superposed:',nsup(iprot),
c     &    ' (from residue',nstart_sup(iprot),' to residue',
c     &    nend_sup(iprot),').'
      ENDDO ! IB
      close (ipdbin)
      return
      end