def get_stem_energies(stems_s1, stems_s2, HB_matrix, G_HB = -1.89, G_stack = -4.5):
    """ determine the folding energies for each stem:
        stems_s1: list of containing a list of residue indices associated with the 5' residues in each stem
        stems_s2: list of containing a list of residue indices associated with the 3' residues in each stem
        HB_matrix: the N x N matrix with counts of the hydrogen-bonds possible between pairs of residues
        G_HB: the energy contribution per hydrogen (default: -1.89 kcal/mol)
        G_stack: the energy contribution per stacking interaction (default: -4.5 kcal/mol)
    """
    stem_energies = []
    for index_a in range(0, len(stems_s1)):
        stem_1 = stems_s1[index_a]
        stem_2 = stems_s2[index_a]        
        # initialize energy
        ene = 0.0        
        # add base-pair energies
        for index_b in range(0, len(stem_1)):
            ene += HB_matrix[stem_1[index_b], stem_2[index_b]]*G_HB    
        # add stacking
        ene += (float(len(stem_1))-1)*G_stack
        stem_energies.append(ene)
    return(stem_energies)