def binding_energy(proton_seq, affinities):
    """Calculate binding energy from proton affinities

    Parameters
    ----------
    proton_seq : ndarray
        protonation state of residues
    affinities : ndarray
        proton affinities for residues

    Returns
    -------
    binding energy : float
        Binding energy in Joules
    """
    return affinities[proton_seq.nonzero()[0]].sum()