read_pdb Subroutine

subroutine read_pdb(self, unit, error)

   !> Instance of the molecular structure data
   type(structure_type),intent(out) :: self

   !> File handle
   integer,intent(in) :: unit

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   integer, parameter :: p_initial_size = 1000 ! this is going to be a protein

   integer :: iatom, jatom, iresidue, try, stat, atom_type, pos, lnum
   real(wp) :: occ, temp, coords(3)
   real(wp), allocatable :: xyz(:,:)
   type(token_type) :: token
   character(len=4) :: a_charge
   character(len=:), allocatable :: line
   character(len=symbol_length), allocatable :: sym(:)
   type(pdb_data), allocatable :: pdb(:)

   allocate(sym(p_initial_size), source=repeat(' ', symbol_length))
   allocate(xyz(3, p_initial_size), source=0.0_wp)
   allocate(pdb(p_initial_size), source=pdb_data())

   iatom = 0
   iresidue = 0

   stat = 0
   do while(stat == 0)
      call next_line(unit, line, pos, lnum, stat)
      if (index(line, 'END') == 1) exit
      if (index(line, 'ATOM') == 1 .or. index(line, 'HETATM') == 1) then
         if (iatom >= size(xyz, 2)) call resize(xyz)
         if (iatom >= size(sym)) call resize(sym)
         if (iatom >= size(pdb)) call resize(pdb)
         iatom = iatom + 1
         pdb(iatom)%het = index(line, 'HETATM') == 1

         if (len(line) >= 78) then
            ! a4: 13:16, a1: 17:17, a3: 18:20, a1: 22:22
            ! a1: 27:27, a4: 73:76, a2: 77:78, a2: 79:80
            pdb(iatom)%name = line(13:16)
            pdb(iatom)%loc = line(17:17)
            pdb(iatom)%residue = line(18:20)
            pdb(iatom)%chains = line(22:22)
            pdb(iatom)%code = line(27:27)
            pdb(iatom)%segid = line(72:74)
            sym(iatom) = line(77:78)
         else
            token = token_type(len(line)+1, len(line)+1)
            call io_error(error, "Too few entries provided in record", &
               & line, token, filename(unit), lnum, "record too short")
            return
         end if
         if (len(line) >= 80) then
            a_charge = line(79:80)
         else
            a_charge = ""
         end if
         if (stat == 0) then
            ! i5: 7-11
            token = token_type(7, 11)
            call read_token(line, token, jatom, stat)
         end if
         if (stat == 0) then
            ! i4: 23-26
            token = token_type(23, 26)
            call read_token(line, token, pdb(iatom)%residue_number, stat)
         end if
         if (stat == 0) then
            ! f8: 31-38
            token = token_type(31, 38)
            call read_token(line, token, coords(1), stat)
         end if
         if (stat == 0) then
            ! f8: 39-46
            token = token_type(39, 46)
            call read_token(line, token, coords(2), stat)
         end if
         if (stat == 0) then
            ! f8: 47-54
            token = token_type(47, 54)
            call read_token(line, token, coords(3), stat)
         end if
         if (stat == 0) then
            ! f6: 55-60
            token = token_type(55, 60)
            call read_token(line, token, occ, stat)
         end if
         if (stat == 0) then
            ! f6: 61-66
            token = token_type(60, 66)
            call read_token(line, token, temp, stat)
         end if
         if (stat /= 0) then
            call io_error(error, "Cannot read coordinates from record", &
               & line, token, filename(unit), lnum, "unexpected value")
            return
         end if

         xyz(:,iatom) = coords * aatoau
         atom_type = to_number(sym(iatom))
         if (atom_type == 0) then
            try = scan(pdb(iatom)%name, 'HCNOSPF')
            if (try > 0) sym(iatom) = pdb(iatom)%name(try:try)//' '
            pdb(iatom)%charge = 0
         else
            read(a_charge(1:1), *, iostat=stat) pdb(iatom)%charge
            if (stat /= 0) then
               stat = 0
               pdb(iatom)%charge = 0
            else
               if (a_charge(2:2) == '-') pdb(iatom)%charge = -pdb(iatom)%charge
            end if
         end if
         if (to_number(sym(iatom)) == 0) then
            call io_error(error, "Cannot map symbol to atomic number", &
               & line, token_type(77, 78), filename(unit), lnum, "unknown element")
            return
         end if
      end if
   end do

   call new(self, sym(:iatom), xyz(:, :iatom))
   self%pdb = pdb(:iatom)
   self%charge = sum(pdb(:iatom)%charge)

end subroutine read_pdb