Calculate a variety of score and distance metrics for the given structures.

def read_and_calculate(workspace, pdb_paths):
    """
    Calculate a variety of score and distance metrics for the given structures.
    """

    # Parse the given restraints file.  The restraints definitions are used to
    # calculate the "restraint_dist" metric, which reflects how well each
    # structure achieves the desired geometry. Note that this is calculated
    # whether or not restraints were used to create the structures in question.
    # For example, the validation runs don't use restraints but the restraint
    # distance is a very important metric for deciding which designs worked.

    restraints = parse_restraints(workspace.restraints_path)

    # Calculate score and distance metrics for each structure.

    from klab.bio.basics import residue_type_3to1_map

    records = []
    metadata = {}
    num_restraints = len(restraints) + 1
    atom_xyzs = {}
    fragment_size = 0

    # It's kinda hard to tell which lines are part of the score table.  The 
    # first column has some pretty heterogeneous strings (examples below) and 
    # all the other columns are just numbers.  My strategy here is to try to 
    # make a regular expression that matches all of these examples, with the 
    # exception of the ligand.  I think the ligand will simply be too 
    # heterogeneous to match, and the purpose of this is to get dunbrack 
    # scores, which the ligand doesn't have.
    #
    #   MET:NtermProteinFull_1
    #   ASN_2
    #   LYS:protein_cutpoint_lower_39
    #   ASP:protein_cutpoint_upper_40
    #   ALA:CtermProteinFull_124
    #   HIS_D_224
    #   pdb_EQU_250

    score_table_pattern = re.compile(
            r'^[A-Z]{3}(?:_[A-Z])?'  # Residue name with optional tautomer.
            r'(?::[A-Za-z_]+)?'      # Optional patch type.
            r'_([0-9]+) '            # Residue number preceded by underscore.
    )                                # The terminal space is important to match
                                     # the full residue number.

    for i, path in enumerate(sorted(pdb_paths)):
        record = {'path': os.path.basename(path)}
        sequence = ""
        sequence_map = {}
        last_residue_id = None
        dunbrack_index = None
        dunbrack_scores = {}

        # Update the user on our progress, because this is often slow.

        sys.stdout.write("\rReading '{}' [{}/{}]".format(
            os.path.relpath(os.path.dirname(path)), i+1, len(pdb_paths)))
        sys.stdout.flush()

        # Read the PDB file, which we are assuming is gzipped.

        try:
            with gzip.open(path) as file:
                lines = file.readlines()
        except IOError:
            print "\nFailed to read '{}'".format(path)
            continue

        if not lines:
            print "\n{} is empty".format(path)
            continue

        # Get different information from different lines in the PDB file.  Some
        # of these lines are specific to different simulations.

        for line in lines:
            line = line.decode('utf8')

            score_table_match = \
                    dunbrack_index and score_table_pattern.match(line)

            if line.startswith('pose'):
                meta = ScoreMetadata(
                        name='total_score',
                        title='Total Score',
                        unit='REU',
                        order=1,
                )
                record['total_score'] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('label'):
                fields = line.split()
                dunbrack_index = fields.index('fa_dun')

            elif score_table_match:
                residue_id = int(score_table_match.group(1))
                for restraint in restraints:
                    if residue_id in restraint.residue_ids:
                        dunbrack_score = float(line.split()[dunbrack_index])
                        dunbrack_scores[residue_id] = dunbrack_score
                        break

            elif line.startswith('rmsd'):
                meta = ScoreMetadata(
                        name='loop_rmsd',
                        title='Loop RMSD (Backbone Heavy-Atom)',
                        unit='Å',
                        guide=1.0, lower=0.0, upper='95%', order=4,
                )
                record[meta.name] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('  all_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats',
                        title='Buried Unsatsified H-Bonds',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('  sc_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats_sidechain',
                        title='Buried Unsatisfied H-Bonds (Sidechain)',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('  bb_heavy_atom_unsats'):
                meta = ScoreMetadata(
                        name='buried_unsats_backbone',
                        title='Buried Unsatisfied H-Bonds (Backbone)',
                        order=5,
                )
                record[meta.name] = float(line.split()[2])
                metadata[meta.name] = meta

            elif line.startswith('time'):
                meta = ScoreMetadata(
                        name='simulation_time',
                        title='Simulation Time',
                        unit='sec',
                        order=5,
                )
                record[meta.name] = float(line.split()[1])
                metadata[meta.name] = meta

            elif line.startswith('FragmentScoreFilter '):
                fragment_size = line.split()[2].split('-')[0]

            elif line.startswith('FSF') or line.startswith('FragmentScoreFilter_metric'):
                splitline = line.split()
                if splitline[1] == 'Max':
                    max_res = 0
                    max_crmsd = 0
                    if splitline[3] == 'res:':
                        max_res = splitline[3]
                        meta = ScoreMetadata(
                                name='max_fragment_crmsd_position',
                                title = 'Max {}-Residue Fragment RMSD \
(C-Alpha) Position'.format(fragment_size),
                                order=7)
                    elif splitline[3] == 'score:':
                        max_crmsd = splitline[3]
                        meta = ScoreMetadata(
                                name='max_fragment_crmsd_score',
                                title = 'Max {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                                order=7)

                elif splitline[1] == 'Min':
                    min_res = 0
                    min_crmsd = 0
                    if splitline[3] == 'res:':
                        min_res = splitline[3]
                        meta = ScoreMetadata(
                                name='min_fragment_crmsd_position',
                                title = 'Min {}-Residue Fragment RMSD \
(C-Alpha) Position'.format(fragment_size),
                                order=8)
                    elif splitline[3] == 'score:':
                        min_crmsd = splitline[3]
                        meta = ScoreMetadata(
                                name='min_fragment_crmsd_score',
                                title = 'Min {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                                order=8)

                elif splitline[1] == 'Avg':
                    meta = ScoreMetadata(
                            name='avg_fragment_crmsd',
                            title='Avg {}-Residue Fragment RMSD \
(C-Alpha)'.format(fragment_size),
                            order=9)
                else:
                    position = splitline[2]
                    crmsd = splitline[4]
                    meta = ScoreMetadata(
                            name='fragment_crmsd_pos_{}'.format(position),
                            title='{}-Residue Fragment RMSD at Res {} \
(C-Alpha)'.format(fragment_size,position),
                            order=6)

                record[meta.name] = float(splitline[4])
                metadata[meta.name] = meta

            elif line.startswith('EXTRA_SCORE'):
                tokens = line[len('EXTRA_SCORE_'):].rsplit(None, 1)
                meta = parse_extra_metric(tokens[0], 5)
                record[meta.name] = float(tokens[1])
                metadata[meta.name] = meta

            elif line.startswith('EXTRA_METRIC'):
                tokens = line[len('EXTRA_METRIC '):].rsplit(None, 1)

                # Ignore the BuriedUnsat filter.  It just reports 911 every 
                # time, and we extract the actual buried unsat information from 
                # some other lines it adds to the PDB.
                if tokens[0] == 'IGNORE':
                    continue
                if tokens[0] == 'Buried Unsatisfied H-Bonds [-|#]':
                    continue

                meta = parse_extra_metric(tokens[0], 5)
                record[meta.name] = float(tokens[1])
                metadata[meta.name] = meta

            elif (line.startswith('ATOM') or line.startswith('HETATM')):
                atom_name = line[12:16].strip()
                residue_id = int(line[22:26].strip())
                residue_name = line[17:20].strip()

                # Keep track of this model's sequence.

                if residue_id != last_residue_id:
                    if line.startswith('ATOM'): 
                        one_letter_code = residue_type_3to1_map.get(residue_name, 'X')
                        sequence += one_letter_code
                        sequence_map[residue_id] = one_letter_code
                        last_residue_id = residue_id
                    elif line.startswith('HETATM'):
                        sequence_map[residue_id] = 'X'
                        last_residue_id = residue_id

                # Save the coordinate for this atom.  This will be used later 
                # to calculate restraint distances.

                atom_xyzs[atom_name, residue_id] = xyz_to_array((
                        line[30:38], line[38:46], line[46:54]))

        # Calculate how well each restraint was satisfied.

        restraint_values = {}
        restraint_values_by_residue = {}
        is_sidechain_restraint = {}
        restraint_units = {
                'dist': 'Å',
                'angle': '°',
        }

        for restraint in restraints:
            d = restraint.distance_from_ideal(atom_xyzs)
            metric = restraint.metric
            backbone_atoms = set(['N', 'C', 'CA', 'O'])
            backbone_restraint = backbone_atoms.issuperset(restraint.atom_names)

            restraint_values.setdefault(metric, []).append(d)
            restraint_values_by_residue.setdefault(metric, {})
            for i in restraint.residue_ids:
                restraint_values_by_residue[metric].setdefault(i, []).append(d)
                is_sidechain_restraint[i] = (not backbone_restraint) \
                        or is_sidechain_restraint.get(i, False)

        for metric, values in restraint_values.items():
            meta = ScoreMetadata(
                    name='restraint_{0}'.format(metric),
                    title='Restraint Satisfaction',
                    unit=restraint_units[metric],
                    guide=1.0, lower=0.0, upper='95%', order=2,
            )
            record[meta.name] = np.max(values)
            metadata[meta.name] = meta

        for metric, values_by_residue in restraint_values_by_residue.items():
            if len(values_by_residue) <= 1:
                continue

            for i in values_by_residue:
                # I want to put the amino acid in these names, because I think 
                # it looks nice, but it causes problems for positions that can 
                # mutate.  So I assume that if a position has a sidechain 
                # restraint, it must not be allowed to mutate.
                aa = sequence_map[i] if is_sidechain_restraint[i] else 'X'
                res = '{0}{1}'.format(aa, i)
                meta = ScoreMetadata(
                        name='restraint_{0}_{1}'.format(metric, res.lower()),
                        title='Restraint Satisfaction for {0}'.format(res),
                        unit=restraint_units[metric],
                        guide=1.0, lower=0.0, upper='95%', order=3,
                )
                record[meta.name] = np.max(values_by_residue[i])
                metadata[meta.name] = meta

        # Finish calculating some records that depend on the whole structure.

        record['sequence'] = sequence
        for i, score in dunbrack_scores.items():
            aa = sequence_map[i] if is_sidechain_restraint[i] else 'X'
            res = '{0}{1}'.format(aa, i)
            meta = ScoreMetadata(
                    name='dunbrack_score_{0}'.format(res.lower()),
                    title='Dunbrack Score for {0}'.format(res),
                    unit='REU',
                    order=5,
            )
            record[meta.name] = score
            metadata[meta.name] = meta

        records.append(record)

    if pdb_paths:
        sys.stdout.write('\n')

    return records, metadata